DEM and Terrain Derivative Features

A digital elevation model is the single richest free covariate available to most geospatial machine learning projects. Elevation alone is rarely the useful signal, though — what a landslide model, a crop-yield model or a flood-risk model actually responds to is the shape of the land: how steeply it falls, which way it faces, whether it curves to collect water or shed it, and whether a pixel sits on a ridge or in a hollow. Those quantities are terrain derivatives, and they are computed from the elevation grid with small finite-difference kernels.

This topic is part of Spatial Feature Engineering for Machine Learning, which covers the full path from raw geodata to model-ready feature tables. Terrain derivatives sit alongside raster band math and index calculation as the two main sources of continuous raster predictors: band math describes what the surface is made of, terrain derivatives describe what shape it has. Both then feed the same downstream steps — zonal statistics for polygon aggregation if your samples are polygons, or direct point sampling if they are coordinates.

Problem Framing

Terrain derivatives look like a solved problem — every GIS ships a slope tool — and that is exactly why they go wrong in production pipelines. The tool hides three decisions that materially change the numbers it returns, and a model trained on one set of decisions will not survive a switch to another.

The first decision is units. Slope is rise / run. The rise comes from the elevation values; the run comes from the pixel size expressed in the raster’s own horizontal units. If the DEM is stored in EPSG:4326, that run is a fraction of a degree while the rise is in metres, so the ratio is not a slope at all — and because a degree of longitude shrinks towards the poles, the error is latitude-dependent. Reprojecting first, as covered in CRS alignment and projection handling, is not optional.

The second decision is kernel and scale. A 3×3 Horn kernel measures the gradient over one pixel width. On a 1 m lidar DEM that is a measurement of surface texture and noise; on a 90 m SRTM DEM it is a measurement of regional landform. The same variable name, slope_deg, means completely different things at the two resolutions, which is why resolution has to be pinned in configuration and recorded with the model artifact.

The third decision is representation. Slope, curvature and roughness are ordinary continuous variables. Aspect is not: it is an angle on a circle, and the discontinuity at north makes the raw degree value hostile to every common estimator. The fix — splitting aspect into northness and eastness — is a small piece of code that changes model performance far more than any hyperparameter, and it belongs in the same family of decisions as encoding categorical geographic features.

Prerequisites & Environment Setup

Terrain derivation needs a raster I/O stack plus SciPy for the moving-window filters. Pin the versions so kernel behaviour and nodata handling are reproducible.

# Pinned requirements for terrain derivatives
rasterio==1.3.10
numpy==1.26.4
scipy==1.13.1
xarray==2024.6.0
rioxarray==0.15.7
geopandas==0.14.4

Install with:

pip install "rasterio==1.3.10" "numpy==1.26.4" "scipy==1.13.1" \
            "xarray==2024.6.0" "rioxarray==0.15.7" "geopandas==0.14.4"

GDAL/PROJ system dependencies (Ubuntu/Debian):

sudo apt-get install -y gdal-bin libgdal-dev libproj-dev proj-data

The proj-data package matters more here than in most raster work: DEMs are frequently distributed on a geoid-based vertical datum (EGM96, EGM2008) while your horizontal CRS is ellipsoidal. If you plan to combine elevation from two sources, install the datum grids so pyproj can perform the vertical shift instead of silently offsetting one source by tens of metres.

Step-by-Step Implementation

Step 1 — Load the DEM and validate its geometry

Everything downstream depends on knowing the pixel size in metres. Read it from the affine transform rather than assuming it, and refuse to continue on an angular CRS.

import numpy as np
import rasterio
from rasterio.crs import CRS


def load_dem(path: str) -> tuple[np.ndarray, dict]:
    """Read a DEM as float32 with nodata converted to NaN.

    Returns:
        (elevation array, profile dict carrying transform, crs and pixel size).
    """
    with rasterio.open(path) as src:
        if src.crs is None:
            raise ValueError("DEM has no CRS; terrain derivatives would be meaningless")
        unit = src.crs.linear_units if src.crs.is_projected else "degree"
        if unit not in ("metre", "meter", "m"):
            raise ValueError(
                f"DEM horizontal unit is '{unit}'. Reproject to a metric CRS "
                "(for example a UTM zone) before computing gradients."
            )
        dem = src.read(1).astype("float32")
        if src.nodata is not None:
            dem[dem == src.nodata] = np.nan
        profile = src.profile.copy()
        xres, yres = abs(src.transform.a), abs(src.transform.e)

    if not np.isclose(xres, yres, rtol=1e-3):
        raise ValueError(f"Non-square pixels ({xres} x {yres}); resample before deriving slope")

    profile.update(cell_size=float(xres))
    return dem, profile


