Using R and r.mapcalc (GRASS) to Estimate Mean Topographic Curvature

Submitted by dylan on Tue, 2010-08-03 20:51.

Recently I was re-reading a paper on predictive soil mapping (Park et al, 2001), and considered testing one of their proposed terrain attributes in GRASS. The attribute, originally described by Blaszczynski (1997), is the distance-weighted mean difference in elevation applied to an n-by-n window of cells:

Equation 4 from (Park et al, 2001)

where n is the number of cells within an (odd-number dimension) square window excluding the central cell, z is the elevation at the central cell, z_{i} is the elevation at one of the surrounding cells i, d_{i} is the horizontal distance between the central cell and surrounding cell i. I wasn't able to get a quick answer using r.neighbors or r.mfilter, so I cooked up a simple R function to produce a solution using r.mapcalc. The results are compared with the source DEM below; concave regions are blue-ish, convex regions are red-ish. The magnitude and range are almost identical to mean curvature derived from v.surf.rst, with a Pearson's correlation coefficient of 0.99. I think that it would be of general interest to add functionality to r.neighbors so that it could perform distance-weighted versions of commonly used focal functions.

Elevation surface (left) and resulting mean curvature estimate (right)

References

• Blaszczynski, J. Landform characterization with geographical information systems Photogrammetric Engineering and Remote Sensing, 1997, 63, 183-191
• Park, S.; McSweeney, K. & Lowery, B. Identification of the spatial distribution of soils using a process-based terrain characterization Geoderma, 2001, 103, 249-272

Implementation and example usage in R

## assumes square pixels
## i.e. EWres = NSres
f <- function(rast.in, rast.out, file=paste(rast.out, '.mapcalc', sep=''), n=5, res=1)
{
# the matrix must be an odd number of rows / cols
if(n %% 2 == 0)
stop('n must be an odd number')

# init a matrix to evaluate our filter on
x <- matrix(1:n^2, nrow=n, ncol=n)

# where is the center of the matrix
x.center.offset <- median(c(1:n))

# generate an offset from the center of the matrix
the.offset <- c(rev(1 - 1:x.center.offset), 1:(x.center.offset-1))

# convert into matrix of row and collumn-wise offsets
# according to GRASS's r.mapcalc syntax
row.offset.mat <- matrix(rep(the.offset, n), nrow=n, byrow=FALSE)
col.offset.mat <- matrix(rep(the.offset, n), nrow=n, byrow=TRUE)

# compute linear distance from center
# we are assuming the EW and NS resolution is constant
dist.mat <- sqrt((row.offset.mat * res)^2 + (col.offset.mat * res)^2)

# generate script
e <- paste('((', rast.in, ' - ', rast.in, '[', col.offset.mat, ',' , row.offset.mat, ']) / ', dist.mat, ')', sep='')

# remove reference to center cell
e <- e[-median(1:n^2)]

# generate sum
e.sum <- paste(e, collapse=' + ')

# add final syntax, along with final division by number of cells - 1
e.final <- paste('Cs = ( ', e.sum, ') / ', signif((n^2)-1, 2), '.0', sep='')

# write out to file
cat(e.final, file=file)
}

# example usage
f(rast.in='rtk_elev', rast.out='Cs', n=5, res=1)

Resulting r.mapcalc input for a 5x5 window

Cs = ( ((rtk_elev - rtk_elev[-2,-2]) / 2.82842712474619) + ((rtk_elev - rtk_elev[-2,-1]) / 2.23606797749979) + ((rtk_elev - rtk_elev[-2,0]) / 2) + ((rtk_elev - rtk_elev[-2,1]) / 2.23606797749979) + ((rtk_elev - rtk_elev[-2,2]) / 2.82842712474619) + ((rtk_elev - rtk_elev[-1,-2]) / 2.23606797749979) + ((rtk_elev - rtk_elev[-1,-1]) / 1.4142135623731) + ((rtk_elev - rtk_elev[-1,0]) / 1) + ((rtk_elev - rtk_elev[-1,1]) / 1.4142135623731) + ((rtk_elev - rtk_elev[-1,2]) / 2.23606797749979) + ((rtk_elev - rtk_elev[0,-2]) / 2) + ((rtk_elev - rtk_elev[0,-1]) / 1) + ((rtk_elev - rtk_elev[0,1]) / 1) + ((rtk_elev - rtk_elev[0,2]) / 2) + ((rtk_elev - rtk_elev[1,-2]) / 2.23606797749979) + ((rtk_elev - rtk_elev[1,-1]) / 1.4142135623731) + ((rtk_elev - rtk_elev[1,0]) / 1) + ((rtk_elev - rtk_elev[1,1]) / 1.4142135623731) + ((rtk_elev - rtk_elev[1,2]) / 2.23606797749979) + ((rtk_elev - rtk_elev[2,-2]) / 2.82842712474619) + ((rtk_elev - rtk_elev[2,-1]) / 2.23606797749979) + ((rtk_elev - rtk_elev[2,0]) / 2) + ((rtk_elev - rtk_elev[2,1]) / 2.23606797749979) + ((rtk_elev - rtk_elev[2,2]) / 2.82842712474619)) / 24.0

Attempting Blaszczynski Cs

Attempting Blaszczynski Cs on a different platform. Not having the expertise to replicate in R, do you have the output Cs range using your input raster as shown?

Thanks!

Excellent information ! I'd

Excellent information ! I'd like to know what color ramp are you using to display your mean curvature estimate ?

Color Ramp

This is a built-in color ramp from GRASS. You can apply it to a raster like this:

r.colors myraster color=curvature

That's what I taught. 0%

That's what I taught.
0% black
-0.01 blue
-0.001 aqua
-0.0001 cyan
0 white
0.0001 yellow
0.001 orange
0.01 red
100% magenta
I tried with Spearfish60 elevation.10m. But the curvature I get is from ~ <-33 to >35. Even scaling (using -e) doesn't make a pretty result.
Did I miss something ?

Curvature Values

So, after running this script, you get values ranging from {-33,35}? I suspect that something isn't right. Do your inputs match the assumptions made by the script: 1m resolution, etc.? If the values are derived from some other software, you could adjust your color table using only percentages.