Terrain attributes

xDEM can derive a wide range of terrain attributes from a DEM.

Attributes are derived in pure Python for modularity (e.g., varying window size) and other uses (e.g., uncertainty), and tested for consistency against gdaldem and RichDEM.

Quick use

Terrain attribute methods can be derived directly from a DEM, using for instance xdem.DEM.slope().

Show the opening of example files. Hide the opening of example files.

import xdem

# Open a DEM from a filename on disk
filename_dem = xdem.examples.get_path("longyearbyen_ref_dem")
dem = xdem.DEM(filename_dem)

# Slope from DEM method
slope = dem.slope()
# Or from terrain module using an array input
slope = xdem.terrain.slope(dem.data, resolution=dem.res)

Coming soon

We are working on further optimizing the computational performance of certain terrain attributes using convolution.

Summary of supported methods

Attribute	Unit (if DEM in meters)	Reference
Slope	Degrees (default) or radians	Horn (1981) or Zevenbergen and Thorne (1987)
Aspect	Degrees (default) or radians	Horn (1981) or Zevenbergen and Thorne (1987)
Hillshade	Unitless	Horn (1981) or Zevenbergen and Thorne (1987)
Curvature	Meters-1 * 100	Zevenbergen and Thorne (1987)
Planform curvature	Meters-1 * 100	Zevenbergen and Thorne (1987)
Profile curvature	Meters-1 * 100	Zevenbergen and Thorne (1987)
Topographic position index	Meters	Weiss (2001)
Terrain ruggedness index	Meters	Riley et al. (1999) or Wilson et al. (2007)
Roughness	Meters	Dartnell (2000)
Rugosity	Unitless	Jenness (2004)
Fractal roughness	Fractal dimension number (1 to 3)	Taud and Parrot (2005)

Note

Only grids with equal pixel size in X and Y are currently supported. Transform into such a grid with xdem.DEM.reproject().

Slope

xdem.DEM.slope()

The slope of a DEM describes the tilt, or gradient, of each pixel in relation to its neighbours. It is most often described in degrees, where a flat surface is 0° and a vertical cliff is 90°. No tilt direction is stored in the slope map; a 45° tilt westward is identical to a 45° tilt eastward.

The slope \(\alpha\) can be computed either by the method of Horn (1981) (default) based on a refined gradient formulation on a 3x3 pixel window, or by the method of Zevenbergen and Thorne (1987) based on a plane fit on a 3x3 pixel window.

For both methods, \(\alpha = arctan(\sqrt{p^{2} + q^{2}})\) where \(p\) and \(q\) are the gradient components west-to-east and south-to-north, respectively, with for Horn (1981):

\[\begin{split} p_{\textrm{Horn}}=\frac{(h_{++} + 2h_{+0} + h_{+-}) - (h_{-+} + 2h_{-0} + h_{--})}{8 \Delta x},\\ q_{\textrm{Horn}}=\frac{(h_{++} + 2h_{0+} + h_{-+}) - (h_{+-} + 2h_{0-} + h_{--})}{8 \Delta y}, \end{split}\]

and for Zevenbergen and Thorne (1987):

\[\begin{split} p_{\textrm{ZevTho}} = \frac{h_{+0} - h_{-0}}{2 \Delta x},\\ q_{\textrm{ZevTho}} = \frac{h_{0+} - h_{0-}}{2 \Delta y}, \end{split}\]

where \(h_{ij}\) is the elevation at pixel \(ij\), where indexes \(i\) and \(j\) correspond to east-west and north-south directions respectively, and take values of either the center (\(0\)), west or south (\(-\)), or east or north (\(+\)):

	West	Center	East
North	\(h_{-+}\)	\(h_{0+}\)	\(h_{++}\)
Center	\(h_{-0}\)	\(h_{00}\)	\(h_{+0}\)
South	\(h_{--}\)	\(h_{0-}\)	\(h_{+-}\)

Finally, \(\Delta x\) and \(\Delta y\) correspond to the pixel resolution west-east and south-north, respectively.

The differences between methods are illustrated in the Slope and aspect methods example.

slope = dem.slope()
slope.plot(cmap="Reds", cbar_title="Slope (°)")

Aspect

xdem.DEM.aspect()

