The term gradient refers to a vector quantity, i.e. an object that has both magnitude and direction. The magnitude or size of the gradient is the slope, whilst the direction in which the maximum value of this magnitude occurs is known as the aspect. Slope is one of the most widely used of surface attributes so understanding how it is generated and what alternatives exist is extremely important.

Slope

The common notion of the slope of a surface or terrain, t, is the amount of rise (e.g. change in elevation) over run in some direction, usually the direction for which the rise over run is greatest, i.e. the steepest path up or down the surface. In practice we would specify some distance over which the slope is to be calculated, such as the amount of rise over a 10 meter distance in the plane. If this was 2 meters, for example, the rise over run would be S=2:10 or 1:5 or 0.2 or 20%. This ratio is actually the tangent function: S=tan(St) where St is the terrain slope in radians. If the run were taken as the surface distance rather than the plane distance then we would have S=sin(St). For gentle slopes (less than 1:4) this ambiguity in the definition of slope is not serious, but in general the tan() version is preferred.

The normal use of the term slope within GIS packages is the tan() value St, usually reported in degrees from the horizontal (11.3deg in this example). To convert between the two one would use an expression such as:

St=degrees(tan‑1(S))= 360tan‑1(S)/2π

For an analytical surface, F(x,y), slope is defined as the magnitude of the first derivative of the surface function:

This formula is not based on a rise-over-run calculation over a fixed interval, but rather assumes that a plane surface can be placed at any point on the surface F(x,y) in such a way that it only just touches the surface — it is a tangential plane and relies on the notion of infinitesimally small distances. It is closely related to the formula for Euclidean distance, and simply shows how much the surface F rises with a small fixed increment in x and y. However, surfaces within GIS are rarely if ever represented by analytic functions — typically they will be modeled as TINs or grids, with a finite resolution. Hence slope calculations will use approximations to the formula above depending on the surface model used, which is itself an approximate representation of the true surface. In almost all instances software packages apply grid-based computation, with TINs being first converted to grid format. Strictly speaking such conversion is not necessary, since TINs can provide interpolated surface values at every point, and use of a suitable (generally non-linear) function could then provide slope (and aspect) values directly.

The slope of a line crossing the surface (e.g. a road or pipeline) is dependent on the direction of the line relative to the underlying surface — it is always less than or equal to the slope of the surface. If the direction in question is an angle α we need to adjust the gradient vector by taking a dot product with a unit vector in the direction of interest. This enables slope values to be computed in any path direction. Such values are known as Directional derivatives. If the differential components for gradient and curvature have been calculated for a surface (e.g. a grid) then directional derivatives can be computed with very little further computation.

Using the grid model described in Section 6.1.3, Raster models, we have surface values represented by triples {x,y,z}, and can compare this with the mathematical formulation above to see that z(x,y) is equivalent to F(x,y). We have described this as an equivalence, but it should not be assumed that this is strictly true — the set of values z(x,y) is a particular representation of a non-mathematical surface or field at a pre-defined set of locations.

With a grid model slope can be determined at any particular location z* by simple approximation methods using the rook’s move or NSEW values for z at that point:

and St computed as before. This is the approach used by Surfer, as described in Section 6.1.3, Raster models. Within ArcGIS the 8-cell method due to Horn (1981), described in Section 6.1.3, Raster models, is used.

Methods based on fitting a plane surface to the unweighted 8-cell neighborhood of z* and using its slope as the surface slope estimate have been found to be less satisfactory. Even less acceptable are slopes computed by examining the maximum difference in cell values in the 8-cell neighborhood in the x- and y-directions and across the diagonals (suitably adjusted for cell size and diagonal length where appropriate).

Vertical Mapper for MapInfo uses yet another approach. It computes the slope and aspect of the 8 triangles defined by the mid-points of the 3x3 window centered on z* (Figure 6‑9). It then averages these to provide its estimates for slope and aspect at the point z*, in a similar manner to one of the methods for computing surface area (as illustrated in Figure 4‑5). The result is a measure that is less susceptible to directional bias, providing a full range of directions (infinitely many are possible). The 8 individual triangles of this local neighborhood each have distinct slope values and aspect (directions). These individual directions can be examined to see which is the steepest, and this slope value used in determining drainage flows across the surface. This approach, with the added element of assigning flows to two adjacent “downslope” cells, provides the basis for the method described by Tarboton (1997) as the D-infinity drainage method (see further, Section 6.4, Watersheds and Drainage). ERDAS uses yet another method, based on the average of the 3 east-west and 3 north-south differences for the estimates of dz/dx and dz/dy and then utilizes these as above for slope and aspect calculation.

