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

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)?

• What language or software are you using? The slope and aspect metrics are available in just about any GIS or remote sensing software you can name. It is also likely coded in a much more efficient way than you could come up with. So, why? If it is just that you want to understand the algorithm then please edit your question to indicate so. Apr 21 at 14:38
• Thank you very much for your response. I am using Python and as you surmised, I am trying to implement and understand the algorithm, not use a library to calculate for me. Apr 21 at 23:01

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.

• Thanks so much, this is super useful. I notice you provided calculations for the center cell. What would slope for an edge cell (say cell with value 2) and corner cell(cell 1 for example) be for, say, the Zevenberg/Thorne algorithm? are missing neighbours assigned value 0 ? Apr 25 at 19:36