The aspect describes the orientation of strongest slope. It is often reported in degrees, where a slope tilting straight north corresponds to an aspect of 0°, and an eastern aspect is 90°, south is 180° and west is 270°. By default, a flat slope is given an arbitrary aspect of 180°.

The aspect \(\theta\) is based on the same variables as the slope, and thus varies similarly between the method of Horn (1981) (default) and that of Zevenbergen and Thorne (1987):

\[ \theta = arctan(\frac{p}{q}), \]

with \(p\) and \(q\) defined in the slope section.

Warning

A north aspect represents the upper direction of the Y axis in the coordinate reference system of the input, xdem.DEM.crs, which might not represent the true north.

aspect = dem.aspect()
aspect.plot(cmap="twilight", cbar_title="Aspect (°)")

Hillshade

xdem.DEM.hillshade()

The hillshade is a slope map, shaded by the aspect of the slope. With a westerly azimuth (a simulated sun coming from the west), all eastern slopes are slightly darker. This mode of shading the slopes often generates a map that is much more easily interpreted than the slope.

The hillshade \(hs\) is directly based on the slope \(\alpha\) and aspect \(\theta\), and thus also varies between the method of Horn (1981) (default) and that of Zevenbergen and Thorne (1987), and is often scaled between 1 and 255:

\[ hs = 1 + 254 \left[ sin(alt) cos(\alpha) + cos(alt) sin(\alpha) sin(2\pi - azim - \theta) \right], \]

where \(alt\) is the shading altitude (90° = from above) and \(azim\) is the shading azimuth (0° = north).

Note, however, that the hillshade is not a shadow map; no occlusion is taken into account so it does not represent “true” shading. It therefore has little analytic purpose, but it still constitutes a great visualization tool.

hillshade = dem.hillshade()
hillshade.plot(cmap="Greys_r", cbar_title="Hillshade")

Curvature

xdem.DEM.curvature()

The curvature is the second derivative of elevation, which highlights the convexity or concavity of the terrain. If a surface is convex (like a mountain peak), it will have positive curvature. If a surface is concave (like a through or a valley bottom), it will have negative curvature. The curvature values in units of m^-1 are quite small, so they are by convention multiplied by 100.

The curvature \(c\) is based on the method of Zevenbergen and Thorne (1987).

\[c = - 100 \frac{(h_{+0} + h_{-0} + h_{0+} + h_{0-}) - 4 h_{00}}{\Delta x \Delta y}.\]

curvature = dem.curvature()
curvature.plot(cmap="RdGy_r", cbar_title="Curvature (100 / m)", vmin=-1, vmax=1, interpolation="antialiased")

Planform curvature

xdem.DEM.planform_curvature()

The planform curvature is the curvature perpendicular to the direction of slope, reported in 100 m^-1, also based on Zevenbergen and Thorne (1987):

\[planc = -2\frac{DH² + EG² -FGH}{G²+H²}\]

with:

\[\begin{split}D &= \frac{h_{0+} + h_{0-} - 2h_{00}} {2\Delta y^{2}}, \\ E &= \frac{h_{+0} + h_{-0} - 2h_{00}} {2\Delta x^{2}}, \\ F &= \frac{h_{--} + h_{++} - h_{-+} - h_{+-}} {4 \Delta x \Delta y}, \\ G &= \frac{h_{0-} - h_{0+}}{2\Delta y}, \\ H &= \frac{h_{-0} - h_{+0}}{2\Delta x}.\end{split}\]

planform_curvature = dem.planform_curvature()
planform_curvature.plot(cmap="RdGy_r", cbar_title="Planform curvature (100 / m)", vmin=-1, vmax=1, interpolation="antialiased")

Profile curvature

xdem.DEM.profile_curvature()

The profile curvature is the curvature parallel to the direction of slope, reported in 100 m^-1, also based on Zevenbergen and Thorne (1987):

\[profc = 2\frac{DG² + EH² + FGH}{G²+H²}\]

based on the equations in the planform curvature section for \(D\), \(E\), \(F\), \(G\) and \(H\).

profile_curvature = dem.profile_curvature()
profile_curvature.plot(cmap="RdGy_r", cbar_title="Profile curvature (100 / m)", vmin=-1, vmax=1, interpolation="antialiased")

