⚙️ How-To Terrain & elevation Maps & visualisation Terrain & Elevation · step 6 of 12 Updated

How to Make a Hillshade from a DEM in Python

Problem statement

A DEM plotted as a colour ramp is nearly unreadable. Elevation bands tell you how high things are and almost nothing about their shape — ridges, valleys, cliffs and gullies all vanish into a smooth gradient.

A hillshade solves it by simulating a light source. It is the single most effective thing you can do to make terrain legible, and it is one formula:

hillshade = 255 * (sin(altitude) * cos(slope)
                   + cos(altitude) * sin(slope) * cos(azimuth - aspect))

The formula is easy. What goes wrong is everything around it: the slope and aspect it depends on, the illumination angles, and the fact that a hillshade of a DSM shows every tree and building rather than the terrain.

Quick answer 

import math

import numpy as np
import rasterio


def hillshade(dem, cell_x, cell_y, *, azimuth=315, altitude=45):
    p = np.pad(dem, 1, mode="edge")
    a, b, c = p[:-2, :-2], p[:-2, 1:-1], p[:-2, 2:]
    d,    f = p[1:-1, :-2],               p[1:-1, 2:]
    g, h, i = p[2:, :-2],  p[2:, 1:-1],  p[2:, 2:]

    dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cell_x)
    dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cell_y)

    slope = np.arctan(np.hypot(dz_dx, dz_dy))
    aspect = np.arctan2(dz_dy, -dz_dx)

    az = math.radians(360.0 - azimuth + 90.0)
    alt = math.radians(altitude)

    shaded = (math.sin(alt) * np.cos(slope)
              + math.cos(alt) * np.sin(slope) * np.cos(az - aspect))
    return np.clip(shaded, 0, 1) * 255


hs = hillshade(dem, cell_x=27.9, cell_y=30.7)
print(f"mean {hs.mean():.1f}  sd {hs.std():.1f}  range {hs.min():.0f}-{hs.max():.0f}")
mean 159.0  sd 55.3  range 0-255

Two conventions and two parameters:

Parameter Default What it does
azimuth 315° (NW) where the light comes from
altitude 45° how high the sun sits

315° is not arbitrary. Lighting from the north-west is a cartographic convention, and departing from it triggers relief inversion in most viewers — hills read as hollows.

A hillside lit from the north-west at 45 degrees altitude, showing the illuminated and shadowed faces.
Two angles define the light. The 315° convention exists because human perception assumes light from above-left.

Step-by-step solution

1. Start from a DTM if you can

A hillshade of a DSM shows the canopy and the rooftops. Over a city that is a striking image of the buildings; over farmland it is a texture of hedgerows; over forest it is a rough surface that hides the landform entirely.

For terrain, use a bare-earth DTM. For an urban visualisation, a DSM is the right choice — just make it deliberately. See DEM, DSM and DTM.

2. Get the cell size in ground units

The same requirement as slope and aspect, and the same failure if you skip it:

if src.crs.is_geographic:
    lat = (src.bounds.bottom + src.bounds.top) / 2
    cell_x = abs(src.transform.a) * 111_320 * math.cos(math.radians(lat))
    cell_y = abs(src.transform.e) * 110_574

With degrees, every cell's slope saturates at 90°, and the hillshade becomes a hard black-and-white pattern with no midtones — high contrast that looks almost intentional. Check hs.std(): a real hillshade sits around 40–60, and the saturated version is much higher.

3. Understand the azimuth conversion 

az = math.radians(360.0 - azimuth + 90.0)

That line looks like magic and is doing two conversions. Compass azimuth runs clockwise from north; the arctan2 output runs anticlockwise from east. 360 - azimuth flips the direction of rotation and + 90 moves the origin.

Getting it wrong does not raise. It rotates the lighting by some fixed amount, which usually reads as "the shadows are on the wrong side" — the inverted relief effect.

4. Choose an altitude that gives contrast 