Figure 6‑9 8-triangle slope computation

It will be apparent that these computations are extremely sensitive to scale. Larger grid intervals will result in smoothing of slope values, which may be desirable in some situations — in others it will be problematic. The availability of suitable datasets and the purposes to which the slope computations are to be put should be the basis of selection of method and resolution.

Some of these issues are illustrated in Figure 6‑10. Here the jagged black line represents the surface we are analyzing, which has some periodicity in its form. A grid resolution represented by the vertical green lines would always result in a slope estimate that was fairly constant and steep, whilst halving the interval (doubling the resolution) would give a rapidly alternating picture of the slopes. However, it may be that for the purposes of analysis a broader picture is required, as illustrated by the straight line, which would be generated if the underlying grid were smoothed (e.g. by a local averaging process) or if a much larger grid spacing were used.

Figure 6‑10 Gradient and sampling resolution

Figure 6‑11B shows the result of computing slope values in degrees for the source map displayed in Figure 6‑11A. The red and yellow areas are those of highest values (up to 65 degrees). The source raster has a resolution of 10 meters, and as the resolution is reduced/cell sizes for analysis are made larger, the maximum slope computed from this particular dataset falls to roughly 50 degrees with a 20m resolution and under 30 degrees with a 100m resolution. Slope computation is also affected to a greater or lesser degree by grid orientation, since it relies on grid values in particular directions.

Aspect

Aspect is defined as the directional component of the gradient vector and is the direction of maximum gradient of the surface at a given point. As with slope, aspect is calculated from estimates of the partial derivatives:

Aspect is computed in degrees from due north, i.e. as an azimuth in degrees not radians. The expression required is:

where the inverse tan function applies to both components — this function is often described, and implemented in software, using the notation atan2(a,b). Since the function is equivalent to atan(b/a) for a>0 the standard inverse tan function may be used if preferred for a>0.

Figure 6‑11 Slope computation output

A. Source terrain map |
B. Slope map |

As with the calculation of slope, Surfer uses simple finite differences for aspect calculation whereas ArcGIS uses the more complex 8-point formulas of Horn. This does not avoid problems of directional bias resulting from the computational procedure, grid structure and/or errors in or rounding of grid values. Figure 6‑12A shows the frequency distribution produced by the ArcGIS Aspect facility when applied to the terrain shown in Figure 6‑11A. Strong peaks in this distribution can be observed at approximately 45deg intervals, which of course does not reflect the true frequency distribution of landscape aspect. Landserf analysis of aspect for a subset of the region covered by Figure 6‑11A is shown in Figure 6‑12B (with zero values excluded). There are still discernible peaks, especially in the 8 major compass-point directions, but the computations appears to be smoother overall. Computations of aspect suffer from edge effects, resulting in zero or undefined values for boundary cells. Zero values are also generated wherever there are flat areas, and these can dominate frequency diagrams and associated analyses. Different packages treat flat areas in different ways, and these should be checked before relying too heavily on the results generated. Exclusion of zeros may be desirable in any further analysis of such data, and checks on the treatment of 360deg should also be made, although these do not seem to be a particular problem with most packages.

In Figure 6‑12C we have generated a circular frequency histogram of the set of aspect values again with zeros excluded. This is a more appropriate way of examining the data but is not available in most GIS packages. In this case the diagram was generated using the rose() function in MATLab with 60 intervals. Note that 0deg and 360deg in Figure 6‑12C are in the usual position for mathematical analysis, i.e. horizontal, and the log10 of frequencies rather than frequencies themselves have been plotted in order to reduce the effect of isolated peaks at certain frequencies.

Clearly there are some problems with this computational process, since areas of sea to the lower right of the output (Figure 6‑13A-C) are shown as having variable aspect — these are linear artifacts or ghost lines, which arise due to artifacts in the underlying source DEM.

Figure 6‑12 Frequency distribution of aspect values

A. ArcGIS Spatial Analyst: Aspect |
B. Landserf: Analyse|Surface parameter|Aspect |

C. Landserf Aspect data — Rose diagram (MATLab), Log freq.

Figure 6‑13 Aspect computation output

A. Source terrain map |
B. Aspect - classified |

C. Aspect — graduated colors |
D. Aspect — NSEW classification |

If a graduated color scale is used to represent this data, as in Figure 6‑13C, the pattern appears clearer, with high values shown in purple in this example, but the circular nature of the data means that this representation is misleading. A more appropriate visualization is one in which categories and colors reflect the circular nature of the data, as in Figure 6‑13D. Here the range 325deg‑ 45deg (broadly north-facing) has been colored pale orange, east-facing yellow, west-facing blue and southerly pale green.