Terrain and Hillshade Automation

Shaded relief is the one map element where a wrong parameter produces something that still looks like a map. A hillshade computed with the wrong vertical exaggeration is not obviously broken — it simply describes terrain that does not exist, with mountains too flat or too jagged, and nothing in the pipeline raises an error. This page covers the automation of terrain derivatives end to end: establishing the DEM contract, computing the z-factor correctly, generating shaded relief that reads as terrain rather than as texture, and compositing it under thematic colour without destroying the data the map exists to show.

Everything here assumes a batch context. A cartographer working interactively can tune an azimuth until the map looks right; a pipeline rendering four hundred atlas sheets cannot, and needs parameters derived from the data rather than chosen by eye.

Prerequisites and Environment Configuration

python==3.11
GDAL==3.8.4
rasterio==1.3.9
numpy==1.26.4
scipy==1.12.0
matplotlib==3.8.4

Three properties of the source DEM must be known before any derivative is computed, and none of them can be safely assumed:

• The horizontal CRS and its units. A DEM in EPSG:4326 has degrees on both horizontal axes; a DEM in a UTM zone has metres. The slope and hillshade algorithms compute a gradient, which is a ratio of vertical change to horizontal distance, so the units of both must be reconciled before the ratio means anything.
• The vertical unit. Metres and feet are both common, and nothing in a GeoTIFF forces the vertical unit to be recorded. A DEM in feet processed as metres yields relief exaggerated by 3.28, which reads as spectacular alpine terrain in a region that has none.
• The nodata value. Ocean cells and voids are typically encoded as a sentinel such as −32768 or −9999. If nodata is not declared, the gradient across the boundary between land and sentinel is enormous, and a bright rim appears along every coastline.

Assert all three at load time and fail loudly when any is missing. A pipeline that guesses these values produces plausible output, which is precisely the failure mode that survives review.

For projection choice, the same reasoning as in Projection Selection Algorithms applies with one addition: terrain derivatives are computed per pixel from neighbouring pixels, so a projection whose scale varies strongly across the extent produces relief whose apparent steepness varies with position. For a national terrain product, reproject the DEM to an equal-area or conformal projection appropriate to the extent before deriving anything.

Conceptual Foundation: What the Shading Algorithm Actually Computes

Hillshade is a single number per pixel: the cosine of the angle between the surface normal and a vector toward the light source, clamped to zero and scaled to the output range. The surface normal comes from the local gradient, estimated over the eight neighbours of each cell using Horn’s method. Everything cartographers argue about — azimuth, altitude, exaggeration, multidirectional blending — is a modification of that one calculation.

The gradient is where units enter. Horn’s estimator divides the elevation difference between neighbours by the horizontal distance between them, and both quantities have to be in the same unit for the resulting slope to be an angle. The z-factor is the conversion constant that makes that true:

import math


def z_factor(crs_is_geographic: bool, centre_latitude: float,
             vertical_unit_metres: float = 1.0) -> float:
    """Vertical exaggeration that reconciles vertical and horizontal units.

    For a projected CRS in metres with elevation in metres the factor is 1.0.
    For a geographic CRS the horizontal unit is degrees, whose ground length
    varies with latitude, so the factor depends on where the raster sits.
    """
    if not crs_is_geographic:
        return vertical_unit_metres

    metres_per_degree = 111320.0 * math.cos(math.radians(centre_latitude))
    if metres_per_degree <= 0:
        raise ValueError("degenerate latitude for z-factor computation")
    return vertical_unit_metres / metres_per_degree

Two consequences follow that catch out most first implementations. First, the factor for a geographic DEM is a very small number, around 1.1e-5 at the equator — not the 0.00001 or 111120 that appear in various forum answers, both of which are the same idea with a different sign of error. Second, because the factor depends on latitude, a single value applied across a continental DEM exaggerates northern relief relative to southern. For a map series spanning more than about ten degrees of latitude, compute the factor per sheet.

Step-by-Step Implementation

Step 1: Validate the DEM before deriving anything

import rasterio


def load_dem(path: str) -> dict:
    """Read a DEM and assert the three properties every derivative depends on."""
    with rasterio.open(path) as src:
        if src.nodata is None:
            raise ValueError(f"{path}: no nodata value declared — coastlines will glow")
        if src.crs is None:
            raise ValueError(f"{path}: no CRS — the z-factor cannot be computed")

        bounds = src.bounds
        centre_lat = (bounds.bottom + bounds.top) / 2 if src.crs.is_geographic else 0.0
        return {
            "path": path,
            "crs": src.crs,
            "is_geographic": src.crs.is_geographic,
            "centre_latitude": centre_lat,
            "nodata": src.nodata,
            "res": src.res,
        }