for altitude in (10, 30, 45, 60, 80):
    hs = hillshade(dem, 27.9, 30.7, altitude=altitude)
    print(f"altitude {altitude:2}°: mean {hs.mean():6.1f}  sd {hs.std():5.1f}  "
          f"range {hs.min():.0f}-{hs.max():.0f}")
altitude 10°: mean   49.2  sd  49.5  range 0-231
altitude 30°: mean  111.3  sd  61.2  range 0-254
altitude 45°: mean  159.0  sd  55.3  range 0-255
altitude 60°: mean  197.6  sd  43.1  range 13-255
altitude 80°: mean  227.9  sd  26.0  range 98-255

The pattern is clear. At 80° the standard deviation collapses to 26 and the range never drops below 98 — the sun is nearly overhead, everything is lit, and the image is washed out. At 10° the mean falls to 49 and most of the map is in shadow.

Contrast actually peaks at 30° here (sd 61.2) rather than at 45° (55.3). 45° remains the convention because it keeps more of the map out of deep shadow while staying close to the peak; 30–35° is the right choice when relief is gentle and needs exaggerating.

5. Combine with elevation colour

A hillshade alone is grey. The standard cartographic product is a colour ramp with the hillshade multiplied over it:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(9, 9))
ax.imshow(dem, cmap="terrain", extent=extent)
ax.imshow(hs, cmap="gray", alpha=0.45, extent=extent)
ax.set_axis_off()

alpha around 0.4–0.5 is the usable range. Higher and the colour is lost; lower and the relief disappears.

Hillshade standard deviation peaking around 30 to 45 degrees altitude and collapsing at 80 degrees.
Contrast peaks in the 30–45° band. At 80° the image is washed out; at 10° it is mostly shadow.

Code examples

Example 1 — a complete hillshade with the inputs validated 

import math

import numpy as np
import rasterio


def cell_size_metres(src):
    if not src.crs.is_geographic:
        return abs(src.transform.a), abs(src.transform.e)
    lat = (src.bounds.bottom + src.bounds.top) / 2
    return (abs(src.transform.a) * 111_320 * math.cos(math.radians(lat)),
            abs(src.transform.e) * 110_574)


def hillshade_from_dem(path, *, azimuth=315.0, altitude=45.0,
                       z_factor=1.0, nodata_below=-100):
    """Standard Horn-based hillshade, 0-255, with NoData handled."""
    with rasterio.open(path) as src:
        dem = src.read(1).astype("float64")
        if src.nodata is not None:
            dem = np.where(dem == src.nodata, np.nan, dem)
        cell_x, cell_y = cell_size_metres(src)
        profile = src.profile.copy()

    invalid = ~np.isfinite(dem) | (dem < nodata_below)
    working = np.where(invalid, np.nanmean(dem[~invalid]), dem) * z_factor

    p = np.pad(working, 1, mode="edge")
    a, b, c = p[:-2, :-2], p[:-2, 1:-1], p[:-2, 2:]
    d,    f = p[1:-1, :-2],               p[1:-1, 2:]
    g, h, i = p[2:, :-2],  p[2:, 1:-1],  p[2:, 2:]
    dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cell_x)
    dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cell_y)

    slope = np.arctan(np.hypot(dz_dx, dz_dy))
    aspect = np.arctan2(dz_dy, -dz_dx)

    az = math.radians(360.0 - azimuth + 90.0)
    alt = math.radians(altitude)
    shaded = (math.sin(alt) * np.cos(slope)
              + math.cos(alt) * np.sin(slope) * np.cos(az - aspect))
    hs = np.clip(shaded, 0, 1) * 255

    from scipy.ndimage import binary_dilation
    halo = binary_dilation(invalid, np.ones((3, 3), bool))
    hs = np.where(halo, np.nan, hs)

    finite = hs[np.isfinite(hs)]
    print(f"  az {azimuth:.0f}° alt {altitude:.0f}° z {z_factor} · "
          f"cell {cell_x:.1f}x{cell_y:.1f} m")
    print(f"  mean {finite.mean():.1f}  sd {finite.std():.1f}  "
          f"range {finite.min():.0f}-{finite.max():.0f}")
    if finite.std() < 30:
        print("  WARNING: low contrast — lower the altitude or raise z_factor")
    if finite.std() > 75:
        print("  WARNING: very high contrast — check the cell size is in metres")
    return hs, profile


