Calculating slope using the 3*3 window approach for elevation data

Question

I have some DEM grid data (25x25, 625 pixels) and I am trying to figure out the correct formula for calculating slope for each pixel. Some googling gave me a bunch of different formulas, but I've figured out that a 3*3 window approach where the value of neighbours of each pixel is used to calculate rise and run is the most common.

From what I can understand, you calculate elevation difference in the x and y directions (taking neighbours into account) and then compute arctan(sqrt(x^2 + y^2)). But I'm trying to figure out the correct formula to calculate elevation difference in the x and y directions. Additionally I'm trying to understand how the formula changes for edge cells which have only 5 neighbours and corner cells which only have 3 neighbours.

As an example, for a 3x3 matrix like

1 2 3
4 5 6
7 8 9

what would slope be for cell 5 (normal cell with 8 neighbours vs cell 2(edge cell with 5 neighbours) and cell 1(corner cell with three neighbours)?

Answer 1

Here is a shotgun answer, detailing four common methods in deriving terrain attributes. Since you do not indicate a language nor software this is provided in a matrix index context. In an application context this would be applied in a kernel (focal) type function.

res = cell resolution of raster

# Evans (1980)
r <- (m[1]+m[3]+m[4]+m[6]+m[7]+m[9]-(2*(m[2]+m[5]+m[8])))/(3*(res^2))
  tx <- (m[1]+m[2]+m[3]+m[7]+m[8]+m[9]-(2*(m[4]+m[5]+m[6])))/(3*(res^2))
    s <- (m[3]+m[7]-m[1]-m[9])/(4*(res^2))
  p <- (m[3]+m[6]+m[9]-m[1]-m[4]-m[7])/(6*res)
q <- (m[1]+m[2]+m[3]-m[7]-m[8]-m[9])/(6*res) 

# Shary (1995)  
r=(m[1]+m[3]+m[7]+m[9]+3*(m[4]+m[6])-2*(m[2]+3*m[5]+m[8]))/(5*(res^2))
  tx=(m[1]+m[3]+m[7]+m[9]+3*(m[2]+m[8])-2*(m[4]+3*m[5]+m[6]))/(5*(res^2))
    s=(m[3]+m[7]-m[1]-m[9])/(4*(res^2))
  p=(m[3]+m[6]+m[9]-m[1]-m[4]-m[7])/(6*res)
q=(m[1]+m[2]+m[3]-m[7]-m[8]-m[9])/(6*res)

# Zevenbergen and Thorne (1987)
p=(m[6]-m[4])/(2*res)
  q=(m[2]-m[8])/(2*res)
    r=(m[4]+m[6]-2*m[5])/(2*(res^2))
  s=(m[3]+m[7]-m[1]-m[9])/(4*(res^2))
tx=(m[2]+m[8]-2*m[5])/(2*(res^2))

# Moore et al. (1993) 
p=(m[6]-m[4])/(2*res)
  q=(m[2]-m[8])/(2*res)
    r=(m[4]+m[6]-2*m[5])/(res^2)
  s=(m[3]+m[7]-m[1]-m[9])/(4*(res^2))
tx=(m[2]+m[8]-2*m[5])/(res^2)

slope = atan(sqrt(p^2+q^2))
aspect = 180-atan2(q,p)+90*(p/abs(p))
planform curvature = -(q^2*r-2*p*q*s+p^2*tx)/((p^2+q^2)*sqrt(1+p^2+q^2))
profile.curvature = -(p^2*r+2*p*q*s+q^2*tx)/((p^2+q^2)*sqrt(1+p^2+q^2)^3) 

References

Evans, I.S. (1980). An Integrated System of Terrain Analysis and Slope Mapping. Zeitschrift für Geomorphologie, Suppl. Bd. 36, 274-295.

Florinsky, I.V. (1998). Accuracy of Local Topographic Variables Derived from Digital Elevation Models. International Journal of Geographical Information Science, 12(1):47-62.

Moore, I.D., P.E. Gessler, G.A. Nielsen, & G.A. Paterson (1993). Soil Attribute Prediction Using Terrain Analysis. Soil Science Society of America Journal, 57:443-452.

Shary, P.A. (1991). The Second Derivative Topographic Method. In: The Geometry of the Earth Surface Structures, edited by Stepanov, I. N., 15-29 (in Russian).

Wilson, J.P., & J.C. Gallant (2000). Terrain Analysis - Principles and Applications. Wiley.

Zevenbergen, L.W. & C.R. Thorne (1987). Quantitative Analysis of Land Surface Topography. Earth Surface Processes and Landforms, 12:47-56.