surfcurv.sml

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");```