hs, profile = hillshade_from_dem("snowdonia_glo30.tif")
  az 315° alt 45° z 1.0 · cell 27.9x30.7 m
  mean 159.0  sd 55.3  range 0-255

The two contrast warnings between them catch both common failures. A standard deviation under 30 means a washed-out image; over 75 means the slope has saturated because the cell size is in degrees.

Example 2 — multidirectional hillshade for complex terrain

A single light source leaves slopes facing directly away from it entirely black, losing all detail there. Blending several azimuths recovers it:

def multidirectional(path, *, azimuths=(225, 270, 315, 360), altitude=45,
                     weights=(0.15, 0.25, 0.35, 0.25)):
    """Weighted blend of several light directions — Swiss-style relief."""
    if len(azimuths) != len(weights):
        raise ValueError("one weight per azimuth")
    if not math.isclose(sum(weights), 1.0):
        raise ValueError(f"weights sum to {sum(weights)}, not 1")

    layers = []
    for az in azimuths:
        hs, profile = hillshade_from_dem(path, azimuth=az, altitude=altitude)
        layers.append(hs)

    blended = np.nansum([w * layer for w, layer in zip(weights, layers)], axis=0)
    blended = np.where(np.isnan(layers[0]), np.nan, blended)

    single = layers[2]                       # the 315° layer, for comparison
    for name, arr in [("single 315°", single), ("multidirectional", blended)]:
        f = arr[np.isfinite(arr)]
        print(f"  {name:18} sd {f.std():5.1f}  cells under 20: {(f < 20).mean():.2%}")
    return blended, profile


blended, profile = multidirectional("snowdonia_glo30.tif")
  single 315°        sd  55.3  cells under 20: 1.36%
  multidirectional   sd  41.5  cells under 20: 0.03%

The trade is explicit in those numbers. Overall contrast drops from 55.3 to 41.5, and the proportion of near-black cells — where all detail is lost — falls from 1.36% to 0.03%, a factor of forty. In rugged terrain that is a good bargain; in gentle terrain the single light is punchier.

Example 3 — writing it as a usable raster 

def write_hillshade(hs, profile, path, *, azimuth, altitude, z_factor=1.0):
    out = np.where(np.isfinite(hs), hs, 0).astype("uint8")
    profile.update(dtype="uint8", count=1, nodata=0,
                   compress="deflate", tiled=True, blockxsize=256, blockysize=256)

    with rasterio.open(path, "w", **profile) as dst:
        dst.write(out, 1)
        dst.update_tags(product="hillshade", azimuth=str(azimuth),
                        altitude=str(altitude), z_factor=str(z_factor),
                        method="horn", range="0-255")
    print(f"  {path}: uint8, az {azimuth}, alt {altitude}")
    return path


write_hillshade(hs, profile, "snowdonia_hillshade.tif", azimuth=315, altitude=45)
  snowdonia_hillshade.tif: uint8, az 315, alt 45

uint8 is the right dtype — the values are already 0–255 and it quarters the file size against float32. Using 0 as NoData costs one shade of black, which is imperceptible.

The tags matter because two hillshades of the same area with different azimuths are visually incompatible: shadows fall on opposite sides, and overlaying or comparing them is meaningless without knowing the parameters.

Explanation

Why 315° is the convention

Human vision resolves shape-from-shading by assuming light comes from above. In a two-dimensional image "above" is interpreted as the top of the frame, and slightly to the left, which corresponds to a north-west azimuth on a north-up map.

Light from the south-east — azimuth 135° — produces a technically correct image in which hills appear as hollows and valleys as ridges. The effect is compelling and hard to un-see once noticed. Nearly every viewer experiences it.

So 315° is not a default to be tuned; it is a perceptual requirement. If a specific ridge is poorly lit, adjust the altitude or use a multidirectional blend rather than moving the azimuth south.

