Terrain Analysis and DEM-Derived Products

Elevation shapes almost everything a satellite sees. It decides how brightly a slope is lit, where water collects, what grows where, and how much of a scene lies in shadow. A digital elevation model and the layers derived from it — slope, aspect, hillshade, curvature — are therefore among the most frequently reused rasters in remote sensing, as features for models, as context for interpretation, and as inputs to correcting imagery itself. This topic belongs to Satellite Processing Workflows & Index Pipelines.

The technical ground rules are few but strict: work in a metre-based projected CRS, fill voids before differentiating, derive at the DEM’s native resolution, and only then resample onto the imagery grid.

---

Prerequisites

pip install "rasterio>=1.3.0" "rioxarray>=0.15" "numpy>=1.23" "scipy>=1.11" "pystac-client>=0.7"

You will need the reprojection and alignment tools from mastering CRS transformations in rasterio and aligning two rasters with reproject_match.

---

The Products and How They Relate

The dependency order matters because errors propagate. A void left unfilled becomes an undefined gradient that erases a patch of slope and aspect; a DEM left in degrees makes every derivative wrong by orders of magnitude. Getting the second box right is most of the work.

---

Step-by-Step Workflow

Step 1 — Assemble a DEM for the area

import rasterio
from rasterio.merge import merge


def mosaic_dem(tile_paths: list[str], bounds=None):
    sources = [rasterio.open(p) for p in tile_paths]
    try:
        arr, transform = merge(sources, bounds=bounds, nodata=-32768)
        profile = sources[0].profile | {"height": arr.shape[1], "width": arr.shape[2],
                                        "transform": transform, "nodata": -32768}
    finally:
        for s in sources:
            s.close()
    return arr[0], profile

Global DEMs ship as one-degree tiles in geographic coordinates, so an area of interest of any size usually needs several. Merging before reprojecting avoids seams at tile edges, and passing bounds limits the mosaic to the area actually needed. The merge mechanics are the same as in merging tiles with rasterio merge.

Step 2 — Project to metres and fill voids

import numpy as np
import rioxarray
from scipy import ndimage


def project_and_fill(dem_path: str, epsg: int, res: float = 30.0):
    dem = rioxarray.open_rasterio(dem_path, masked=True).squeeze("band")
    dem = dem.rio.reproject(f"EPSG:{epsg}", resolution=res, resampling=1)   # bilinear
    vals = dem.values.astype("float32")
    missing = ~np.isfinite(vals)
    if missing.any():
        # nearest-valid fill: the value of the closest non-missing cell
        idx = ndimage.distance_transform_edt(missing, return_distances=False, return_indices=True)
        vals = vals[tuple(idx)]
    return dem.copy(data=vals), int(missing.sum())

Nearest-valid filling is the simplest approach and adequate for small voids. Large voids — radar shadow in mountains, water bodies in some products — deserve interpolation from the surrounding surface or a fill from a secondary DEM, both covered in filling voids in DEM tiles.

Step 3 — Derive slope and aspect at native resolution

import numpy as np


def slope_aspect(dem: np.ndarray, res: float) -> tuple[np.ndarray, np.ndarray]:
    gy, gx = np.gradient(dem.astype("float64"), res)
    slope = np.degrees(np.arctan(np.hypot(gx, gy)))
    aspect = (np.degrees(np.arctan2(-gx, gy)) + 360) % 360      # 0 = north, clockwise
    return slope.astype("float32"), aspect.astype("float32")

The sign conventions in the aspect formula are the part most often got wrong. Row index increases southward in a north-up raster, so the northward gradient is the negative of gy; with the arrangement above, aspect is measured clockwise from north, which is what every GIS and every downstream formula expects. The alternatives and their trade-offs are in computing slope and aspect from a DEM.

Step 4 — Align with the imagery

slope_da = dem.copy(data=slope)
slope_on_s2 = slope_da.rio.reproject_match(s2_red, resampling=1)   # bilinear

Resampling the derivative rather than the DEM is the rule that matters here; slope computed from an upsampled DEM is systematically too gentle.

Step 5 — Use the layers

Terrain layers serve three distinct purposes: as predictors in a model, where they capture environmental gradients, as described in adding terrain derivatives as predictors; as display context, where hillshade makes relief legible under other layers; and as inputs to correcting imagery for the uneven illumination slopes receive.

---

Why Terrain Illumination Matters for Imagery

In mountainous terrain this illumination effect can exceed the difference between land-cover classes. A classifier trained without accounting for it learns to split forest into a bright class and a dark class along the ridge line, and a change detector compares images whose shadows fell differently because the sun was lower. Topographic correction models the local illumination angle from slope, aspect and solar geometry and normalises the imagery against it, as covered in applying a terrain illumination correction.

---

Choosing a DEM Source

The right elevation model depends on the resolution of the imagery it will be paired with, whether bare-ground or surface heights are wanted, and what coverage exists for the area.

The surface-model caveat deserves attention in forested areas. A global 30 m DEM measures roughly the top of the canopy, so a forest edge appears as a cliff tens of metres high and the slope layer shows steep terrain where the ground is flat. For illumination correction that is arguably correct — the canopy is what reflects the light — but for anything about the ground itself, such as water flow or erosion, it is simply wrong, and a bare-earth product is needed.

Matching the DEM’s resolution to the imagery is the other practical rule. Pairing metre-scale lidar with 10 m imagery multiplies the processing without adding information the imagery can use; the slope is then aggregated back to 10 m anyway. Choose the DEM at or slightly finer than the imagery resolution, and derive at that resolution.

---

Terrain Products at Scale

