Slope Algorithms
Znw
|
Zn
|
Zne
|
Zw
|
Z
|
Ze
|
Zsw
|
Zs
|
Zse
|
The slope is taken as a dZ value divided by a horizontal
distance. The horizontal distance is the data spacing
(east-west, north-south, or diagonal) for the one point and nine
point methods, and twice that distance for the four and eight
point methods. For DEMs like the USGS 10 m and 30, the NS
and EW spacings are the same; they differ for data like the USGS
NED, USGS 1:250K DEM, SRTM, and DTED with geographic spacing.
Aspect, the downhill direction, can be computed as a byproduct of the slope
computation.
Slope and aspects, the magnitude and direction of the vector
tangent to the topographic surface pointing downhill at a point,
have been computed with a multitude of methods; Eyton (1991) and Carter (1992) list a
number of references. The following discussion centers on the
methods implemented in MICRODEM; the references point to clear
expositions of the method, and not necessarily the first use of
the method, many of which are in hard to location publications in the gray
literature. The algorithms fall into several categories,
depending on the number of neighboring points considered to
compute the slope and aspect.
In general the
choice of slope algorithm does not make much difference, and you should just
stick with the default.
Two neighbors:
- Batson & others (1975):
calculates component of slope in direction of one of the
nearest neighbors. Strictly speaking this value of slope
should be assigned to a spot halfway between the point
and the neighbor used. This is the fastest algorithm,
with no calculation for aspect, and has limited
applicability. MICRODEM used it at one time for modeling
radar reflectivity from the DEM, as the look direction
for SLAR is typically to the east or west.
Three neighbors:
- O'Neill & Mark (1987):
three points uniquely define a plane; this algorithm uses
the point and its neighbors to the north and east (Z, Zn,
and Ze). The simplest measure, it has the lowest
correlations with all other methods. Other equally
arbitrary choices of points to use would give different
slopes and aspect directions. This method would be used
for TINs, where the slope and aspect apply to the
triangular facet and not the point at which the slope
will change. Called the simple method by
Jones (1998) who chose the points to the north and
west. This method assigns the values from the
triangle to the point in the center and assumes it covers
a rectangular grid. A more sophisticated, but
harder to implement, version would compute a different
slope for each of the eight triangular regions about the
point.
- Four contiguous right triangles (Onorati and others, 1992)
Sub-pixel calculations: using two neighbors about the
point, four subpixels can be computed. The four use
points (Z, Zn, Zw), (Z, Zn, Ze), (Z, Zs,
Zw), and (Z, Zs, Ze).
Four neighbors to N, S, E, and W (excluding point itself) (rook's case):
- Four neighbors (e.g. Ritter, 1987;
Zevenbergen and Thorne, 1987; Band, 1989; Eyton, 1991;
Carter, 1992) Called the Fleming & Hoffer (1979)
method by Jones (1998) for a Purdue University technical
report. Evans (1972) may have described this
method, and amplified it in Durham University technical
reports (Evans, 1979; Young, 1978); Shary (2002) calls
this the Evans-Young method. This gets a dx from Ze-Zw,
and a dy from Zn-Zs, and calculates the slope and aspect
from them. The value is assigned to the central
point, even though its elevation is not used in the
calculation.
Four neighbors to NW, SE,N E, and SW (excluding point itself) (bishop's case):
- Diagonal Ritters: proposed by Jones (1998), using the 4
diagonal neighbors. This gets a "dx" from
Znw-Zse, and a "dy" from Zne-Zsw, and
calculates the slope and aspect from them. This
effectively averages slopes over a longer distance, since
the diagonals are longer than the distances to the
nearest neighbors. The value is assigned to the central
point, even though its elevation is not used in the
calculation.
Eight neighbors (excluding point itself) (Queen's case):
- Horn method (Horn, 1981): nearest points weighted more
than diagonal neighbors. This method is also known as the
Sobel operator (Richards, 1986). The point itself has no
influence on the calculated slope (Guth,
1995). Jones (1998) discusses two variants, using 1/r
and 1/r² weightings. The value is assigned to the
central point, even though its elevation is not used in
the calculation.
- Sharpnack & Akin (Sharpnack and Akin, 1969; also
"mean method" in ELAS, cited in Ritter, 1987):
all neighbors weighted equally. The point itself has no
influence on the calculated slope (Guth,
1995) The value is assigned to the central
point, even though its elevation is not used in the
calculation.
- Local Trend Surface (Heerdegen and Beran, 1982): this
method is equivalent to fitting a second order trend
surface (z = a + bx + cy + dx² + ey² + fxy) to the nine
points using the algorithm in Davis (1973), but an order
of magnitude faster. The algorithm was modified for
MICRODEM to allow different x and y data spacing, but
gives the same results as the Sharpnack & Akin
algorithm. The point itself has no influence on the
calculated slope (Guth, 1995).The
value is assigned to the central point, even though its
elevation is not used in the calculation.
- Constrained quadratic surface method (Wood, 1996 thesis,
cited in Jones, 1998). The trend surface is constrained
to go through the point itself. The value is
assigned to the central point.
Nine points:
- Steepest Adjacent Neighbor (Collins,
1975; Travis et al. 1975, Dept of Agriculture Forest
Service Technical Report cited in Jones 1998; ELAS
software cited in Ritter, 1987): this method
consistently gives steeper results than all others. It
picks the steepest of the eight adjacent points,
considering that the diagonal distance is longer and that
the x and y spacings might be different. The aspect is
only available to the nearest 45° (approximate if the x
and y spacings differ), and if the steepest slope is
uphill is assigned as the opposite of the slope
direction, which may not be the steepest downhill
slope. The value is assigned to the central point.
- slope1 = (Zn-Z) / (NS Spacing)
- slope2 = (Zs-Z) / (NS Spacing)
- slope3 = (Ze-Z) / (EW Spacing)
- slope4 = (Zw-Z) / (EW Spacing)
- slope5 = (Znw-Z) / (Diagonal Spacing)
- slope6 = (Zne-Z) / (Diagonal Spacing)
- slope7 = (Zse-Z) / (Diagonal Spacing)
- slope8 = (Zsw-Z) / (Diagonal Spacing)
- Pick largest absolute value from eight candidates
as the slope
- Steepest downhill neighbor (D8): similar to the steepest
neighbor, but the point must be lower in elevation than
the central point. The value is assigned to the
central point. Most drainage basin algorithms use a
variant of this algorithm, looking only for the neighbor
into which water will flow.
- Average neighbor: averages slopes to the eight
nearest neighbors. Its estimate is about 50% of the
Steepest neighbor, because the slopes at 90° to the
maximum will typically be nearly horizontal. The
value is assigned to the central point.
Table. Selected Slope and Aspect Algorithms
expanded and adapted from Guth (1995).
Name |
Short Name |
Neighbors Used |
Neighbors |
Comments |
Slope
references
|
Steepest Down Hill |
SDN |
9 |
|
Also called D8; only 8 aspects |
O'Callahan and Mark, 1987 |
Steepest Adjacent Neighbor |
SAN |
9 |
|
Only 8 aspects |
Sharpnack & Akin, 1969 |
Average Neighbor |
AVN |
9 |
|
Only 8 aspects |
|
Guth Hybrid (Steepest + 8 even for aspect) |
HYB |
9 |
|
Preferred in Guth, 1995 |
Guth, 1995 |
Four Neighbors |
N4 |
4 |
zn,zs,ze,zw |
"Rook's case": Preferred in Jones, 1998 Second
order finite difference
|
Fleming and Hoffer, 1979; Ritter, 1987; Zevenbergen
and Thorne, 1987; O'Neill & Mark, 1987 |
3 Neighbors |
N3 |
3 |
z,zn,zne |
In a grid, two choices for the triangles |
O'Neill and Mark, 1987 |
8 neighbors, even weighting Evans or Evans-Young or Sharpnack and
Akin |
N8E |
8 |
|
"Queen's case": Equivalent to fitting a
second order trend surface. Preferred by Evans, 1998;
Florinsky, 1998. Identical weights. Third order finite difference
|
Sharpnack and Akin, 1969; Evans, 1980; Heerdegen and Beran, 1982; Wood, 1996;
Evans, 1998; Florinsky, 1998; Pike and others, 2009; Olaya, 2009 |
8 neighbors weighted Horn |
N82 |
8 |
|
Third order finite difference weighted by reciprocal
of twice distance |
Differential weights (Horn, 1981); Sobel operator (Richards, 1986)
|
8 neighbors weighted |
N8R |
8 |
|
Third order finite difference weighted by reciprocal of sqrt2 * distance |
Unwin, 1981 |
Frame Finite Difference |
FFD |
4 |
|
Difference above, below, right, and left |
|
Simple Difference |
SD |
3 |
|
|
The other choice of points from O'Neill and Mark, 1987 |
Two neighbors |
|
2 |
z,ze |
No aspect |
Batson & others (1975) |
Program |
Evans-Young Algorithm |
Horn Algorithm |
Zevenbergen and Thorne Algorithm |
Works with Geographic DEMs |
References |
Sharpnack and Akin, 1969; Evans, 1980; Heerdegen and Beran, 1982; Wood, 1996;
Evans, 1998; Florinsky, 1998; Pike and others, 2009; Olaya, 2009 |
Horn, 1981; Richards, 1986 |
Fleming and Hoffer, 1979; Ritter, 1987; Zevenbergen
and Thorne, 1987;O'Neill & Mark, 1987 |
Guth and Kane, 2021 |
MICRODEM |
Yes |
Yes |
Yes |
Yes |
GDAL |
|
Yes |
Yes |
No |
GRASS |
|
Yes |
|
? Strobl and others (2021) reported yes, but cannot verify now |
WhiteBox Tools |
|
Yes |
|
Yes |
ArcGIS |
|
Yes |
|
Geodetic option only |
Maps showing where slope algorithms differ.
Assess the effect of region size on
computed slope.
Slope
references
Slope units, in percent.
Slope
algorithms for arc second DEMs
Three slope methods appear to have the broadest support:
- Four closest neighbors (FCU):
- Eight neighbors unweighted (ENU): present preferred method in MICRODEM.
Users must still be aware that this method produces peaks that are too flat,
but its smoothing providew a "better" slope surface.
The Steepest Adjacent Neighbor (SAN): requires adjustments to
get reliable aspect distribution, and does not give a
smooth slope distribution. It was the preferred method in MICRODEM
until the summer of 2003, and there are still reasons to prefer it. The Eight Neighbors Unweighted method became the default in
MICRODEM the summer of 2003, primarily because it produces the smoothest,
most realistic slope histograms, and reasonable aspect distributions It still has problems with
smoothing of slopes in valleys, peaks, and ridges.
Guth (1995) reviewed all the
independent slope and aspect algorithms he could find. He
suggested that the steepest adjacent neighbor provided the best
estimate of slope, and the 8 neighbors with even weights provided
the best estimate of aspect. More recent comparisons of slope algorithms include Hodgson (1998) and Jones (1998), but neither cited Guth (1995).
Justification to use largest slope of 8 to adjacent elevations
in grid:
- Each point is surrounded by 8 neighbors. Four are to the
N,S,E,& W, and 4 approximately NE,SE,SW,& NW. How
close those are depend upon the type of DEM and latitude;
it is exact at the equator where x and y spacing is equal
or using a UTM based DEM. Withe DTED, the worst cases will be in the
vicinity of 50°, where just to the south the spacing is
60 and 90 m, and just to the north the spacing is 90 and
120 m. There the diagonal values are at 53 and 34°
instead of 45° (arc tan (120/90) and arc tan (60/90)).
- The direction of maximum slope can thus occur no more
than 28° away from the largest of the eight directional
slopes to nearest neighbors in the DTED grid (half of 56,
the complement of 34).
- If we assume a plane for the surface in the vicinity of
the greatest slope, we can use vector arithmetic to
project the slope in various directions. The projection
will be a function of the cosine. For the worst case, the
estimated slope to the DTED value in the closest
direction will be 87.9% of the true slope (cos(28) =
0.879), and that occurs only when the true maximum lies
directly between two grid values. At the equator, with
maximum error in estimating, the estimated slope is 92.4%
of its true value.
- Since maximum error is only about 10%, and given the
overall accuracy of the data to start with, we do not
feel it necessary to obtain the added rigor of fitting a
surface to the data points in the vicinity and then
computing the slope of the tangent surface.
- The other methods effectively compute the slope are a
region twice the size of the data spacing.
- See Guth (1995) for additional
details and a comparison on the effect of algorithm on
calculated slope.
Discussions of slope on arc-second DEMs
Last revised 2/9/2021