Generating Multidirectional Hillshade with gdaldem

Run gdaldem hillshade with -multidirectional when your terrain contains ridges running parallel to the conventional north-west light, because a single light source renders both flanks of such a ridge identically and it disappears from the map entirely.

Core Algorithm and Workflow

Single-source hillshade computes one number per cell: the cosine of the angle between the surface normal and a vector toward the light. A ridge running north-west to south-east under a north-west light presents both flanks at the same angle to that light, so both receive the same value and the ridge renders as a flat band. The terrain is in the data; it is simply not in the picture.

Multidirectional shading, as implemented in GDAL following Mark’s 1992 method, replaces the single source with four at 225, 270, 315 and 360 degrees. Each cell’s final value is a weighted combination, with weights derived from the cell’s own aspect so that the light most nearly perpendicular to the slope contributes most. The result is a shading in which no orientation is systematically invisible.

The workflow is otherwise unchanged from single-source shading:

  1. Validate the DEM’s CRS, vertical unit and nodata value.
  2. Compute the z-factor from the CRS, exactly as before — the lighting model does not affect the gradient.
  3. Invoke gdaldem hillshade with -multidirectional, omitting -az because it is ignored.
  4. Compare against a single-source render on terrain that runs north-west to south-east.
The ridge a single light source cannot show Two panels of the same terrain. Under a single north-west source, a ridge running north-west to south-east presents both flanks at the same angle to the light, so both render at the same value and the ridge is invisible. Under four sources at 225, 270, 315 and 360 degrees, the flanks receive different weighted illumination and the ridge is clearly defined. single source at 315° both flanks identical — the ridge is not visible four sources, aspect-weighted flanks separate — the ridge reads correctly Neither panel has more data than the other. The first simply cannot express what it contains.
This is not a subtle difference in tone. Under the wrong light angle an entire ridge system is absent from the map, and nothing in the pipeline reports it.

Production-Ready Python Implementation

import math
import subprocess
import rasterio


def multidirectional_hillshade(dem_path: str, out_path: str,
                               vertical_unit_metres: float = 1.0) -> str:
    """Render multidirectional shaded relief from a validated DEM.

    Returns the output path. Raises if the DEM lacks the metadata every
    terrain derivative depends on.
    """
    with rasterio.open(dem_path) as src:
        if src.crs is None:
            raise ValueError(f"{dem_path}: no CRS, so the z-factor is undefined")
        if src.nodata is None:
            raise ValueError(f"{dem_path}: no nodata value; coastlines will glow")
        geographic = src.crs.is_geographic
        centre_lat = (src.bounds.bottom + src.bounds.top) / 2

    if geographic:
        cos_lat = math.cos(math.radians(centre_lat))
        if abs(cos_lat) < 1e-6:
            raise ValueError("extent reaches a pole; reproject before shading")
        z = vertical_unit_metres / (111320.0 * cos_lat)
    else:
        z = vertical_unit_metres

    subprocess.run(
        [
            "gdaldem", "hillshade", dem_path, out_path,
            "-multidirectional",        # supersedes -az; do not pass both
            "-z", f"{z:.10g}",
            "-compute_edges",           # no nodata hairline at tile borders
            "-co", "COMPRESS=DEFLATE",
            "-co", "TILED=YES",
        ],
        check=True,
    )
    return out_path

Passing -az alongside -multidirectional is accepted and silently ignored, which is worth knowing because a script that sets an azimuth and sees no change will otherwise be debugged in the wrong place. Leaving it out documents the behaviour.

-compute_edges is not optional in a tiled workflow. Without it the outermost pixel row and column are nodata, and mosaicking tiles produces a visible one-pixel grid across the whole product.