The vertical unit cannot be read from the file in the general case, so it belongs in the dataset’s own manifest alongside its licence and vintage. Treat it as required metadata rather than something to infer.

Step 2: Generate the shaded relief

The gdaldem utility implements Horn’s method and is faster than any pure-Python equivalent. Drive it through subprocess rather than reimplementing the gradient:

import subprocess


def hillshade(dem: dict, out_path: str, azimuth: float = 315.0,
              altitude: float = 45.0, multidirectional: bool = False) -> str:
    """Render shaded relief from a validated DEM."""
    zf = z_factor(dem["is_geographic"], dem["centre_latitude"])

    cmd = ["gdaldem", "hillshade", dem["path"], out_path,
           "-z", f"{zf:.10g}", "-compute_edges",
           "-co", "COMPRESS=DEFLATE", "-co", "TILED=YES"]

    if multidirectional:
        cmd.append("-multidirectional")   # azimuth is ignored in this mode
    else:
        cmd += ["-az", str(azimuth), "-alt", str(altitude)]

    subprocess.run(cmd, check=True)
    return out_path

Two flags matter more than they appear to. -compute_edges fills the one-pixel border that would otherwise be nodata, which prevents a hairline seam at every tile boundary in a mosaicked product. -multidirectional replaces the single light source with a weighted combination of four, and its effect is specific: terrain whose ridges run parallel to the light direction is invisible under single-source illumination, because both flanks receive identical light. Multidirectional shading recovers those ridges at the cost of slightly softer overall contrast.

Step 3: Choose the azimuth deliberately

The conventional north-west azimuth of 315 degrees is not arbitrary. Human depth perception assumes light from above, and a light source in the southern half of the compass inverts the reading: ridges appear as valleys and craters as domes. This is the hollow-mountain illusion, and it is strong enough that readers do not merely find the map ambiguous — they confidently read it backwards.

Keep the primary azimuth between about 270 and 360 degrees. When a particular ridge system runs north-west to south-east and disappears under that lighting, use multidirectional mode rather than rotating the light southward.

Step 4: Composite under the thematic layer

This is where most terrain automation goes wrong. Placing a grey hillshade over a thematic map at partial alpha mixes grey into every class, desaturating the palette and pulling the whole map toward the middle of the value range. The legend swatches, drawn without the overlay, then no longer match anything on the map.

The correct composite puts relief underneath and uses a blend that modulates value while preserving hue:

import numpy as np


def multiply_blend(thematic_rgb: np.ndarray, shade: np.ndarray,
                   strength: float = 0.55) -> np.ndarray:
    """Composite shaded relief under a thematic layer, preserving hue.

    thematic_rgb : float array (h, w, 3) in 0..1
    shade        : float array (h, w) in 0..1, 0.5 = neutral illumination
    """
    if thematic_rgb.shape[:2] != shade.shape:
        raise ValueError("shade and thematic layer must share a grid")

    # Pull the shade toward neutral by `strength`, so 0 leaves colour untouched.
    modulated = 1.0 - strength * (1.0 - shade[..., None])
    return np.clip(thematic_rgb * modulated, 0.0, 1.0)

Because multiply only ever darkens, the palette’s lightest class stays recognisable and every class keeps its hue. A strength around 0.5 gives visible relief without compromising the contrast ratios validated under Accessibility Sync in Cartography — although the validation must be re-run on the composited output, because a class that passed against white may not pass against its own shaded form.

Complete Working Code Example

import math
import subprocess
import numpy as np
import rasterio


def terrain_basemap(dem_path: str, thematic_rgb: np.ndarray,
                    out_path: str, strength: float = 0.55) -> str:
    """Derive shaded relief from a DEM and composite it under a thematic layer.

    thematic_rgb must already be on the DEM's grid; resample it beforehand
    rather than here, so the resampling method is an explicit decision.
    """
    with rasterio.open(dem_path) as src:
        if src.nodata is None or src.crs is None:
            raise ValueError(f"{dem_path}: DEM is missing nodata or CRS")
        profile = src.profile
        centre_lat = (src.bounds.bottom + src.bounds.top) / 2
        geographic = src.crs.is_geographic

    zf = 1.0
    if geographic:
        zf = 1.0 / (111320.0 * math.cos(math.radians(centre_lat)))

    shade_path = out_path + ".shade.tif"
    subprocess.run(
        ["gdaldem", "hillshade", dem_path, shade_path,
         "-z", f"{zf:.10g}", "-multidirectional", "-compute_edges",
         "-co", "COMPRESS=DEFLATE", "-co", "TILED=YES"],
        check=True,
    )

    with rasterio.open(shade_path) as shd:
        shade = shd.read(1).astype("float32") / 255.0

    if shade.shape != thematic_rgb.shape[:2]:
        raise ValueError("thematic layer is not on the DEM grid — resample first")

    modulated = 1.0 - strength * (1.0 - shade[..., None])
    composite = np.clip(thematic_rgb * modulated, 0.0, 1.0)

    profile.update(count=3, dtype="uint8", nodata=None,
                   compress="deflate", tiled=True)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write((composite * 255).astype("uint8").transpose(2, 0, 1))

    return out_path

