Computing Slope, Aspect and TRI from a Forest DTM in Python

Terrain derivatives are the cheapest ecological predictors you will ever produce: one bare-earth raster yields slope, aspect, ruggedness, curvature, and a wetness proxy, and between them they explain a large share of the variation in soil moisture, insolation, windthrow risk, and species occurrence. They are also the easiest layers to get quietly wrong, because every one of them is silently sensitive to cell size, to the CRS, and — for aspect — to the fact that it is a circular quantity that cannot be averaged. This guide derives them correctly from the surface produced by Digital Terrain Model Generation, inside Canopy Height Modeling & Terrain Extraction.

When to use each derivative

For habitat modelling, the usual minimum set is slope, northness, eastness, and TPI at one or two neighbourhood sizes — a combination that feeds directly into Environmental Predictor Stacking.

The gradient, and why the CRS decides everything

Slope and aspect both come from the surface gradient, estimated over a 3 × 3 window. With cell size $c$ in the same units as elevation:

$$\text{slope}=\arctan \sqrt{{\left(\frac{\partial z}{\partial x}\right)}^2+{\left(\frac{\partial z}{\partial y}\right)}^2},\text{aspect}=\arctan 2\left(\frac{\partial z}{\partial y},-\frac{\partial z}{\partial x}\right)$$

The requirement hidden in “the same units” is absolute. If the DTM is in a geographic CRS the horizontal step is degrees while the vertical is metres, so the ratio is off by a factor of roughly 111,000 and every slope comes back near zero. Reproject to a metre-based projected CRS first — the discipline covered in Coordinate Reference Systems for Forestry.

Minimal reproducible example

The Horn method — the same one behind gdaldem and most GIS software — is a dozen lines of numpy and avoids a dependency for the two derivatives you will always need.

import numpy as np
import rasterio


def slope_aspect(dtm_path: str):
    """Horn slope (degrees) and aspect (degrees clockwise from north)."""
    with rasterio.open(dtm_path) as src:
        z = src.read(1, masked=True).astype("float64")
        c = abs(src.transform.a)
        profile = src.profile
        assert not src.crs.is_geographic, "reproject to a metre-based CRS first"

    z = z.filled(np.nan)
    # Horn's 3x3 kernel, applied by shifting rather than convolving so NaNs stay local.
    def sh(dy, dx):
        return np.roll(np.roll(z, dy, axis=0), dx, axis=1)

    dzdx = ((sh(-1, -1) + 2 * sh(0, -1) + sh(1, -1))
            - (sh(-1, 1) + 2 * sh(0, 1) + sh(1, 1))) / (8 * c)
    dzdy = ((sh(-1, -1) + 2 * sh(-1, 0) + sh(-1, 1))
            - (sh(1, -1) + 2 * sh(1, 0) + sh(1, 1))) / (8 * c)

    slope_deg = np.degrees(np.arctan(np.hypot(dzdx, dzdy)))
    aspect_deg = np.degrees(np.arctan2(dzdy, -dzdx))
    aspect_deg = (450.0 - aspect_deg) % 360.0
    aspect_deg[np.hypot(dzdx, dzdy) < 1e-9] = -1.0     # flat: aspect undefined
    return slope_deg, aspect_deg, profile

Flat cells get -1 rather than an arbitrary bearing. Leaving them at 0 tells every downstream model that every flat pixel faces due north, which is a fabricated north-facing signal across valley bottoms and lake margins.

Aspect cannot enter a model as a number, so decompose it into two linear components before use:

def northness_eastness(aspect_deg: np.ndarray, slope_deg: np.ndarray):
    """Aspect as two continuous predictors, weighted by slope so flats read as zero."""
    rad = np.radians(aspect_deg)
    weight = np.sin(np.radians(slope_deg))       # flat ground carries no aspect signal
    north = np.cos(rad) * weight
    east = np.sin(rad) * weight
    flat = aspect_deg < 0
    north[flat] = 0.0
    east[flat] = 0.0
    return north, east

Weighting by the sine of slope is what makes these behave: a 1° slope facing north is ecologically almost identical to flat ground, and an unweighted northness would give both the value 1.

Ruggedness and position

