Multidirectional Hillshade for LiDAR DTMs

TL;DR: gdaldem hillshade -multidirectional dtm.tif hs_multi.tif combines illumination from several azimuths so linear features at any orientation are visible, unlike a single 315° light that hides features running parallel to it. For interpretation work, compute your own weighted multi-azimuth hillshade in NumPy, then blend it with a slope layer to separate relief from shading artefacts.

# Context and Motivation

This guide is part of Hillshade, Slope and Aspect. A standard hillshade lights the terrain from one direction, conventionally the north-west at 315°. It is intuitive, but it has a blind spot: a ditch, bank, field boundary or fault scarp running north-west to south-east is lit along its length and casts no shadow, so it nearly disappears. On LiDAR DTMs, which are used precisely to find subtle linear features — relict field systems, hollow ways, drainage lines, landslide scarps — that blind spot matters. Multidirectional hillshading averages illumination from several directions so no orientation is favoured, while keeping the familiar top-left-lit look.

# Prerequisites and Assumptions

• A bare-earth DTM GeoTIFF in a projected CRS with metres for both horizontal and vertical units.
• GDAL 2.2+ for -multidirectional (the option has been available for years in gdaldem).
• For the custom version: NumPy, rasterio.

# Step-by-Step Implementation

# Step 1 — Standard multidirectional hillshade

gdaldem hillshade -multidirectional combines hillshades from four azimuths — 225°, 270°, 315° and 360° — weighted by the local aspect, following the method of Mark (1992) as implemented in GDAL.

# Step 2 — Combined shading for low relief

-combined blends slope and oblique shading, which brings out subtle relief on flat ground without the overall grey cast of a standard hillshade.

# Step 3 — Custom azimuth sets

For interpretation, compute hillshades at 8 or 16 azimuths and average them, optionally weighting a preferred direction.

# Step 4 — Vertical exaggeration

A z-factor of 2–3 emphasizes low relief; keep it at 1 for honest cartography.

# Step 5 — Blend with slope

Multiplying hillshade by a slope-derived layer (dark on steep slopes) separates true steepness from lighting effects.

# Complete Working Example

Command line:

gdaldem hillshade -multidirectional -z 1.0 -compute_edges -co COMPRESS=DEFLATE \
  mosaic/dtm_1m.tif out/hs_multi.tif
gdaldem hillshade -combined -z 2.0 -compute_edges mosaic/dtm_1m.tif out/hs_combined.tif
gdaldem hillshade -az 315 -alt 45 -compute_edges mosaic/dtm_1m.tif out/hs_315.tif

A custom 16-azimuth hillshade blended with slope in NumPy:

"""Custom multi-azimuth hillshade blended with slope."""
from __future__ import annotations

import numpy as np
import rasterio


def hillshade(z: np.ndarray, res: float, azimuth: float, altitude: float = 45.0,
              zfactor: float = 1.0) -> np.ndarray:
    dzdy, dzdx = np.gradient(z * zfactor, res)
    slope = np.arctan(np.hypot(dzdx, dzdy))
    aspect = np.arctan2(-dzdx, dzdy)                     # radians, clockwise from north
    az, alt = np.radians(azimuth), np.radians(altitude)
    hs = np.sin(alt) * np.cos(slope) + np.cos(alt) * np.sin(slope) * np.cos(az - aspect)
    return np.clip(hs, 0, 1)


def multi(z: np.ndarray, res: float, n: int = 16, zfactor: float = 1.0) -> np.ndarray:
    stack = [hillshade(z, res, az, zfactor=zfactor) for az in np.linspace(0, 360, n, endpoint=False)]
    return np.mean(stack, axis=0)


with rasterio.open("mosaic/dtm_1m.tif") as ds:
    z = ds.read(1, masked=True).filled(np.nan)
    profile, res = ds.profile, ds.res[0]