Why altitude trades contrast against shadow

The formula's first term, sin(altitude) * cos(slope), is the base illumination that every cell receives. At high sun angles it dominates: sin(80°) is 0.98, so even steep slopes are brightly lit and the differences between them shrink.

The second term, cos(altitude) * sin(slope) * cos(az - aspect), carries the directional information. It scales with cos(altitude), so it is strongest when the sun is low.

Hence the measured pattern: standard deviation 26 at 80° and 61 at 30°. Lower sun means more relief, until the shadows become so deep that whole slopes go to zero and detail is lost there instead.

The same ridge lit from the north-west reading as a ridge, and lit from the south-east reading as a valley.
Identical terrain, opposite perception. This is why the azimuth is a convention rather than a preference.

Why z_factor exists

z_factor multiplies the elevations before the gradient is computed, exaggerating relief. It has two legitimate uses.

The first is gentle terrain. A floodplain with 5 m of relief over 10 km produces slopes under 0.1° and a hillshade that is uniformly grey. A z_factor of 5 or 10 makes the landform visible.

The second is unit mismatch. If the elevations are in feet and the cell size in metres, z_factor=0.3048 corrects it. This is the classic use in ArcGIS documentation and it is a fix for a bug rather than a styling choice.

Note that exaggerated relief is a visualisation, not a measurement — never compute slope from a z-exaggerated DEM.

Why the hillshade inherits every DEM problem

A hillshade is slope and aspect combined into an image, so everything that corrupts those corrupts it:

  • Degrees as cell size — slopes saturate, and the image becomes hard black and white with no midtones.
  • NoData in the window — a one-cell ring of hard-lit and hard-shadowed pixels around every void.
  • A DSM — every hedgerow and building rendered in sharp relief, obscuring the terrain.
  • Striping or steps in the source — amplified dramatically, because the derivative of a step is a spike.

That last one makes a hillshade an excellent quality-control tool. Artefacts invisible in an elevation ramp — tile seams, interpolation stripes, void-fill patches — are unmistakable in a hillshade. Render one of every new DEM before using it for anything.

Edge cases or notes

  • Azimuth is compass degrees, clockwise from north. The 360 - az + 90 conversion is required, and getting it wrong rotates the lighting silently.
  • A standard deviation over about 75 means saturated slopes — check the cell size is in metres.
  • clip(0, 1) before scaling is what stops negative values (cells facing directly away from the sun) wrapping to bright.
  • Do not compute slope from a z-exaggerated DEM. Exaggeration is for looking at, not measuring.
  • uint8 with 0 as NoData is the standard output. It costs one shade of black.
  • Multidirectional blends lose global contrast to recover detail on shadowed faces. Worth it in rugged terrain, not in gentle terrain.
  • Overlay alpha around 0.4–0.5 balances colour against relief. Above 0.6 the colour ramp disappears.
  • A hillshade is the fastest DEM quality check there is. Tile seams and interpolation artefacts jump out.

FAQ

What azimuth and altitude should I use?

315° and 45°. The azimuth is a perceptual convention, not a preference — lighting from the south-east makes hills look like valleys.

Why does my hillshade look washed out?

The altitude is too high. At 80° the standard deviation collapses to about 26 and nothing is in shadow. Drop to 30–45°.

Why is my hillshade harsh black and white with no midtones?

The slopes have saturated, almost always because the cell size is in degrees. Convert it to metres first.

What is z_factor for?

Exaggerating gentle relief so it is visible, or correcting a unit mismatch between vertical feet and horizontal metres. Never use a z-exaggerated DEM for measurement.

Should I use a DSM or a DTM?

A DTM for terrain. A DSM renders every tree and building in relief, which is striking for a city and obscuring for a landscape.

What is a multidirectional hillshade?

A weighted blend of several azimuths. It recovers detail on faces the single light leaves black, at the cost of overall contrast — 41.5 against 55.3 in the example here.

How do I combine it with elevation colour?

Draw the colour ramp, then the hillshade over it at alpha around 0.45. Higher loses the colour, lower loses the relief.