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


with rasterio.open("snowdonia_glo30.tif") as src:
    dem = src.read(1).astype("float64")
    extent = (src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top)

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.1  sd 55.3  range 0-255

Two conventions and two parameters:

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 DTM, ground units, azimuth conversion, altitude

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:

with rasterio.open("snowdonia_glo30.tif") as src:
    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
    degrees = src.res

for label, (cx, cy) in [("metres", (cell_x, cell_y)), ("degrees", degrees)]:
    h = hillshade(dem, cx, cy)
    print(f"{label:8} sd {h.std():5.1f}   black {np.mean(h == 0):6.1%}   "
          f"mid-tones 64-192 {np.mean((h > 64) & (h < 192)):6.1%}")

metres   sd  55.3   black   0.4%   mid-tones 64-192  60.0%
degrees  sd  70.5   black  51.7%   mid-tones 64-192  39.2%

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

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}  black {np.mean(hs == 0):4.0%}  white {np.mean(hs == 255):4.0%}")

altitude 10°: mean   49.3  sd  49.5  range 0-231  black  32%  white   0%
altitude 30°: mean  111.4  sd  61.2  range 0-254  black   7%  white   0%
altitude 45°: mean  159.1  sd  55.3  range 0-255  black   0%  white   0%
altitude 60°: mean  197.7  sd  43.1  range 13-255  black   0%  white   0%
altitude 80°: mean  228.0  sd  25.9  range 98-255  black   0%  white   0%

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()
fig.savefig("snowdonia_relief.png", dpi=100, bbox_inches="tight")

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

Single and multidirectional hillshades of Snowdon

The data used on this page

The Copernicus DEM (GLO-30) over the Snowdon massif, 0.1° × 0.07° around the summit, read from its public bucket. The window crosses the 4°W tile boundary, so two tiles are merged. It is a DSM, but on bare mountain there is little standing on it.

import rasterio
from rasterio.merge import merge

BBOX = (-4.10, 53.03, -4.00, 53.10)
COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
       "Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
tiles = [COP.format("W005"), COP.format("W004")]
mosaic, transform = merge(tiles, bounds=BBOX)
with rasterio.open(tiles[0]) as src:
    profile = src.profile | {"driver": "GTiff", "width": mosaic.shape[2], "height": mosaic.shape[1],
                             "transform": transform, "compress": "deflate"}
with rasterio.open("snowdonia_glo30.tif", "w", **profile) as dst:
    dst.write(mosaic)
print(f"snowdonia_glo30.tif {mosaic.shape[2]} × {mosaic.shape[1]}, {mosaic.min():.1f} to {mosaic.max():.1f} m")

snowdonia_glo30.tif 240 × 252, 54.0 to 1074.5 m

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 np.mean(finite == 0) > 0.2:
        print(f"  WARNING: {np.mean(finite == 0):.0%} of cells pure black — 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.1  sd 55.3  range 0-255

The two warnings between them catch both common failures. A standard deviation under 30 means a washed-out image. The degree mistake is not caught by the standard deviation — step 2's saturated version had an sd of 70.5, well inside what real terrain can produce — but by the black share: 51.7% of its cells were pure black, against 0.4% for the real hillshade.

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")

  az 225° alt 45° z 1.0 · cell 27.9x30.7 m
  mean 169.1  sd 53.7  range 0-255
  az 270° alt 45° z 1.0 · cell 27.9x30.7 m
  mean 163.8  sd 53.6  range 0-255
  az 315° alt 45° z 1.0 · cell 27.9x30.7 m
  mean 159.1  sd 55.3  range 0-255
  az 360° alt 45° z 1.0 · cell 27.9x30.7 m
  mean 157.8  sd 54.9  range 0-255
  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.

The 315° convention and the altitude trade

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.

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.

Compass azimuths, black maps and overlay alpha

• Azimuth is compass degrees, clockwise from north. The 360 - az + 90 conversion is required, and getting it wrong rotates the lighting silently.
• Half the map pure black means saturated slopes — check the cell size is in metres. The standard deviation alone will not tell you.
• 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.

Internal links

• Hillshade looks flat, washed out or inverted — diagnosing a hillshade that fails
• How to calculate slope and aspect from a DEM in Python — the two inputs
• Slope and aspect explained — the conventions behind them
• Digital elevation models explained — DTM versus DSM for terrain rendering
• How to plot multiple layers in GeoPandas — overlaying the hillshade under vector data
• How to add a basemap to a GeoPandas plot — the alternative backdrop
• The raster data model explained — dtype and NoData for the output
• How to generate contour lines from a DEM in Python — the classic hillshade companion
• How to compute building shadows for a given time in Python — real sun positions rather than a fixed azimuth
• How to compute a horizon profile from a DEM — shading over a whole year, not one sun position

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.