Adding Terrain Derivatives as Predictors

Derive on the DEM’s own grid, encode aspect circularly, then resample onto the imagery grid:

import numpy as np
import rioxarray

dem = rioxarray.open_rasterio("dem_30m.tif", masked=True).squeeze("band")
res = abs(dem.rio.resolution()[0])
gy, gx = np.gradient(dem.values, res)                       # metres per metre

slope = np.degrees(np.arctan(np.hypot(gx, gy)))
aspect = np.arctan2(-gx, gy)                                # radians, circular
features = {"slope": slope, "aspect_sin": np.sin(aspect), "aspect_cos": np.cos(aspect)}

Order matters more than the formulas: derive first, resample second. This page belongs to feature engineering for pixel-based models in Raster Machine Learning & Model Inference.

Why Terrain Carries Signal Satellites Cannot

Two pixels with identical reflectance can be entirely different surfaces if one sits on a steep north-facing slope and the other on a valley floor. Terrain tells a model about water availability, insolation, land use pressure and soil development — all of which drive what grows where, and none of which are legible in a single date of reflectance.

Terrain separates surfaces that look alike A hillside cross-section shows a shaded north-facing slope, a dry south-facing slope and a valley floor. Pixels on the shaded slope and in the valley can share a reflectance value while being different cover types; slope, aspect and elevation give the model the information needed to tell them apart. One reflectance value, three different surfaces shaded slope aspect 340 deg valley floor slope 2 deg dry slope aspect 160 deg All three green markers can share a reflectance vector; only terrain distinguishes them.

The practical consequence is that terrain features usually rank in the middle of an importance list — not the strongest predictors, but the ones that resolve the cases the spectral features cannot.

Environment & Setup

Package Version pin Used for
rioxarray >=0.15 DEM loading and reproject_match onto the imagery grid
numpy >=1.23 Gradients, curvature and the circular encoding
rasterio >=1.3.0 Writing the derivative stack
scipy >=1.11 Optional smoothing before differentiating a noisy DEM 
pip install "rioxarray>=0.15" "numpy>=1.23" "rasterio>=1.3.0" "scipy>=1.11"

Complete Working Example

import numpy as np
import rioxarray
import xarray as xr


def terrain_features(dem_path: str, reference: xr.DataArray,
                     *, smooth: bool = True) -> dict[str, xr.DataArray]:
    """Slope, aspect and curvature derived natively, then matched to the imagery grid."""
    dem = rioxarray.open_rasterio(dem_path, masked=True).squeeze("band")

    if dem.rio.crs != reference.rio.crs:
        # Derivatives need a projected CRS in metres, not degrees
        dem = dem.rio.reproject(reference.rio.crs)

    if smooth:
        from scipy import ndimage
        vals = ndimage.uniform_filter(np.nan_to_num(dem.values, nan=0.0), size=3)
        vals[~np.isfinite(dem.values)] = np.nan
    else:
        vals = dem.values

    res = abs(dem.rio.resolution()[0])
    gy, gx = np.gradient(vals, res)

    slope_deg = np.degrees(np.arctan(np.hypot(gx, gy))).astype("float32")
    aspect_rad = np.arctan2(-gx, gy)

    # Second derivatives give profile curvature: convex ridges vs concave hollows
    gyy, _ = np.gradient(gy, res)
    _, gxx = np.gradient(gx, res)
    curvature = (gxx + gyy).astype("float32")

    native = {
        "elevation": dem.astype("float32"),
        "slope": dem.copy(data=slope_deg),
        "aspect_sin": dem.copy(data=np.sin(aspect_rad).astype("float32")),
        "aspect_cos": dem.copy(data=np.cos(aspect_rad).astype("float32")),
        "curvature": dem.copy(data=curvature),
    }
    # Resample the DERIVATIVES, never the DEM first
    return {k: v.rio.reproject_match(reference) for k, v in native.items()}

Smoothing before differentiating is a judgement call that depends on the DEM. A global 30 m product carries enough vertical noise that raw gradients produce a speckled slope raster; a 3×3 mean removes most of it at the cost of flattening genuinely sharp features such as terraces and escarpments. For a lidar-derived surface, skip the smoothing entirely.

Derive First, Resample Second

