Syntax Highlighing:
comments, key words, predefined symbols, class members & methods, functions & classes
### surfcurve.sml
### Computes curvatures (profile = vertical and plan = horizontal)
### for each DEM raster cell. Profile curvature is positive for
### convex upward surfaces; plan curvature is positive for convex
### downhill surfaces.
### Curvatures are computed for each DEM cell using a 3 x 3 kernel
### by computing a local best-fit quadratic surface for the nine-cell
### neighborhood. The equation of this surface can be written:
### z = a*x^2 + b*y^2 + c*xy + d*x + ee*y + f
### The best-fit surface is not constrained to exactly fit any of the
### DEM values in the neighborhood; it provides a smoothed surface
### that minimizes the impact of elevation errors in the DEM.
### Algorithm modified from:
### Evans, I.S. (1979) An integrated system of terrain analysis and slope mapping.
### Final report on grant DA-ERO-591-73-G0040, University of Durham, England.
### http://www.soi.city.ac.uk/~jwo/phd/04param.php3
### See also Schmidt, J., Evans, I.S., and Brinkmann, J. (2003), Comparison of
### polynomial models for land surface curvature calculation. International
### Journal of Geographical Information Science, v.17, no.8, p. 797-814.
#### Global variables #########################################
numeric scale;
### With scale = 1, units of curvature are radians / meter.
### With scale = 100, units of curvature are radians / 100 m.
numeric a, b, c, d, ee, f; # coefficients of quadratic equation
raster DEM; # input elevation surface raster
raster Profile, Plan; # output curvature rasters
numeric lin, col; # line and column counters
numeric cellX, cellY; # DEM cell sizes in column and line directions
numeric prfdenom, plndenom; # values of denominators of the two curvature expressions
numeric doSmooth; # flag for whether or not to used averaged DEM value
array numeric kernel[3,3];
# Procedure to read 3x3 kernel of DEM values
# surrounding current cell. Notation of cells in kernel is as follows:
# ----------------
# | z1 | z2 | z3 |
# ----------------
# | z4 | z5 | z6 |
# ----------------
# | z7 | z8 | z9 |
# ----------------
proc readkernel (class raster K, numeric line, numeric colm){
numeric z1 = K[line-1, colm-1];
numeric z2 = K[line-1, colm];
numeric z3 = K[line-1, colm+1];
numeric z4 = K[line, colm-1];
numeric z5 = K[line, colm];
numeric z6 = K[line, colm+1];
numeric z7 = K[line+1, colm-1];
numeric z8 = K[line+1, colm];
numeric z9 = K[line+1, colm+1];
}
# Function to check surrounding area of Kernel. If no difference,
# there is no need to do some computations.
# Also used as a safety against divide by zero.
func levelcheck() {
if (z5 != z1) return true;
if (z5 != z2) return true;
if (z5 != z3) return true;
if (z5 != z4) return true;
if (z5 != z6) return true;
if (z5 != z7) return true;
if (z5 != z8) return true;
if (z5 != z9) return true;
return false;
}
# Function to check kernel values for null cells.
# Returns false if any cell is null.
func nullcheck() {
if (IsNull(z1) ) return false;
if (IsNull(z2) ) return false;
if (IsNull(z3) ) return false;
if (IsNull(z4) ) return false;
if (IsNull(z5) ) return false;
if (IsNull(z6) ) return false;
if (IsNull(z7) ) return false;
if (IsNull(z8) ) return false;
if (IsNull(z9) ) return false;
return true;
}
## procedure to compute coefficients for local quadratic surface
## dx and dy are line and column cell sizes
proc coeff(numeric dx, numeric dy) {
a = (z1 + z3 + z4 + z6 + z7 + z9)/(6 * dx^2) - (z2 + z5 + z8)/(3 * dx^2);
b = (z1 + z2 + z3 + z7 + z8 + z9)/(6 * dy^2) - (z4 + z5 + z6)/(3 * dy^2);
c = (z3 + z7 - z1 - z9)/(4 * dx * dy);
d = (z3 + z6 + z9 - z1 - z4 - z7)/(6 * dx);
ee = (z1 + z2 + z3 - z7 - z8 -z9)/(6 * dy);
}
proc makeCurv(class raster D) {
for each D[lin,col] {
if (lin != 1 && lin != D.$Info.NumLins && col != 1 && col != D.$Info.NumCols) { ## don't process edge cells
readkernel(D, lin, col); # read 3x3 kernel values and check for null cells
if ( nullcheck() ) {
if( levelcheck() ) { ## don't process cell if 3x3 kernel is level
coeff(cellX, cellY); ## call procedure to compute coefficients for quadratic surface
# compute denominators for profile and plan curvature expressions to check for division by zero
prfdenom = round(10^7 * ( (sqr(ee) + sqr(d) ) * ( (1 + sqr(d) + sqr(ee) ) ^ 1.5)) ) * 10^-7;
plndenom = round(10^7 * ( (sqr(ee) + sqr(d) )^1.5) ) * 10^-7;
if (prfdenom == 0 || plndenom == 0) { ## avoid division by zero (can occur for summits, depressions, and saddles)
if ( (a > 0 && b > 0) || (a < 0 && b < 0) ) {
Profile[lin, col] = - (a + b) * scale; ## case for summit or depression
Plan[lin, col] = -32768;
}
else {
Profile[lin,col] = 0; ## case for saddles
Plan[lin, col] = 0;
}
} # end if (XXXdenom == 0)
else { # curvature computation for typical cells
Profile[lin, col] = scale * -2 * ( a * sqr(d) + b * sqr(ee) + c * d * ee) / prfdenom;
Plan[lin, col] = scale * -2 * (b * sqr(d) + a * sqr(ee) - c * d * ee) / plndenom;
}
}
else {
Plan[lin, col] = -32768; ## case for flat kernel area
Profile[lin, col] = 0;
}
} # end if ( levelcheck() )
else {
# case for kernels including null cells, not enough data to compute curvature
Plan[lin,col] =-32768;
Profile[lin,col] =-32768;
}
} # end if ( nullcheck() )
else {
# case for edge cells, not enough data to compute curvature
Plan[lin,col] = -32768;
Profile[lin,col] = -32768;
}
}
}
#############################################################
### Main program
GetInputRaster(DEM);
cellX = ColScale(DEM);
cellY = LinScale(DEM);
GetOutputRaster(Profile, DEM.$Info.NumLins, DEM.$Info.NumCols, "32-bit float");
GetOutputRaster(Plan, DEM.$Info.NumLins, DEM.$Info.NumCols, "32-bit float");
SetNull(Profile,-32768);
SetNull(Plan,-32768);
scale = PopupNum("Enter 1 to compute curvature in radians/m or \nenter 100 for curvature in radians/100 m", 1, 100, 0);
doSmooth = PopupYesNo("Use averaged DEM values for curvature?", 1);
if (doSmooth == 1) {
class RASTER Temp;
# pop up dialog window to get size of filter window for averaging DEM values
local numeric filtsize = PopupNum("Size of averaging window in cells? \n(3, 5, 7, or 9)", 5, 3, 9, 0);
if ( filtsize % 2 == 0 ) then ## check for even number entered as window size
filtsize = filtsize + 1;
CreateTempRaster(Temp, DEM.$Info.NumLins, DEM.$Info.NumCols, DEM.$Info.Type);
## compute averaged DEM value for non-null cell, else assign null value to Temp raster
numeric i, j;
for i = 1 to DEM.$Info.NumLins {
for j = 1 to DEM.$Info.NumCols {
if ( IsNull(DEM[i,j]) ) then
Temp[i,j] = -32768;
else
Temp[i,j] = FocalMean(DEM, filtsize, filtsize, i, j);
}
}
SetNull(Temp, -32768);
makeCurv(Temp);
CloseRaster(Temp);
} # end of if (doSmooth() )
else
makeCurv(DEM);
for each Profile {
if (Profile < -32768) then
Profile = -32768;
if (Plan < -32768) then
Plan = -32768;
}
SetNull(Profile, -32768);
SetNull(Plan, -32768);
CreateHistogram(Profile);
CreateHistogram(Plan);
CreatePyramid(Profile);
CreatePyramid(Plan);
CopySubobjects(DEM, Profile, "georef");
CopySubobjects(DEM, Plan, "georef");
Download