Topographic position index

xdem.DEM.topographic_position_index()

The topographic position index (TPI) is a metric of slope position, described in Weiss (2001), that corresponds to the difference of the elevation of a central pixel with the average of that of neighbouring pixels. Its unit is that of the DEM (typically meters) and it can be computed for any window size (default 3x3 pixels).

\[ tpi = h_{00} - \textrm{mean}_{i\neq 0, j\neq 0} (h_{ij}) . \]

tpi = dem.topographic_position_index()
tpi.plot(cmap="Spectral", cbar_title="Topographic position index (m)", vmin=-5, vmax=5)

Terrain ruggedness index

xdem.DEM.terrain_ruggedness_index()

The terrain ruggedness index (TRI) is a metric of terrain ruggedness, based on cumulated differences in elevation between a central pixel and its surroundings. Its unit is that of the DEM (typically meters) and it can be computed for any window size (default 3x3 pixels).

For topography (default), the method of Riley et al. (1999) is generally used, where the TRI is computed by the squareroot of squared differences with neighbouring pixels:

\[ tri_{\textrm{Riley}} = \sqrt{\sum_{ij}(h_{00} - h_{ij})^{2}}. \]

For bathymetry, the method of Wilson et al. (2007) is generally used, where the TRI is defined by the mean absolute difference with neighbouring pixels:

\[ tri_{\textrm{Wilson}} = \textrm{mean}_{ij} (|h_{00} - h{ij}|) . \]

tri = dem.terrain_ruggedness_index()
tri.plot(cmap="Purples", cbar_title="Terrain ruggedness index (m)")

Roughness

xdem.DEM.roughness()

The roughness is a metric of terrain ruggedness, based on the maximum difference in elevation in the surroundings, described in Dartnell (2000). Its unit is that of the DEM (typically meters) and it can be computed for any window size (default 3x3 pixels).

\[ r_{\textrm{D}} = \textrm{max}_{ij} (h{ij}) - \textrm{min}_{ij} (h{ij}) . \]

roughness = dem.roughness()
roughness.plot(cmap="Oranges", cbar_title="Roughness (m)")

Rugosity

xdem.DEM.rugosity()

The rugosity is a metric of terrain ruggedness, based on the ratio between planimetric and real surface area, described in Jenness (2004). It is unitless, and is only supported for a 3x3 window size.

\[ r_{\textrm{J}} = \frac{\sum_{k \in [1,8]} A(T_{00, k})}{\Delta x \Delta y}. \]

where \(A(T_{00,k})\) is the area of one of the 8 triangles connected from the center of the center pixel \(00\) to the centers of its 8 neighbouring pixels \(k\), cropped to intersect only the center pixel. This surface area is computed in three dimensions, accounting for elevation differences.

rugosity = dem.rugosity()
rugosity.plot(cmap="YlOrRd", cbar_title="Rugosity")

Fractal roughness

xdem.DEM.fractal_roughness()

The fractal roughness is a metric of terrain ruggedness based on the local fractal dimension estimated by the volume box-counting method of Taud and Parrot (2005). The fractal roughness is computed by estimating the fractal dimension in 3D space, for a local window centered on the DEM pixels. Its unit is that of a dimension, and is always between 1 (dimension of a line in 3D space) and 3 (dimension of a cube in 3D space). It can only be computed on window sizes larger than 5x5 pixels, and defaults to 13x13.

fractal_roughness = dem.fractal_roughness()
fractal_roughness.plot(cmap="Reds", cbar_title="Fractal roughness (dimensions)")

Generating multiple attributes at once

Often, one may seek more terrain attributes than one, e.g. both the slope and the aspect. Since both are dependent on the gradient of the DEM, calculating them separately is computationally repetitive. Multiple terrain attributes can be calculated from the same gradient using the xdem.DEM.get_terrain_attribute() function.

Elevation error map

Elevation error map

Terrain attributes

Terrain attributes

Spatial propagation of elevation errors

Spatial propagation of elevation errors

Estimation and modelling of heteroscedasticity

Estimation and modelling of heteroscedasticity

Standardization for stable terrain as error proxy

Standardization for stable terrain as error proxy