Performance Tuning and Cartographic Best Practices

  • Derive at the delivery resolution. Shading a one-metre lidar DEM for a 1:50 000 sheet computes 625 samples for every output pixel. Resample the DEM close to the output pixel size first; the shading is indistinguishable and the run is two orders of magnitude cheaper.
  • Tile with a two-pixel overlap. The gradient estimator reads each cell’s eight neighbours, so a tile boundary without overlap produces a band of wrong values along every internal edge.
  • Stretch afterwards, not by exaggeration. Multidirectional output is flatter by construction. If it reads as washed out, apply a contrast stretch to the shading raster; raising the z-factor to compensate misrepresents the terrain and reintroduces terracing on gentle slopes.
  • Blend with a single-source render when both are wanted. A weighted average of roughly 60 per cent multidirectional and 40 per cent single-source keeps most of the ridge recovery while restoring directional drama. Average in linear light — averaging the 8-bit outputs biases the result toward the darker input.
  • Cache the derivative. Shading depends only on the DEM and the parameters, none of which vary per sheet, so it belongs in a build step rather than inside the per-sheet render loop.
What averaging four sources does to the tonal range Three histograms of shading values. Single-source shading spreads across the full zero to 255 range with peaks near both ends. Multidirectional shading is compressed into the middle third, which reads as flat. The stretched multidirectional histogram re-expands to the full range while keeping the ridge information the averaging recovered. single source full range, strong separation multidirectional compressed — reads as flat stretched range restored, ridges kept Stretch the shading raster; do not raise the z-factor to compensate for flatness. Raising the exaggeration changes what the map says about the terrain. Stretching changes only how the same shading is mapped onto the output range, which is a presentation decision.
The flatness is a property of the averaging, not of the terrain, so it should be corrected where it arises — in the tonal mapping, not in the geometry.
The four sources, and how a cell's aspect weights them A compass with four light sources marked at 225, 270, 315 and 360 degrees. Beside it, three example cells with different aspects show the weight each source contributes. A cell facing north-west takes most of its light from the 315 degree source. A cell facing south-west takes most from 225. A cell facing north takes most from 360, which is the source a conventional single-light render omits entirely. 360° 315° 270° 225° four fixed sources cell aspect dominant source weight north-west 315° 0.54 south-west 225° 0.51 north 360° 0.58 The third row is the case a single 315° light cannot express, because that cell faces the source edge-on and returns a mid value. Weights sum to one per cell, so the mode changes the distribution of light, not its total.
The weighting is per cell and derived from aspect, which is why the mode recovers detail rather than simply brightening everything by an average.

Integration and Next Steps

The shading raster produced here is an input to compositing rather than a finished map. Place it underneath the thematic layer with a multiply blend so hue is preserved, as described under Terrain and Hillshade Automation, and re-run contrast validation on the composite because a palette validated against white can fail against its own shaded form.

In a tile pipeline the derivative belongs to the build, not the request path: shade once at each zoom band’s working resolution, store the results, and let tile seeding distribute them. In an atlas pipeline it belongs to the resolve-once stage described in Atlas and Map Series Automation, so that every sheet composites against identical terrain.

Frequently Asked Questions

Does -multidirectional still respect the azimuth argument?

No. The mode uses four fixed sources at 225, 270, 315 and 360 degrees, weighted per cell by aspect, and any -az value is accepted and ignored. Scripts that set both and then investigate why the azimuth has no effect are debugging the wrong layer. Omit the flag in multidirectional mode so the command states what it actually does.

Why is multidirectional output lower contrast than single-source?

Averaging four directions necessarily pulls values toward the middle: a slope fully lit from one direction is partly shadowed from another, and the weighted combination lands between. That flattening is the price of the mode’s benefit. Where it matters, apply a contrast stretch to the shading raster afterwards. Raising the z-factor produces a superficially similar increase in contrast while also exaggerating the terrain and amplifying any quantisation terracing in the source DEM.

Can I combine multidirectional and single-source shading?

Yes, and it is a reasonable default for large-format products: render both and take a weighted average, around 60 per cent multidirectional to 40 per cent single-source. The blend keeps most of the ridge recovery while restoring some directional drama. Do the averaging in linear light rather than on the 8-bit outputs, because the transfer curve biases a naive average toward the darker of the two inputs.

Is multidirectional shading appropriate for a physical relief map?

Usually not. A small-scale physical map communicates the shape and arrangement of major landforms, and strong directional lighting is what makes a mountain range legible at a glance. Multidirectional shading earns its place on large-scale maps covering terrain that runs in many directions within one sheet, where the risk of losing a whole ridge system to an unlucky light angle is real and the loss of drama is affordable.

Back to Terrain and Hillshade Automation