Deriving Contours and Hillshade with GDAL

A DEM is rarely the deliverable. What the client opens is a contour set they will draft against, a hillshade they will place under an orthomosaic, or a slope raster feeding a stability analysis. All three are derived from the same elevation grid by short, well-defined operations — and all three amplify whatever noise that grid contains, because every one of them is a derivative.

This page covers the derivation chain: how much to smooth and where, how to choose a contour interval the data actually supports, why hillshade is illustration while slope is measurement, and how to keep the whole set consistent so the contours a client drafts against agree with the surface they are drawn from.

The DEM this consumes is the one produced in generating DSM and DTM from point clouds with PDAL, and the export conventions are those of exporting Cloud-Optimized GeoTIFFs with rasterio.

Audience and prerequisites. Python 3.10+, GDAL ≥ 3.6 with its Python bindings or the CLI on PATH, and a DEM whose vertical datum you can state.

# Prerequisites

# Conceptual architecture

Every product on this page is a spatial derivative of elevation, and derivatives amplify high-frequency content. A DEM with 3 cm of per-cell noise at a 5 cm grid has an implied slope noise of over 30° between adjacent cells — which is why an unsmoothed hillshade of a photogrammetric DEM looks like sandpaper and why contours drawn from one wander.

Figure 1 — The derivation tree. Smoothing belongs above the branch, once, so the four products describe the same surface rather than four slightly different ones.

# Step 1: Smooth deliberately, and record how much

The smoothing kernel is a survey decision, not an aesthetic one: it sets the smallest terrain feature that survives into the deliverable. A Gaussian with a standard deviation of one cell removes most per-cell noise and preserves features a few cells across; three cells produces a visibly cleaner product and erases a kerb.

import numpy as np
import rasterio
from scipy.ndimage import gaussian_filter


def smooth_dem(src_path: str, dst_path: str, sigma_cells: float = 1.0) -> None:
    """Gaussian-smooth a DEM, preserving NoData rather than smearing it."""
    with rasterio.open(src_path) as src:
        band = src.read(1, masked=True)
        profile = src.profile

    filled = band.filled(np.nan)
    valid = np.isfinite(filled)

    # Normalised convolution: smooth values and the mask, then divide, so
    # cells beside a void are not pulled toward zero by the void itself.
    num = gaussian_filter(np.where(valid, filled, 0.0), sigma_cells)
    den = gaussian_filter(valid.astype(float), sigma_cells)
    out = np.where(den > 1e-6, num / np.maximum(den, 1e-6), profile["nodata"])
    out[~valid] = profile["nodata"]           # voids stay voids

    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(out.astype(profile["dtype"]), 1)
        dst.update_tags(SMOOTH_SIGMA_CELLS=str(sigma_cells))

The normalised convolution matters at every void edge. A plain Gaussian treats a NoData cell as a zero elevation, which drags the cells around every void several metres downward and produces a ring of false depression that contours will faithfully draw. Recording the sigma in the raster’s tags means a later comparison between two epochs can check that both were smoothed the same way.

# Step 2: Choose an interval the data supports

A contour interval finer than the DEM’s vertical accuracy draws lines the data cannot justify: the contours are then mapping noise, and they wander in a way that looks like terrain detail. The conventional rule is an interval of at least twice the vertical RMSE, and at least three times where the client will draft against them.

Figure 2 — The interval follows from the measured vertical accuracy. Quoting a finer interval than the band allows is a claim about the data that the accuracy report contradicts.

import subprocess


def contours(dem_path: str, out_path: str, interval_m: float,
             index_every: int = 5, layer: str = "contours") -> None:
    """Generate contours with an index attribute for cartographic weighting."""
    subprocess.run(
        ["gdal_contour", "-a", "elev", "-i", str(interval_m),
         "-f", "GPKG", "-nln", layer, "-amin", "elev_min", "-amax", "elev_max",
         dem_path, out_path],
        check=True,
    )
    # Index contours are every Nth line; tag them so the map can weight them.
    subprocess.run(
        ["ogrinfo", out_path, "-dialect", "SQLite", "-sql",
         f"UPDATE {layer} SET is_index = "
         f"(CAST(elev / {interval_m} AS INTEGER) % {index_every} = 0)"],
        check=True,
    )

# Step 3: Hillshade is illustration; slope is measurement

gdaldem hillshade combines slope and aspect with an assumed light direction to produce an 8-bit image. It is invaluable for reading terrain and it is not a measurement: change the azimuth and the same landscape looks different, and features aligned with the light direction disappear entirely.

That last property is worth stating to clients. A hillshade lit from the north-west — the conventional default, because human perception reads north-west lighting as convex — makes north-west-trending features nearly invisible. On a site whose main structures run that way, a second hillshade at a different azimuth, or a multi-directional hillshade, is not decoration but the difference between seeing the site and not.

import subprocess


def hillshade(dem_path: str, out_path: str, azimuth: float = 315.0,
              altitude: float = 45.0, z_factor: float = 1.0,
              multidirectional: bool = False) -> None:
    cmd = ["gdaldem", "hillshade", dem_path, out_path,
           "-az", str(azimuth), "-alt", str(altitude), "-z", str(z_factor),
           "-compute_edges", "-co", "COMPRESS=DEFLATE"]
    if multidirectional:
        cmd.append("-multidirectional")       # ignores -az; lights from four
    subprocess.run(cmd, check=True)


def slope_and_aspect(dem_path: str, slope_path: str, aspect_path: str) -> None:
    """Slope in degrees and aspect in degrees clockwise from north."""
    for mode, out in (("slope", slope_path), ("aspect", aspect_path)):
        subprocess.run(
            ["gdaldem", mode, dem_path, out, "-compute_edges",
             "-co", "COMPRESS=DEFLATE", "-co", "PREDICTOR=3"],
            check=True,
        )

