ETHZ_Logo RAMSES_Logo_Right   RAMSES   RAMSES_Logo_Left Systems Ecology  
Start    search button      Modules:   A-Z   Function   Layer        QuickRefs:   DM   AuxLib   AuxLibE   SciLib   EasyMW   MW   ISIS   RMSLib

DEFINITION MODULE GridDSOp;

  (*******************************************************************

    Module  GridDSOp     (Version 1.0)

      Copyright (c) 2001-2006 by Dimitrios Gyalistras, Andreas Fischlin
      and ETH Zurich.

    Purpose   Operations on gridded data sets, plus
              auxilary procedures for distance calculation
              and interpolation between gridpoints.

    Remarks   For general management of gridded data sets
              see module GridDataSets.


    Programming

      o Design
        Dimitrios Gyalistras      27/09/2001

      o Implementation
        Dimitrios Gyalistras      27/09/2001
        Andreas Fischlin          24/05/2002


    ETH Zurich
    Systems Ecology
    CHN E 35.1
    Universitaetstrasse 16
    8092 Zurich
    SWITZERLAND

    URLs:
        <mailto:RAMSES@env.ethz.ch>
        <http://www.sysecol.ethz.ch>
        <http://www.sysecol.ethz.ch/SimSoftware/RAMSES>


    Last revision of definition:  24/05/2002  AF

  *******************************************************************)


  FROM DMFiles IMPORT TextFile;
  IMPORT Errors;

  FROM GridDataSets IMPORT GridDataSet;

  CONST
    notDone = Errors.onlyAnInsert;
    (*
      All procedures exported below return errCode = Errors.allOk if the
      call to the procedure was succesfull, and resCode = notDone if not.
    *)


  (*******************************************************************)
  (*#####   Distance calculation and interpolation procedures   #####*)
  (*******************************************************************)


  PROCEDURE Distance( lon1, lat1: REAL; lon2, lat2: REAL ): REAL;
    (* Returns geographical distance (in km) between two points given in degrees. *)

  PROCEDURE EuclidDistance( x1, y1: REAL; x2, y2: REAL ): REAL;
    (* Returns Euclidian distance between two points given in arbitrary units. *)


  PROCEDURE XGridpointDist( gds    : GridDataSet;
                            xCoord : REAL;
                            yCoord : REAL ): REAL;
    (*
      Returns the gridpoint distance at location "xCoord"/"yCoord" in
      direction of the X coordinate (i.e., in longitudinal direction)
      according to the grid definition of the gridded data set "gds".  If the
      data set is "inDegrees" (see procedure DefineGridDataSet) the procedure
      will return the gridpoint distance in km.  Otherwise it will always
      return the data set's cell size, independent of the location
      considered.
    *)

  PROCEDURE YGridpointDist( gds    : GridDataSet;
                            xCoord : REAL;
                            yCoord : REAL ): REAL;
    (*
      Same as XGridpointDist, but in direction of the Y coordinate
      (i.e. in latitudinal direction).
    *)


  PROCEDURE GetZValue( gds         : GridDataSet;
                       dataMatNr   : INTEGER;
                       xCoord      : REAL;
                       yCoord      : REAL;
                       VAR zValue  : REAL;
                       VAR resCode : INTEGER;
                       VAR errTxt  : ARRAY OF CHAR );
    (*
      The procedure estimates the value of the continuous field Z(x,y) that
      underlies the gridded data set "gds" at location ("xCoord","yCoord").
      The location of interest must be within the sector covered by the
      gridded data set "gds", otherwise resCode will be returned <>
      Errors.allOk.  If the location corresponds to a gridpoint of the data
      set, the corresponding gridpoint value is returned.  Otherwise, the
      field value is estimated by fitting a quadratic surface to the 9
      closest gridpoints surrounding the location.  If no such gridpoints
      exist, or if the gridpoint data are missing for one or several of these
      gridpoints a missing value will be returned, and resCode will be
      returned as Errors.allOk.

      Note, for fitting of the quadratic surface it is assumed that all 9
      surrounding gridpoints are ordered on a regular grid (fixed distance
      between any two gridpoints).  In case that the data grid is NOT defined
      "inDegrees" the gridpoint distance is indeed constant and it is given
      by the data set's "cellSize" attribute.  Otherwise the gridpoint
      distance will be estimated for the location of interest as
      0.5*(XGridpointDist+YGridpointDist).  This approximation is probably
      quite inaccurate if the data set extends to high latitudes since
      XGridpointDist decreases with latitude and becomes zero at the Earth's
      poles.  (In contrast, the YGridpointDist is constant and depends only
      on the data set's cell size).

      Reference for the fitting of the quadratic surface:
        Zevenbergen, L.W. & Thorne, C.R. (1987): Quantitative analysis of
        land surface topography. Earth Surface Processes and Landforms,
        12: 47-56.
    *)


  PROCEDURE GetSurfaceParams( gds          : GridDataSet;
                              dataMatNr    : INTEGER;
                              xCoord       : REAL;
                              yCoord       : REAL;
                              VAR slope    : REAL;
                              VAR aspect   : REAL;
                              VAR profCurv : REAL;
                              VAR planCurv : REAL;
                              VAR resCode  : INTEGER;
                              VAR errTxt   : ARRAY OF CHAR );
    (*
      The procedure estimates the surface parameters of the continuous field
      Z(x,y) that underlies the gridded data set "gds" at location
      ("xCoord","yCoord").  The parameters are obtained by fitting a
      quadratic surface to the 9 closest gridpoints surrounding the location
      (see description of procedure GetZValue above).

      Given that the analyzed data matrix contains elevation data (in m), the
      computed surface parameters are as follows:

      slope     = dZ/dS, where S is the distance in the aspect direction,
                  (in %, i.e. 100*[m]/[m])
      aspect    = direction of maximum slope, 0-360 (in [degrees])
                  (i.e., N = 0[deg], S = 180[deg])
      profCurv  = ddZ/ddS in direction of slope, i.e. profile curvature,
                  (in 1/1000[m])
      planCurv  = ddZ/ddS transverse to the slope, i.e. planform curvature
                  (in 1/1000[m])
    *)


  (*************************************************)
  (*#####   Operations on gridded data sets   #####*)
  (*************************************************)


  PROCEDURE SelectGridDataSetRowsCols( parentGDS   : GridDataSet;
                                       minRow      : INTEGER;
                                       maxRow      : INTEGER;
                                       minCol      : INTEGER;
                                       maxCol      : INTEGER;
                                       VAR newGDS  : GridDataSet;
                                       VAR resCode : INTEGER;
                                       VAR errTxt  : ARRAY OF CHAR );
    (*
      Parameters *Row and *Col can be outside the range of
      rows/columns available in parentGDS. The corresponding
      rows/columns will then contain missing values.
    *)


  PROCEDURE SelectGridDataSetSector( parentGDS   : GridDataSet;
                                     minXCoord   : REAL;
                                     maxXCoord   : REAL;
                                     minYCoord   : REAL;
                                     maxYCoord   : REAL;
                                     VAR newGDS  : GridDataSet;
                                     VAR resCode : INTEGER;
                                     VAR errTxt  : ARRAY OF CHAR );
    (*
      Parameters *XCoord and *Lat can be outside the range of
      longitudes or latitudes available in parentGDS. The
      corresponding rows/columns of newGDS will then contain
      missing values.
    *)


  PROCEDURE SmoothGridDataSet( parentGDS    : GridDataSet;
                               smoothFactor : INTEGER;
                               VAR newGDS   : GridDataSet;
                               VAR resCode  : INTEGER;
                               VAR errTxt   : ARRAY OF CHAR );
    (*
      Smoothing produces a new data set by assigning to each
      gridpoint the average of all gridpoints within a quadratic
      area of size "smoothFactor x smoothFactor" surrounding that
      gridpoint.  "smoothFactor" must be odd and > 1.  The output
      grid has the same dimensions as the input grid.  Smoothing
      starts in the top left of the grid (cell with row=col=
      smoothFactor DIV 2) and stops at the bottom right. A missing
      value is assigned to the output grid if the number of valid
      gridpoint values for averaging is < 0.75*smoothFactor^2.
    *)


  PROCEDURE AggregateGridDataSet( parentGDS    : GridDataSet;
                                  aggregFactor : INTEGER;
                                  VAR newGDS   : GridDataSet;
                                  VAR resCode  : INTEGER;
                                  VAR errTxt   : ARRAY OF CHAR );
    (*
      Aggregation reduces the grid size by taking every
      "aggregFactor" rows and columns averages over
      "aggregFactor x aggregFactor" cells.  "aggregFactor" must
      be odd and > 1.  Aggregation starts in the top left of
      the grid (cell with row=col=aggregFactor DIV 2) and stops
      at the bottom right.  A missing value is assigned to
      the output grid if the number of valid gridpoint values
      for averaging is < 0.75*aggregFactor^2.
    *)


  PROCEDURE RegridDataSet( parentGDS   : GridDataSet;
                           refXCoord   : REAL;
                           refYCoord   : REAL;
                           cellSize    : REAL;
                           VAR newGDS  : GridDataSet;
                           VAR resCode : INTEGER;
                           VAR errTxt  : ARRAY OF CHAR );
    (*
      Produces from data set "gds" a new data set that is
      defined on a grid with reference gridpoint ("refXCoord",
      "refYCoord") and cell size "cellSize".  The reference
      gridpoint is assumed to be in the same units as the grid
      of data set "gds" (i.e., either in degrees longitude/
      latitude, or in arbitrary units).  The gridpoint values
      of the new data set are determined as described in
      procedure GetZValue above.
    *)


  PROCEDURE MakeCHCoordsGridDataSet( gds         : GridDataSet;
                                     refXCoord   : REAL;
                                     refYCoord   : REAL;
                                     cellSize    : REAL;
                                     VAR newGDS  : GridDataSet;
                                     VAR resCode : INTEGER;
                                     VAR errTxt  : ARRAY OF CHAR );
    (*
      Produces from gridded data set "gds" which must be given in degree
      coordinates a new data set which is defined on the Swiss (CH)
      coordinate system.  The grid of the new data set is defined by
      parameters "refXCoord", "refYCoord" and "cellSize".  Gridpoint
      values of the new grid are estimated according to procedure
      GetZValue described above.
    *)


  PROCEDURE MakeDegreeCoordsGridDataSet( gds         : GridDataSet;
                                         refXCoord   : REAL;
                                         refYCoord   : REAL;
                                         cellSize    : REAL;
                                         VAR newGDS  : GridDataSet;
                                         VAR resCode : INTEGER;
                                         VAR errTxt  : ARRAY OF CHAR );
    (*
      Produces from gridded data set "gds" which must be given in Swiss
      (CH) coordinates a new data set in degree coordinates.  The grid of
      the new data set is defined by parameters "refLon", "refLat" and
      "cellSize".  Gridpoint values of the new grid are estimated
      according to procedure GetZValue described above.
    *)


  (**********************************************************)
  (*#####   Operations on pairs of gridded data sets   #####*)
  (**********************************************************)


  PROCEDURE GridsAreIdentical( gds1, gds2: GridDataSet;
                               VAR errTxt: ARRAY OF CHAR ): BOOLEAN;
    (*
      Returns TRUE if the two gridded data sets gds1 and gds2 exist and
      if their grids have the same cell size and overlap.  Note, overlapping
      does not mean that the two gridded data sets cover the same geographical
      sectors - it simply means that all gridpoints of the two (infinite)
      underlying grids coincide.
      In case that FALSE is returned, errTxt will contain an explanation
      why this is the case.
    *)


  PROCEDURE MergeGridDataSets( gds1, gds2             : GridDataSet;
                               allowOverwrWithMisVals : BOOLEAN;
                               VAR newGDS             : GridDataSet;
                               VAR resCode            : INTEGER;
                               VAR errTxt             : ARRAY OF CHAR );
    (*
      Merges the two data sets "gds1" and "gds2" into a new
      data set.  Precondition is that both data sets are
      defined on the same grid.  Merging is done such, that
      "gds2" overwrites all values of "gds1" at common
      gridpoints.  However, if flag "allowOverwrWithMisVals" is
      found false, any missing values from "gds2" are not used
      to overwrite non-missing values from "gds1".
    *)



END GridDSOp.

  Contact RAMSES@env.ethz.ch Last updated: 25-Jul-2011 [Top of page]