Slope Algorithms

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

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:

1. 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)).
2. 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).
3. 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.
4. 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.
5. The other methods effectively compute the slope are a region twice the size of the data spacing.
6. See Guth (1995) for additional details and a comparison on the effect of algorithm on calculated slope.

Discussions of slope on arc-second DEMs

• https://gis.stackexchange.com/questions/14750/using-srtm-global-dem-for-slope-calculation/229186#229186),
• Equi-7
  • https://cartography.tuwien.ac.at/eurocarto/wp-content/uploads/2015/09/3_6_ppt.pdf
  • https://github.com/TUW-GEO/Equi7Grid/blob/master/docs/doc_files/EQUI7_V13_grid_and_tiling_definition_document_issue06.pdf
• https://github.com/r-barnes/dggridR/

Last revised 2/9/2021