-compute_edges deserves a mention: without it, gdaldem leaves a one-cell NoData border on every output because the 3 × 3 window has no neighbours there. Across a tiled workflow that border becomes a visible grid of missing lines, and it is one flag away.

The z_factor is only ever 1.0 when horizontal and vertical units match. On a DEM in a geographic CRS — degrees horizontally, metres vertically — slope computed with the default is meaningless by roughly a factor of 100,000. Reproject to a projected CRS before deriving anything.

# Contours across a void, and other places the lines lie

A contour generator has no notion of confidence. Given a filled DEM it will draw through interpolated ground with exactly the same weight it gives measured ground, and the resulting lines are indistinguishable on paper. On a site with a large water body or a deep occlusion behind a building, that means a client drafts against contours that were invented by a fill algorithm.

Three habits keep this honest. Contour the unfilled surface where the deliverable permits it, so voids appear as gaps and the reader can see them. Clip the contour set to the fill mask described in fixing holes and voids in drone DSMs when the deliverable must be continuous, and ship the mask alongside. And never contour across water at all: a water surface returns no reliable elevation, the fill invents one, and contours over a lake are pure fiction that looks like bathymetry.

The same argument applies at the survey boundary. Contours generated to the raster’s edge continue into the region where the block had control on one side only, which is where the reconstruction is weakest. Clipping the contour set to the area of interest, inset from the flown extent rather than coincident with it, discards exactly the lines least supported by the data.

# Parameter deep-dive

# Cartographic weight is part of the deliverable

A contour set with every line the same weight is technically complete and practically hard to read. The convention that makes contours legible is hierarchy: index contours every fourth or fifth interval, drawn heavier and carrying an elevation label, with intermediate contours lighter and unlabelled. That is a property of the data — an attribute on each feature — not of the styling, because a client who opens the layer in their own software will not inherit your symbology.

Two further attributes are worth writing at generation time. The source surface, so a DSM contour set is never confused with a DTM one. And the smoothing sigma, so a later revision can be produced identically rather than approximately. Both are single fields, both are known at the moment of generation, and both are impossible to recover afterwards from the geometry alone.

Labels themselves are best left to the client’s cartography. Placing them at generation time bakes in a scale assumption — a label spacing that reads well at 1:500 is unreadable at 1:2000 — and the elevation attribute is all any mapping package needs to place its own.

# Verification and output inspection

The check that matters is consistency: contours must lie on the surface they claim to describe. Sampling the DEM along each contour vertex and comparing against the contour’s own elevation attribute catches a whole class of mistakes — contours generated from a different (unsmoothed, or differently-projected) DEM than the one shipped.

import fiona
import numpy as np
import rasterio


def assert_contours_match_dem(gpkg: str, dem_path: str, tol_m: float = 0.02,
                              sample_every: int = 10) -> None:
    """Every contour vertex should sit at its stated elevation on the DEM."""
    worst = 0.0
    with rasterio.open(dem_path) as dem, fiona.open(gpkg) as src:
        for feat in src:
            elev = float(feat["properties"]["elev"])
            coords = feat["geometry"]["coordinates"]
            pts = [tuple(p[:2]) for p in coords[::sample_every]]
            if not pts:
                continue
            for val in dem.sample(pts, indexes=1):
                if np.isfinite(val[0]):
                    worst = max(worst, abs(float(val[0]) - elev))
    assert worst <= tol_m, (
        f"contour vertices depart from the DEM by up to {worst:.3f} m — "
        "the contours were generated from a different raster")

A tolerance of half a cell’s worth of interpolation is generous; anything larger means the two products disagree, which a client will discover by drafting against the contours and finding the surface elsewhere.

# Troubleshooting

Contours are dense, wandering squiggles rather than smooth lines. The interval is finer than the DEM’s vertical noise. Either smooth more, or widen the interval to at least twice the checkpoint RMSE. Post-smoothing the lines instead is the wrong fix — it makes them look right while leaving them off the surface, which the verification above will catch.

A ring of concentric contours surrounds every void. The smoothing treated NoData as zero elevation and pulled the surrounding cells down. Use the normalised convolution above, or mask before smoothing.

Slope values are absurd — thousands of degrees. The DEM is in a geographic CRS, so horizontal units are degrees and vertical are metres. Reproject to a projected CRS; z_factor is a workaround that is only correct at one latitude.

Hillshade has a grid of thin dark lines. Each tile was shaded independently without -compute_edges, leaving a NoData border. Add the flag, or shade the mosaic rather than the tiles.

The client says the terrain “looks flat” in one direction. Features aligned with the light azimuth are invisible by construction. Supply a second azimuth or use -multidirectional; this is a property of hillshading, not a defect in the DEM.

Contours disagree with the DSM but agree with the DTM, or vice versa. Two surfaces are in play and the contour set was generated from the other one. Name the source surface in the contour layer’s metadata — a contour file that does not say which surface it came from is ambiguous the moment both are delivered.

Ship the smoothed DEM alongside the contours, not only the raw one. It is the surface the lines were actually drawn from, and a client who wants to reproduce or extend the contour set needs it rather than the noisier original.

# Related

• Generating DSM and DTM from point clouds with PDAL
• Generating smooth contours without staircase artifacts
• Computing slope and aspect rasters in Python
• Exporting Cloud-Optimized GeoTIFFs with rasterio
• Fixing holes and voids in drone DSMs

← DEM/DSM Generation & Raster Export Automation

Figure 3 — The systematic blindness of directional shading. Worth knowing before a client reports that a road or a pipeline trench “is not in the data”.