TRI is the mean absolute elevation difference between a cell and its eight neighbours; TPI is the difference between a cell and the mean of a wider neighbourhood. Both are in metres and both depend entirely on the window you choose.

from scipy.ndimage import generic_filter, uniform_filter


def tri(z: np.ndarray) -> np.ndarray:
    """Terrain Ruggedness Index: mean absolute difference to the 8 neighbours."""
    def _tri(w):
        centre = w[4]
        return np.nanmean(np.abs(np.delete(w, 4) - centre))
    return generic_filter(z, _tri, size=3, mode="nearest")


def tpi(z: np.ndarray, radius_px: int = 10) -> np.ndarray:
    """Topographic Position Index: cell elevation minus its neighbourhood mean.

    Positive = ridge or convex slope; negative = valley or hollow.
    """
    size = 2 * radius_px + 1
    return z - uniform_filter(z, size=size, mode="nearest")

Report the radius with every TPI layer. A TPI at 50 m separates individual gullies; the same index at 500 m separates whole hillslopes from valley floors, and the two are different predictors that happen to share a name.

Parameter reference

Expected output and verification

Three assertions catch nearly every real error in a terrain-derivative stack:

import numpy as np

slope, aspect, profile = slope_aspect("dtm_1m.tif")

assert np.nanmax(slope) <= 90.0, "slope above 90 degrees is arithmetically impossible"
assert np.nanmedian(slope) > 0.5, "near-zero slope everywhere — the DTM is probably geographic"
finite = aspect[aspect >= 0]
assert finite.min() >= 0 and finite.max() < 360, "aspect outside the 0-360 range"

# Aspect must be roughly uniform over a large area of mixed terrain; a spike at one
# bearing means flats were assigned a direction instead of being flagged.
hist, _ = np.histogram(finite, bins=36, range=(0, 360))
assert hist.max() < 4 * np.median(hist), "aspect histogram spike — check the flat-cell handling"

The aspect histogram check is the one that catches the subtle bug. A stack that fails it will still produce plausible maps and a model that fits well, while carrying a fabricated directional signal wherever the ground is level.

For the wetness index, fill depressions first — an unfilled DTM has interior sinks where flow accumulation stops and TWI collapses:

import richdem as rd


def twi(dtm_path: str) -> np.ndarray:
    """Topographic Wetness Index from a depression-filled DTM."""
    dem = rd.LoadGDAL(dtm_path)
    rd.FillDepressions(dem, epsilon=True, in_place=True)
    acc = rd.FlowAccumulation(dem, method="Dinf")
    slope_rad = np.radians(np.asarray(rd.TerrainAttribute(dem, attrib="slope_degrees")))
    tan_beta = np.tan(np.maximum(slope_rad, np.radians(0.1)))   # guard against flats
    cell = abs(rd.LoadGDAL(dtm_path).geotransform[1])
    return np.log((np.asarray(acc) * cell + cell) / tan_beta)

The np.maximum guard is not cosmetic: tan(0) is zero, and dividing by it turns every flat cell into infinity, which then propagates through any normalisation applied to the whole stack.

Common pitfalls

• Deriving slope from a geographic CRS. The single most common error, and it produces a plausible-looking near-flat raster rather than an obvious failure.
• Averaging aspect. Resampling, zonal means, or a focal filter over raw aspect degrees produces nonsense at the 0/360 wrap. Decompose to northness and eastness first, average those, and recombine only if you must.
• Comparing TRI across resolutions. TRI at 1 m and TRI at 10 m are different quantities with the same name. Resample every input to a common grid before stacking.
• Deriving terrain from the CHM by mistake. The canopy surface produces slope values driven by crown edges rather than by ground. Confirm the input is bare earth.
• Forgetting to fill depressions before TWI. Real LiDAR DTMs are full of small sinks, some genuine and some artifacts, and every one truncates flow accumulation upstream of it.

Frequently Asked Questions

Related

• Digital Terrain Model Generation — produces the bare-earth surface these derivatives require
• Comparing TIN, IDW and Kriging for DTM Interpolation — the interpolation choice that sets how noisy your slope raster will be
• Environmental Predictor Stacking — where these layers become model covariates
• Coordinate Reference Systems for Forestry — why a projected CRS is a hard requirement here

Up: Digital Terrain Model Generation · Canopy Height Modeling & Terrain Extraction