Generating Hillshade Rasters in Python

Combine slope and aspect with a light direction, then scale for display:

import numpy as np

zenith = np.radians(90 - altitude_deg)                    # 45° altitude by convention
azimuth = np.radians(azimuth_deg)                          # 315° = north-west light
shade = (np.cos(zenith) * np.cos(slope)
         + np.sin(zenith) * np.sin(slope) * np.cos(azimuth - aspect))
hillshade = np.clip(shade, 0, 1) * 255

Slope and aspect in radians, aspect clockwise from north — and the rest is one formula. This page belongs to terrain analysis and DEM-derived products in Satellite Processing Workflows & Index Pipelines.

The Illumination Formula

Brightness is the cosine of the incidence angle A slope's surface normal points away from the ground perpendicular to it. The light arrives from a direction set by azimuth and altitude. Hillshade brightness is the cosine of the angle between the normal and the light: one when the slope faces the light directly, zero when the light grazes the surface or comes from behind. cos(incidence) = brightness surface normal light, 315° at 45° brightness 1.0 facing the light 0.7 flat ground at 45° 0.0 facing away Flat ground is not white: it receives light at the sun's altitude, not head-on.

The formula is the cosine of the angle between the slope’s surface normal and the light direction, expanded in terms of slope, aspect, and the light’s zenith and azimuth. Flat ground gets the cosine of the zenith — about 0.7 for a 45° light — which is why hillshade images have a mid-grey base with bright and dark slopes on either side, rather than white plains.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Trigonometry
scipy >=1.11 Gradient kernels
rasterio >=1.3.0 Reading the DEM and writing uint8 output
matplotlib >=3.8 Optional: blending hillshade under a colour layer 
pip install "numpy>=1.23" "scipy>=1.11" "rasterio>=1.3.0" "matplotlib>=3.8"

Complete Working Example

import numpy as np
import rasterio
from scipy import ndimage

KX = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]) / 8.0
KY = np.array([[1, 2, 1], [0, 0, 0], [-1, -2, -1]]) / 8.0


def hillshade(dem: np.ndarray, res: float, *, azimuth: float = 315.0,
              altitude: float = 45.0, z_factor: float = 1.0) -> np.ndarray:
    dzdx = ndimage.convolve(dem, KX, mode="nearest") / res * z_factor
    dzdy = ndimage.convolve(dem, KY, mode="nearest") / res * z_factor
    slope = np.arctan(np.hypot(dzdx, dzdy))
    aspect = np.arctan2(-dzdx, -dzdy)                        # downslope, clockwise from north

    zen = np.radians(90.0 - altitude)
    az = np.radians(azimuth)
    shade = (np.cos(zen) * np.cos(slope)
             + np.sin(zen) * np.sin(slope) * np.cos(az - aspect))
    return np.clip(shade, 0, 1).astype("float32")


def write_hillshade(dem_path: str, out_path: str, **kwargs) -> None:
    with rasterio.open(dem_path) as src:
        dem = src.read(1, masked=True).filled(np.nan).astype("float64")
        res = abs(src.transform.a)
        profile = src.profile | {"dtype": "uint8", "nodata": 0, "count": 1,
                                 "compress": "deflate", "tiled": True}
    hs = hillshade(np.nan_to_num(dem, nan=np.nanmean(dem)), res, **kwargs)
    out = (hs * 254 + 1).astype("uint8")
    out[~np.isfinite(dem)] = 0
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(out, 1)
        dst.update_tags(azimuth=str(kwargs.get("azimuth", 315)),
                        altitude=str(kwargs.get("altitude", 45)))

Storing hillshade as uint8 with 0 reserved for nodata keeps it tiny and makes it directly usable as a display layer. It is a picture of relief rather than a measurement, so eight bits lose nothing that matters — the same reasoning as in converting rasters to 8-bit for display.

Multidirectional Shading

One light hides some ridges With a single light from the north-west, ridges running north-west to south-east are lit evenly on both sides and nearly disappear. Averaging hillshades from several azimuths — for example 225, 270, 315 and 360 degrees — lights every orientation from at least one direction, so all ridges remain visible. single light at 315° four azimuths averaged NW–SE ridge almost invisible both ridges defined Multidirectional shading costs four times the compute and removes orientation bias. 
def multidirectional(dem, res, azimuths=(225, 270, 315, 360), altitude=45.0):
    return np.mean([hillshade(dem, res, azimuth=a, altitude=altitude) for a in azimuths], axis=0)