dem, profile = load_dem("dem_utm33n.tif")
assert np.isfinite(dem).any(), "DEM is entirely nodata"
print(f"DEM {dem.shape} at {profile['cell_size']} m, "
      f"{np.isnan(dem).mean():.2%} nodata")

Step 2 — Fill voids before differentiating

A finite-difference kernel reads nine cells. If one of them is NaN, the output is NaN, so a handful of voids in the source DEM smears into a much larger hole in every derivative band. Fill small voids first; leave large ones masked and let the model see them as missing.

from scipy import ndimage


def fill_small_voids(dem: np.ndarray, max_void_px: int = 500) -> np.ndarray:
    """Interpolate voids smaller than max_void_px, leave larger ones as NaN."""
    filled = dem.copy()
    void = np.isnan(dem)
    if not void.any():
        return filled

    labels, n = ndimage.label(void)
    sizes = ndimage.sum(void, labels, index=np.arange(1, n + 1))
    small = np.isin(labels, np.flatnonzero(sizes <= max_void_px) + 1)

    # Nearest-neighbour fill is deterministic and preserves the elevation range.
    idx = ndimage.distance_transform_edt(void, return_distances=False, return_indices=True)
    filled[small] = dem[tuple(i[small] for i in idx)]
    return filled


dem_filled = fill_small_voids(dem)
assert np.isnan(dem_filled).sum() <= np.isnan(dem).sum(), "void fill increased nodata"

Step 3 — Slope and aspect from the Horn kernel

Horn’s method weights the eight neighbours of each cell and is the formulation used by both GDAL and ArcGIS, so results stay comparable with the rest of the geospatial toolchain.

def slope_aspect(dem: np.ndarray, cell_size: float) -> tuple[np.ndarray, np.ndarray]:
    """Return (slope in degrees, aspect in degrees clockwise from north)."""
    # Horn 3x3 weights
    kx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype="float32")
    ky = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype="float32")

    dz_dx = ndimage.convolve(dem, kx, mode="nearest") / (8.0 * cell_size)
    dz_dy = ndimage.convolve(dem, ky, mode="nearest") / (8.0 * cell_size)

    slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
    aspect = np.degrees(np.arctan2(dz_dy, -dz_dx))
    aspect = (450.0 - aspect) % 360.0          # convert to compass bearing
    aspect[slope < 1e-6] = np.nan              # aspect is undefined on flat ground
    return slope.astype("float32"), aspect.astype("float32")


slope, aspect = slope_aspect(dem_filled, profile["cell_size"])
assert np.nanmax(slope) <= 90.0, "slope exceeded 90 degrees — check cell size units"
assert np.nanmin(aspect) >= 0.0 and np.nanmax(aspect) < 360.0

Step 4 — Make aspect model-safe

Feed the model the two circular components, never the raw bearing. Flat cells get zero for both, which correctly says no preferred direction.

