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?
Best Answer
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
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
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
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.)