The function deliberately refuses to resample. Choosing between nearest-neighbour and bilinear for a thematic raster is a cartographic decision — nearest preserves class boundaries, bilinear invents intermediate classes that appear in no legend — and hiding it inside a compositing helper is how a categorical layer silently acquires values that do not exist.

Performance Optimization Patterns

Derive once, reuse across sheets. Hillshade depends only on the DEM, the z-factor and the light parameters. None of those vary per atlas sheet, so the derivative belongs in a build step that runs once, not inside the per-sheet loop. This is the same hoisting argument made about font metrics and spatial indexes in the performance section of the fundamentals overview, and terrain is usually the largest instance of it.

Compute at the delivery resolution, not the source resolution. A one-metre lidar DEM composited into a 1:50 000 sheet contributes roughly one useful sample per 25 output pixels; deriving hillshade at full source resolution and downsampling afterwards costs 625 times the memory for output that is very slightly softer. Resample the DEM to approximately the output pixel size first, then derive.

Tile large derivations with an overlap. A national DEM will not fit in memory. Process it in tiles with an overlap of at least two pixels on every side, because the gradient estimator reads the eight neighbours of each cell; without the overlap, a one-pixel band of wrong values appears along every internal tile edge and mosaics into a visible grid.

Use the compressed, tiled GeoTIFF options. COMPRESS=DEFLATE with TILED=YES typically reduces a hillshade to a quarter of its uncompressed size and makes windowed reads cheap, which matters when the compositing step reads the raster back one block at a time.

Common Pitfalls and Debugging

Glowing coastlines. An undeclared nodata value makes the gradient between a −32768 sentinel and a 3-metre coastal cell enormous, and the shading algorithm renders that as a bright rim. Declare nodata; if the source genuinely lacks one, mask the sentinel explicitly before deriving.

Terracing on gentle slopes. An integer DEM has a one-unit vertical quantum. Where the terrain rises by less than one unit per pixel, the derived slope snaps between discrete values and the hillshade shows contour-like steps. Convert to float and apply a light Gaussian smoothing — a sigma of about one pixel — before deriving. Raising the z-factor makes this worse, not better, which is why it is often misdiagnosed.

Relief that varies in strength across a large map. This is the latitude-dependent z-factor problem. If a single factor was used across a continental extent, northern terrain will read as steeper than southern terrain of identical slope. Recompute per tile or reproject to a projected CRS first.

The thematic layer looks flat after compositing. Almost always a sign that the shade raster was applied over rather than under, or that alpha was used instead of multiply. Check the composite order first; it is a one-line fix that is easy to overlook because the output still looks like a reasonable map.

Contrast failures that appear only on shaded areas. A palette validated against a white background can fail against its own shaded form, because multiply darkens the fill while the label colour stays fixed. Run the contrast validation on the composited raster, sampling in both the brightest and darkest shaded regions.

Conclusion

Terrain shading is a deterministic transform with a small number of parameters, and every one of those parameters can be derived from the data rather than chosen by eye: the z-factor from the CRS and vertical unit, the azimuth from a perceptual constraint, the blend mode from whether a thematic layer is present, the resolution from the delivery scale. Automating it well means resolving each of those from measurable inputs and asserting the DEM contract loudly enough that a missing nodata value fails the build rather than the map. From here, the compositing decisions connect directly to Color Theory for GIS, and the per-sheet reuse pattern to the batch machinery in Batch Queue Orchestration.

Related

• Generating Multidirectional Hillshade with gdaldem — the four-source illumination model and when it beats a single light.
• Blending Hillshade with Landcover Without Muddy Colour — multiply, soft light and the contrast re-validation each requires.
• Debugging Terracing Artefacts in Derived Terrain Rasters — quantisation, smoothing kernels and how to tell terracing from real structure.
• Color Theory for GIS — why a multiply blend preserves the correspondence between map and legend.

Back to Automated Cartographic Design Fundamentals