A single light direction hides features aligned with it: a ridge running parallel to the light is lit equally on both flanks and nearly vanishes. Averaging several azimuths removes that bias, at the cost of a slightly flatter overall look. For base maps covering varied terrain, multidirectional shading is usually the better default; for a small area with one dominant grain, a single well-chosen azimuth can be more striking.

Hillshade as a Base Layer

Hillshade is rarely shown alone. Its main use is underneath other data — a land-cover map, an NDVI layer, an elevation colour ramp — where it adds a sense of relief without competing for attention. The standard technique is to blend it multiplicatively with the colour layer, so bright slopes leave the colour unchanged and shaded slopes darken it.

import matplotlib.pyplot as plt
from matplotlib.colors import LightSource

ls = LightSource(azdeg=315, altdeg=45)
rgb = ls.shade(elevation, cmap=plt.cm.gist_earth, blend_mode="soft",
               vert_exag=1.5, dx=res, dy=res)

Keep the hillshade subtle in such blends. At full strength it overwhelms the colour and makes classes on shaded slopes hard to read; a soft blend, or the hillshade drawn at around 40% opacity beneath a translucent data layer, gives relief without distortion. The layer-ordering advice in interactive raster exploration in Jupyter applies directly.

Verification

import numpy as np

hs = hillshade(dem, res=30.0)
flat = np.cos(np.radians(90 - 45))
assert 0.0 <= hs.min() and hs.max() <= 1.0
assert abs(np.median(hs) - flat) < 0.15, "median far from the flat-ground value — check the DEM units"
A healthy hillshade histogram Over mixed terrain the hillshade histogram peaks near the flat-ground value, the cosine of the light's zenith, with tails toward bright sun-facing and dark shaded slopes. A histogram piled at zero and one instead means every slope is being treated as vertical, the signature of a DEM still in degrees. Expected distribution of shade values flat ground ≈ 0.71 0 shaded 1 facing light Spikes at 0 and 1 mean the DEM is in degrees, not metres.

The median of a hillshade over mixed terrain sits close to the flat-ground value, the cosine of the light’s zenith. A median far from it — every pixel near white or near black — usually means the DEM is still in degrees, so every slope is effectively vertical. Rendering the result and checking that valleys look like valleys is the final, most reliable test.

Common Errors

Relief looks inverted

The light comes from the south-east, or the aspect convention is flipped. Use azimuth 315 and check aspect against a known slope before changing anything else.

The whole image is near white or near black

The DEM is in geographic degrees. Reproject to metres before shading.

Hillshade is noisy and speckled

The DEM carries vertical noise that the gradient amplifies. Smooth lightly first, or use a coarser DEM.

Flat regions show no relief at all

Relief is genuinely subtle. Apply a vertical exaggeration of two to five, and record the factor used in the output tags.

Frequently Asked Questions

Q: Why is hillshade lit from the north-west? Convention and perception. Light from the upper left of a north-up map makes relief read correctly for most viewers; light from the south-east makes valleys look like ridges, an illusion known as relief inversion.

Q: Should hillshade use the real sun position? For display, no — a fixed conventional light direction reads best. For correcting imagery, yes — the illumination must match the sun at the moment of acquisition, which is a different product with a different purpose.

Q: What is vertical exaggeration for? Making subtle relief visible. Multiplying the gradients by a factor of two to five brings out low hills in flat regions; in mountains it saturates and should be left at one.

Q: Can hillshade be served as web tiles? Yes, and it is an ideal candidate: small, static, and reused under many layers. Build it once, convert it to a COG with overviews, and serve it like any other layer, as covered in serving raster tiles with TiTiler.

Q: Does hillshade need to be recomputed for each map? No. For a given DEM and light direction it never changes, so compute it once, store it as a tiled, compressed layer with overviews, and reuse it under every map of that area. Recomputing it per figure is pure waste, and storing it guarantees every map shares the same relief shading.