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.
RAMSES
Modules:
A-Z
Function
Layer
QuickRefs:
DM
AuxLib
AuxLibE
SciLib
EasyMW
MW
ISIS
RMSLib