hs = multi(z, res, n=16, zfactor=2.0)
dzdy, dzdx = np.gradient(z, res)
slope_deg = np.degrees(np.arctan(np.hypot(dzdx, dzdy)))
slope_term = 1 - np.clip(slope_deg / 45.0, 0, 1) * 0.5      # darken steep slopes by up to half
out = (np.nan_to_num(hs * slope_term) * 255).astype("uint8")

profile.update(dtype="uint8", nodata=0, count=1, compress="deflate")
with rasterio.open("out/hs_multi16_slope.tif", "w", **profile) as o:
    o.write(out, 1)

Averaging many azimuths removes directional bias entirely, producing a flatter-looking image; that is ideal for detecting features, less so for maps meant to look natural. GDAL’s multidirectional mode keeps a dominant north-west light and is the better choice for cartography.

# Choosing a Rendering for the Job

Different audiences need different renderings of the same DTM, and it is worth producing more than one rather than arguing for a single “best” hillshade.

For general maps and reports, GDAL’s multidirectional hillshade at z-factor 1 reads naturally: light still appears to come from the upper left, which readers expect, and no terrain orientation is hidden. Overlay it semi-transparently on imagery or a colour-ramped elevation layer.

For feature detection and interpretation — archaeological survey, mapping landslide scarps, tracing drainage — an equal-weight multi-azimuth hillshade with mild exaggeration, blended with slope, is more useful. It looks flatter and less natural, but relief of every orientation carries equal contrast, which is what a searching eye needs. Local relief models and sky-view factor are further options in the same family, and interpreters often switch between several.

For very flat terrain, -combined or a hillshade with a low sun altitude (20–30°) brings out centimetre-scale relief that a standard rendering flattens into uniform grey.

Whatever the choice, record the parameters — azimuths, altitude, z-factor — with the output, because a feature that is visible under one rendering and not another is itself a finding.

# Key Parameter Table

# Verification

• Orientation test. Find a known linear feature parallel to 315° and confirm it is visible in the multidirectional output but faint in the single-azimuth one.
• No NoData border. With -compute_edges, edge cells are valid.
• Units. If the output is nearly uniformly white or black, horizontal and vertical units differ; see below.

# Gotchas and Edge Cases

Geographic rasters. A DTM in degrees with heights in metres needs a z-factor or -s 111120 (scale) so slopes are computed correctly; otherwise the hillshade is flat grey. Better, reproject to a projected CRS first.

Tile seams. Hillshading tiles separately leaves a one-pixel seam at every edge. Hillshade a VRT mosaic, or use -compute_edges on buffered tiles.

Interpretation bias. Every hillshade is a rendering choice. For archaeological or geomorphological interpretation, inspect several renderings — multidirectional, slope, local relief — rather than trusting one.

Output bit depth. Hillshades are usually 8-bit; computing in float and scaling once avoids banding when blending several layers.

# Frequently Asked Questions

What is a multidirectional hillshade?

A hillshade that combines illumination from several azimuths rather than one, so terrain features of every orientation cast visible shading. GDAL’s gdaldem supports it with the multidirectional option.

Why do some features disappear in a normal hillshade?

Features running parallel to the light direction are lit along their length and cast no shadow, so they show little contrast. With the default light from the north-west, north-west to south-east features are hardest to see.

What z-factor should I use?

One for faithful cartography in projected coordinates with metric heights. Two or three emphasizes low relief for interpretation. Geographic rasters need a scale conversion rather than an arbitrary z-factor.

Should I use -combined or -multidirectional?

Multidirectional for general terrain visualization with no directional blind spots; combined for very flat terrain, where blending slope shading brings out subtle relief. Try both on a sample area.

# Related

• Hillshade, Slope and Aspect — terrain derivatives overview
• Exporting Hillshade from a LiDAR DTM — the standard single-azimuth product
• Generating Slope and Aspect Rasters with gdaldem — the slope layer used in blending
• Computing Terrain Roughness and TPI — relief metrics without lighting
• Building a Seamless DTM Mosaic from Tiles — hillshading without seams