Terrain layers are ideal candidates for computing once and reusing indefinitely. Elevation does not change on the timescales of most satellite analyses, so a slope, aspect and hillshade set derived for a country serves every project in that country for years. The efficient pattern is to derive them once at the native DEM resolution, store them as Cloud-Optimized GeoTIFFs with overviews, register them in the same catalog as the imagery, and read them by window from every downstream pipeline.

That reuse is where consistency pays off. If two projects derive slope separately, with different smoothing or different resampling, their results disagree for reasons unrelated to the question either is asking. A single reference set of terrain layers removes that variable. The authoring rules are those in writing and validating Cloud-Optimized GeoTIFFs, and the catalog registration follows the same approach as any other product.

Computing the derivatives over a large area is embarrassingly parallel with one caveat: neighbourhood operations such as gradients need a one-pixel overlap at tile edges, or seams appear along every boundary. Processing in windows with a small halo — the same pattern described for model inference in tiled inference with overlapping windows — removes them, and a one-pixel halo is all a 3×3 gradient needs.

---

Terrain and Other Pipelines

Terrain products connect to several other workflows on this site. Cloud-shadow detection uses the DEM to predict where shadows can and cannot fall, which sharpens masks in mountains. Mosaicking in steep terrain benefits from illumination-corrected inputs, so that seams between scenes acquired at different sun angles are less pronounced. And change detection in mountainous areas should either correct for illumination or restrict comparisons to anniversary dates with similar sun geometry, since a lower sun lengthens shadows and mimics loss of vegetation on shaded slopes.

In each case the DEM enters as supporting data rather than as the subject, which is why a consistent, well-documented reference set matters: every pipeline that uses it inherits its accuracy and its conventions.

---

Parameter Reference

---

Beyond Slope and Aspect

Slope and aspect are the workhorses, but several other derivatives earn their place in specific analyses, and all follow the same rules about projection and native-resolution derivation.

Curvature — the second derivative of elevation — separates convex ridges from concave hollows, which matters for soil moisture, erosion and where water accumulates. Profile curvature follows the steepest slope; plan curvature runs across it. Both are noisy on a global DEM and usually need light smoothing before they carry useful signal.

The topographic position index compares each cell’s elevation with the mean of its neighbourhood, classifying terrain into ridges, slopes, valleys and flats at a chosen scale. It is simple to compute with a focal mean and remarkably informative as a model feature, because it encodes landscape position in a single number.

Flow direction and accumulation, which trace where water would go, belong to hydrological analysis rather than imagery processing, and they are unusually sensitive to DEM quality: a single spurious pit can redirect a whole catchment. For that work a hydrologically conditioned DEM is essential, and a raw global surface model is not a suitable input.

Each of these can be computed with the same array tools used above, and each should be computed once, stored, and reused rather than re-derived per project. The feature-engineering view of several of them is in adding terrain derivatives as predictors.

---

Verification & Testing

import numpy as np

slope, aspect = slope_aspect(dem_filled.values, res=30.0)
assert np.nanmax(slope) < 90, "slope ≥ 90° means the DEM is not in metres"
assert 0 <= np.nanmin(aspect) and np.nanmax(aspect) < 360
assert np.isfinite(dem_filled.values).all(), "voids remain after filling"
print(f"median slope {np.nanmedian(slope):.1f}°, 99th pct {np.nanpercentile(slope, 99):.1f}°")

A quick visual check is often more informative than any assertion: render the hillshade and look at it. A correct hillshade looks like relief lit from the north-west; an inverted aspect convention makes valleys look like ridges, and an unfilled void shows as a flat grey hole. Both are obvious at a glance and subtle in the numbers.

---

Troubleshooting

Slope values in the thousands

The DEM is in geographic degrees. Reproject to a metre-based CRS before differentiating.

Hillshade looks inside-out

The aspect or the gradient sign is inverted. Check with a known north-facing slope.

Stripes or blocks in the slope layer

The DEM was resampled with nearest-neighbour, or derived after upsampling. Derive at native resolution with bilinear resampling.

Holes in slope and aspect

Voids in the DEM propagate to their neighbours through the gradient. Fill before deriving.

Terrain layers misalign with imagery by a pixel

They were resampled onto a slightly different grid. Use reproject_match against the imagery itself.

---

Frequently Asked Questions

Q: Which global DEM should I use? A modern 30 metre global DEM is the default choice for most satellite work at 10 to 30 metre resolution. Use a national lidar product where one exists and the analysis needs finer detail, and be aware that global DEMs measure the top of vegetation and buildings rather than bare ground.

Q: Why must a DEM be in a projected CRS before computing slope? Because slope divides a change in elevation, in metres, by a horizontal distance. In a geographic CRS that distance is in degrees, and the result is off by a factor of about a hundred thousand. Reproject to metres first.

Q: Is a DEM a digital terrain model or a surface model? Most global products are surface models: they include the heights of forest canopy and buildings. For hydrology or ground slope under forest, that bias matters and a terrain model, usually from lidar, is needed instead.

Q: How much vertical accuracy does a global DEM have? Typically a few metres in open terrain and considerably worse on steep slopes and under dense canopy. That is ample for slope classes and illumination correction at 30 metres, and inadequate for anything depending on sub-metre relief such as drainage in flat floodplains.

Q: Can terrain layers be computed on the fly instead of stored? They can, and for a small area the cost is negligible. For anything reused across projects or computed over large extents, storing them once is both cheaper and more consistent, since every consumer then sees exactly the same slope and aspect rather than slightly different re-derivations.

---

Related

• Computing Slope and Aspect from a DEM — algorithms, conventions and edge handling.
• Generating Hillshade Rasters in Python — relief shading for display.
• Filling Voids in DEM Tiles — from nearest-valid to secondary-source fills.
• Applying a Terrain Illumination Correction — normalising imagery for slope lighting.
• Adding Terrain Derivatives as Predictors — terrain as model features.