def aspect_components(aspect_deg: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Split a compass bearing into northness and eastness in [-1, 1]."""
    rad = np.radians(aspect_deg)
    northness = np.cos(rad)
    eastness = np.sin(rad)
    flat = np.isnan(aspect_deg)
    northness[flat] = 0.0
    eastness[flat] = 0.0
    return northness.astype("float32"), eastness.astype("float32")


northness, eastness = aspect_components(aspect)
assert np.nanmax(np.abs(northness)) <= 1.0 and np.nanmax(np.abs(eastness)) <= 1.0

Step 5 — Curvature, topographic position and roughness

Curvature separates water-shedding convex ground from water-collecting concave ground. Topographic position index compares each cell to the mean of its neighbourhood and is the cheapest way to encode ridge-versus-valley. Roughness is the local elevation range and captures surface texture.

def curvature(dem: np.ndarray, cell_size: float) -> np.ndarray:
    """Laplacian curvature: positive = convex (shedding), negative = concave (collecting)."""
    lap = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]], dtype="float32")
    return (ndimage.convolve(dem, lap, mode="nearest") / (cell_size ** 2)).astype("float32")


def tpi(dem: np.ndarray, window: int = 11) -> np.ndarray:
    """Topographic position index: cell elevation minus the neighbourhood mean."""
    if window % 2 == 0:
        raise ValueError("TPI window must be odd so the cell sits at the centre")
    mean = ndimage.uniform_filter(dem, size=window, mode="nearest")
    return (dem - mean).astype("float32")


def roughness(dem: np.ndarray, window: int = 3) -> np.ndarray:
    """Local elevation range within the window."""
    hi = ndimage.maximum_filter(dem, size=window, mode="nearest")
    lo = ndimage.minimum_filter(dem, size=window, mode="nearest")
    return (hi - lo).astype("float32")


curv = curvature(dem_filled, profile["cell_size"])
tpi_11 = tpi(dem_filled, window=11)
tpi_33 = tpi(dem_filled, window=33)
rough = roughness(dem_filled, window=3)
assert np.isfinite(np.nanmean(tpi_11)), "TPI is all NaN — check the void fill step"

Step 6 — Stack and write a multi-band feature raster

Write one file, not eight, so the feature set travels as a single versioned artifact and the band order can never drift between training and inference.

BAND_NAMES = ["elevation", "slope_deg", "northness", "eastness",
              "curvature", "tpi_11", "tpi_33", "roughness"]


def write_terrain_stack(path: str, bands: list[np.ndarray], profile: dict) -> None:
    out = profile.copy()
    out.pop("cell_size", None)
    out.update(count=len(bands), dtype="float32", nodata=np.nan,
               compress="deflate", tiled=True, blockxsize=512, blockysize=512)
    with rasterio.open(path, "w", **out) as dst:
        for i, (arr, name) in enumerate(zip(bands, BAND_NAMES), start=1):
            dst.write(arr.astype("float32"), i)
            dst.set_band_description(i, name)


write_terrain_stack(
    "terrain_features.tif",
    [dem_filled, slope, northness, eastness, curv, tpi_11, tpi_33, rough],
    profile,
)

Verification & Testing

Terrain code fails quietly, so test it against surfaces whose answer you already know. A synthetic tilted plane has a slope you can compute by hand, and that single test catches every units mistake.

import numpy as np


def test_slope_on_synthetic_plane():
    """A plane rising 10 m over 100 m must yield slope = atan(0.1) = 5.71 degrees."""
    cell = 10.0
    rows, cols = 50, 50
    yy, xx = np.mgrid[0:rows, 0:cols]
    plane = (xx * cell) * 0.1          # 10% grade towards +x

    slope, aspect = slope_aspect(plane.astype("float32"), cell)
    interior = slope[2:-2, 2:-2]

    assert np.allclose(interior, np.degrees(np.arctan(0.1)), atol=1e-3)
    # A surface rising towards +x faces west, i.e. a bearing of 270 degrees.
    assert np.allclose(aspect[2:-2, 2:-2], 270.0, atol=1e-3)


def test_tpi_sums_to_zero_on_a_plane():
    """On any planar surface every cell equals its neighbourhood mean."""
    plane = np.fromfunction(lambda r, c: 3.0 * r + 2.0 * c, (60, 60), dtype="float32")
    assert np.allclose(tpi(plane, window=9)[6:-6, 6:-6], 0.0, atol=1e-3)

For a visual sanity check, render a hillshade from your slope and aspect arrays and compare it against a basemap. Shadows that fall on the wrong side of every ridge mean the dz_dy sign convention is inverted — a common error when a DEM is read with a north-up versus south-up transform.

def hillshade(slope_deg, aspect_deg, azimuth=315.0, altitude=45.0):
    """Standard hillshade for visual QA of slope and aspect correctness."""
    az, alt = np.radians(360.0 - azimuth + 90.0), np.radians(altitude)
    s, a = np.radians(slope_deg), np.radians(aspect_deg)
    shaded = np.sin(alt) * np.cos(s) + np.cos(alt) * np.sin(s) * np.cos(az - a)
    return np.clip(shaded * 255.0, 0, 255).astype("uint8")

Troubleshooting & Common Errors

ValueError: DEM horizontal unit is 'degree' — raised by the loader above. The DEM is geographic. Reproject it with rasterio.warp.reproject to a metric CRS covering your study area before any gradient work; see reprojecting a raster with rasterio warp reproject.

Slope maxes out at exactly 90 degrees along tile edges — the DEM was mosaicked from tiles with mismatched nodata values, so a cliff of -9999 sits next to real elevation. Convert nodata to NaN at read time (Step 1) rather than relying on the sentinel value surviving arithmetic.

Every derivative band is entirely NaN — one NaN anywhere propagates through ndimage.convolve. Confirm the void fill ran, and confirm the fill did not itself introduce NaN by nearest-neighbour sampling from another void.

Aspect is rotated by 90 or 180 degrees compared with a reference hillshade — the dz_dy sign convention depends on whether row 0 is the north or the south edge. Rasterio uses a north-up transform with a negative e term, which the code above assumes. If your array was flipped with np.flipud, negate dz_dy.

Curvature values are astronomically large — curvature scales with 1 / cell_size², so a 1 m DEM produces numbers a thousand times larger than a 30 m DEM. This is correct behaviour, not a bug, but it means curvature must be rescaled with the rest of the feature matrix; see feature scaling for geospatial inputs.

Model accuracy collapses when a new DEM tile is added — a different source (SRTM versus a national lidar product) has a different vertical datum and noise floor, so slope and roughness distributions shift. Treat a DEM source change as a covariate shift event and check it with model drift detection for geospatial inference.

Performance Optimisation

A national DEM at 10 m resolution is tens of gigabytes, well past what fits in memory as float32 with eight derivative bands. Three techniques carry the workload.

Process in overlapping blocks. Every kernel here has a finite footprint: 3×3 for slope, window for topographic position index. Read a block plus a halo of max_window // 2 pixels, compute, then discard the halo before writing. Rasterio’s block_windows() gives you the native tiling, and matching it avoids re-reading compressed strips.

HALO = 33 // 2 + 1

with rasterio.open("dem_utm33n.tif") as src:
    for _, window in src.block_windows(1):
        padded = window.round_offsets().round_lengths()
        padded = rasterio.windows.Window(
            max(0, padded.col_off - HALO), max(0, padded.row_off - HALO),
            padded.width + 2 * HALO, padded.height + 2 * HALO,
        ).intersection(rasterio.windows.Window(0, 0, src.width, src.height))
        tile = src.read(1, window=padded).astype("float32")
        # ... derive, then slice the halo away before writing

Keep everything float32. SciPy promotes to float64 whenever an input is float64, doubling both memory and time for no accuracy benefit at DEM precision. Cast once at read and assert the dtype after each derivative.

Reach for dask only when blocking is not enough. rioxarray.open_rasterio(..., chunks={"x": 2048, "y": 2048}) gives lazy arrays that xarray’s rolling operations can traverse in parallel, and the same graph runs on a multi-core laptop or a distributed cluster without code changes — the pattern described in scaling batch inference with Dask and Ray. Set chunks to a multiple of the GeoTIFF’s internal block size or every read will straddle two compressed tiles.

Cache the stack, never the code path. Terrain does not change between training runs. Compute the derivative stack once as a cloud-optimized GeoTIFF, register it with spatial dataset versioning with DVC and lakeFS, and have both training and inference read that artifact rather than recomputing.

FAQ

Why does slope look wrong when the DEM is in EPSG:4326?

Slope is a ratio of vertical change to horizontal distance. In EPSG:4326 the horizontal units are degrees while elevation is in metres, so the finite-difference denominator is meaningless — the values come out orders of magnitude off, and the error changes with latitude because a degree of longitude shrinks towards the poles. Reproject the DEM to a metric CRS such as the appropriate UTM zone before computing any gradient.

Should aspect be fed to a model as degrees?

No. Aspect is circular: 359° and 1° describe almost the same direction but sit at opposite ends of the numeric range, so a decision tree will happily split between them and a linear model will read a 358-unit jump. Decompose aspect into northness (cos) and eastness (sin) and feed those two bounded, continuous columns instead.

How large should the moving window be for topographic position index?

The window sets the landform scale the feature responds to. At 10 m resolution a 3×3 window captures gullies and micro-relief while a 33×33 window captures valley-bottom versus hilltop position. Since the informative scale is rarely known in advance, compute two or three window sizes and let feature importance — see model explainability for spatial predictions — tell you which scale the model actually uses.

Do terrain features need to be recomputed at inference time?

They should not be. Terrain is static on the timescale of a model’s life, so compute the derivative stack once, version it, and read the same artifact in training and serving. Recomputing at inference from a different DEM source, resolution or void-fill setting is a textbook training-serving skew: the model sees slope and roughness distributions it was never fitted on.

Related

• Raster Band Math and Index Calculation — the spectral counterpart to terrain derivatives
• Computing Slope and Aspect from a DEM with NumPy — the kernel implementation in detail
• Deriving Topographic Wetness Index for Flood Models — combining slope with upslope contributing area
• Zonal Statistics for Polygon Aggregation — reducing the terrain stack to per-polygon predictors
• Feature Scaling for Geospatial Inputs — normalising bands whose ranges differ by orders of magnitude

Part of: Spatial Feature Engineering for Machine Learning