Why upsampling the DEM first flattens the terrain A 30 metre elevation profile has a sharp break of slope. Bilinear upsampling to 10 metres rounds the break into a curve, and slope computed from the upsampled surface underestimates the true gradient by roughly a third at the steepest point. Computing slope on the native grid and then upsampling the slope preserves the peak. elevation profile slope from each route solid: native 30 m · dashed: bilinear to 10 m the break of slope is rounded away green: derive then resample — peak preserved red: resample then derive — peak cut by a third

The bias is systematic, not random: upsampling always smooths, so slope from an upsampled DEM is always too low, and always most wrong exactly where the terrain is steepest. A model trained on such features will systematically under-weight steep ground, and the error is invisible in the training metrics because the same bias applies at inference.

There is a second, subtler reason to derive natively. Slope depends on the pixel size used in the gradient, so the same DEM differentiated at 30 m and at 10 m gives different numbers for the same hillside — the 10 m version reports local micro-relief, the 30 m version reports the hillside. Deriving once, at a documented resolution, makes the feature mean one thing across the whole project. The general alignment mechanics are in aligning two rasters with reproject_match, and the resampling choices in choosing the right resampling method for Sentinel-2.

Verification

Terrain features fail quietly, so check them against physical expectations rather than against a schema.

Aspect on a circle, not on a line On a number line, north at 359 degrees and north at 1 degree sit at opposite ends despite describing the same direction. Encoded as a point on the unit circle through sine and cosine, the two values are adjacent, and a model splitting on either coordinate splits a meaningful half of the compass. The same two slopes, two encodings 1 deg 359 deg maximally distant as numbers raw degrees adjacent, as they should be sine and cosine pair 
import numpy as np

assert np.nanmax(slope) <= 90, "slope above vertical means the CRS is in degrees"
assert np.nanpercentile(slope, 99) < 60, "99th percentile slope implausibly steep"
assert np.isclose(np.nanmax(aspect_sin ** 2 + aspect_cos ** 2), 1.0, atol=1e-5)
print("mean slope:", float(np.nanmean(slope)), "deg")

The first assertion catches the single most common defect: differentiating a DEM that is still in geographic coordinates, where the pixel size is roughly 0.00027 degrees and np.gradient divides by that, producing slopes in the thousands. The unit circle check confirms the aspect encoding survived resampling — bilinear resampling of a sine and a cosine band separately does not exactly preserve the norm, and a large deviation indicates nearest-neighbour would have been the better choice.

Because the encoding is two bands rather than one, both must travel together through every later step; dropping aspect_sin while keeping aspect_cos leaves a feature that cannot distinguish east from west. Treat the pair as a single feature for pruning purposes.

A final sanity check that costs nothing: correlate slope against elevation. In most landscapes they are weakly positively correlated; a strong negative correlation usually means the DEM’s vertical sign convention is inverted, which happens with some bathymetry-merged products.

Common Errors

Slope values in the hundreds or thousands

The DEM is in a geographic CRS, so np.gradient divided metres of elevation by degrees of distance. Reproject to a metre-based CRS before differentiating, using the approach in reprojecting a raster from UTM to WGS84 in reverse.

Terrain features are blocky at the imagery resolution

Nearest-neighbour resampling was used on a continuous surface. Use bilinear for elevation, slope and curvature; reserve nearest for categorical terrain classes.

Aspect features make the model worse

Raw degrees were used instead of the sine and cosine pair, so the model is trying to split a circular variable on a linear axis. Replace the single aspect band with the two-band encoding.

Frequently Asked Questions

Q: Should I resample the DEM first or compute slope first? Compute first. Slope derived from an upsampled surface is smoothed twice and systematically underestimates steep ground. Derive on the native grid, then resample the derivative, which also makes the units unambiguous.

Q: Why encode aspect as sine and cosine? Because aspect is circular: 359 degrees and 1 degree are adjacent on the ground but maximally distant as numbers. A tree splitting on raw aspect wastes splits recovering that, while sine and cosine place neighbouring directions next to each other.

Q: Do terrain features help in flat landscapes? Rarely for slope and aspect, but elevation and a wetness index can still carry real signal, because a metre of relief in a floodplain separates wet from dry ground. Check the feature importance rather than assuming either way.