6

I am currently developing a program that requires that I calculate the slope for each cell in a DEM. However when I try to understand the algorithm to be used I get stuck trying to understand exactly what is going on. The slope algorithm in question is the one detailed in the article Hillshading and the reflectance map (Horn 1981 p. 20-21).

The algorithm is also explained here

I understand how the average height factors for the x and y direction gets calculated and I also understand that the factors have to be combined if the slope of the surface in an arbitrary direction is to be established. I dread that the part which keeps troubling me is probably the easiest.

The part I do not understand is why the Pythagorean theorem is used When the two factors are combined into a single average value? Both dx and dy will be values rising towards z, so how can the Pythagorean theorem be applied to two lines along the same axis? I've tried to visualize the dx and dy values as three-dimensional vectors with dx being (1, 0, dx) and dy being (0, 1, dy). Adding these two vectors together produces the two sides in a right triangle where the hypotenuse will be equal to the answer the equation is supposed to give me. I am then left with the hypotenuse vector as the rise and that just feels wrong when I attempt to perceive it as the rise. Am I misunderstanding something here?

8

Summary

The slope of a surface is determined by the angle made by its normal (perpendicular) vector to the "straight up" direction. The relevant vectors are the 3D gradient of the surface (thought of as the zero set of a differentiable function) and the "straight up" vector (0,0,1). Angles are determined by dot products, leading to straightforward formulas. (The use of dot products is equivalent to the Pythagorean Theorem, which therefore seems to pop up everywhere.) Because gradients are involved, the DEM is used to provide discrete approximations to the derivatives. Different slope algorithms result from different approximations to the derivatives.

These ideas are easily extended to compute directional slopes, leading to algorithms for all the first-order topographic calculations performed by GISes, including contouring, flow tracking, and much more. Understanding these concepts gives you the ability to find efficient ways of solving many problems of topographic analysis and to create new algorithms where needed.

The conceptual setting

A DEM is a discrete representation of a surface which we may abstractly think of as the graph of a function f. This means f is thought of as associating a unique real number to each point within the raster's extent, not just the centers of its cells. To compute a slope for a DEM cell, we begin (at least conceptually) with this functional model to derive a formula and then hope to use information in the DEM to approximate the formula's results.

A technical requirement

Letting x and y be Cartesian coordinates for the DEM locations, we may think of f as being a function of those two variables. It has a slope if and only if it is differentiable with respect to (x,y): this rules out locations where folds, sharp peaks, cusps, cliffs, and edges appear.

How slopes are defined

The surface represented by f is, strictly speaking, a set of points in three dimensions. They are all of the form (x, y, f(x,y)). Equivalently, these points all have (Cartesian) coordinates (x,y,z) which are solutions to the equation

z - f(x,y) = 0. 

Let's call the left hand side F(x,y,z). Its gradient is a vector with components given by the partial derivatives of F. From the definition of F, these components are (-fx, -fy, 1) where fx is the partial derivative of f with respect to x (conceptually, the slope in the x direction) and fy is the partial derivative of f with respect to y (the slope in the y direction). The gradient is always perpendicular to the surface and the slope is defined to be the angle between the gradient and "straight up", which is parallel to (and in the direction of) the unit vector (0,0,1).

Interpretation

All slope information is contained in this normal vector (-fx, -fy, 1). The third component represents a unit vertical change (of the normal vector, not of the surface!): its "rise". The first two components represent "runs" along the two coordinate axes: a distance of fx backwards along the x-axis and a distance of fy backwards along the y-axis. The Pythagorean Theorem tells us the total run within the XY plane must equal sqrt(fx^2 + fy^2).

When the slope is zero, the normal vector points straight up: it will be parallel to (0,0,1). It has no "run" at all. When the slope is extremely steep--almost a cliff--the normal vector will point almost sideways. To do so, the "run" sqrt(fx^2 + fy^2) must be much larger than the "rise" of 1. Evidently--because we are looking at the normal vector rather than at vectors lying within the surface itself--the usual definition of slope must be reversed: it equals the ratio of the run to the rise for the normal vector. Because we have a unit rise, we can read the slope right off the gradient, obtaining

slope = sqrt(fx^2 + fy^2).

That's the key result. The slope as an angle is the inverse tangent of this quantity (as proclaimed in the ArcGIS documentation, for instance).

Generalization to directional slopes

More generally, the slope in some direction given by a horizontal unit vector (u,v,0) in the XY plane is determined by the complement of the angle it makes to the normal vector (-fx, -fy, 1). The sine of this angle is given by the normalized dot product

sin(directional slope angle) = (u*fx + v*fy) / sqrt(1 + fx^2 + fy^2),

because (u,v,0) * (-fx, -fy, 1) = -(u*fx + v*fy) and the length of (-fx, -fy, 1) is sqrt(1 + fx^2 + fy^2). From this it follows immediately that the maximum slope is in the direction parallel to (fx, fy), the minimum is in the opposite direction (parallel to (-fx, -fy)), and halfway between (parallel to (-fy, fx) in either direction) the slope is always zero. Whence (-fy, fx) always points along a contour line and (fx, fy) always points outwards from a contour line (perpendicular to it) towards higher values of f. Given formulas for the derivatives fx and fy (see the next section), these results provide ways to track along contours or down/up from contours, as is done with contouring algorithms, stream flow and watershed calculations, pathdistance calculations based on slopes encountered along routes, and particle-tracking algorithms.

Computing derivatives

All methods of computing slopes of DEMs ultimately use the key result, typically by estimating the partial derivatives (that is, directional slopes) fx and fy by means of finite differences. For instance, a crude and obvious estimate of fx is to compute the ratio rise:run by looking from the immediate right neighbor of a cell to its immediate left neighbor, subtracting the values of the DEM at those points, and dividing by their distance apart (two cell lengths). Similarly, fy can be estimated by looking from the upper neighbor to the lower neighbor, subtracting the values of the DEM at those points, and dividing by two cell widths.

(ArcGIS uses a more refined estimate. It can be understood either as the slope obtained by approximating f by a biquadratic surface determined by all nine points in a 3 by 3 neighborhood, performing a least-squares fit of a plane to those points, or by using a weighted average of fx (as computed crudely) across the top, middle, and lower horizontal strips and, similarly, a weighted average of fy across the left, middle, and right vertical strips.)

  • 2
    Apparently I had misunderstod most of it. I have read through this now a few times and I understand most of it. It is a shame I can only upvote this once. Truly, thank you for going great lengths to assist me in this matter. – Steinin Aug 15 '13 at 12:41
  • 1
    You can hardly be blamed for misunderstanding the ArcGIS documentation, Steinin. As far as I can tell, most of the technical parts of those help pages are filtered through committees of professional re-writers who seem not to understand adequately much of what they present. As a result, an alarming proportion of those pages are ambiguous, useless, or downright incorrect. I firmly believe the programmers are, by and large, writing correct code, but I would like to suggest that you (and just about anyone else) seek information about algorithms and concepts elsewhere. – whuber Aug 15 '13 at 14:19
  • Just a small question in relation to slope. Horn's method is equivalent with the Sobel Operator from image processing? – Steinin Nov 11 '13 at 13:43
  • Up to the sign and a constant of proportionality given by the cell size, they are the same. – whuber Nov 11 '13 at 15:45
  • What you call 'crude and obvious' estimate of gradient is actually pretty accurate mathematically. It's also identical to the slope obtained from the biquadratic surface interpolating 9 points. – mkadunc Sep 13 '16 at 23:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.