|
|
|
|
||
|
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.
|
||
|
|
|