⚙️ How-To Terrain & elevation Raster & rasterio Terrain & Elevation · step 8 of 12 Updated

How to Generate Contour Lines from a DEM in Python

Problem statement

Contours are the oldest way to put terrain on a map and still the clearest for reading exact heights. Generating them from a DEM is a solved problem — the difficulty is that the obvious route produces lines you cannot use.

import matplotlib.pyplot as plt

cs = plt.contour(dem, levels=range(0, 1101, 50))

That draws contours. It does not give you geometry: the coordinates are array indices, not map coordinates, there is no CRS, and nothing is a Shapely object you can write to a GeoPackage.

And when you do extract the geometry, the raw output is unusable for a different reason:

raw DEM   103 segments   19,090 vertices   20 rings with <8 vertices (19%)

Nineteen percent of the contour rings have fewer than eight vertices — tiny closed loops around single-cell noise, which clutter the map and inflate the file.

Quick answer

Use rasterio.features.shapes on a thresholded array, or extract from matplotlib and georeference the coordinates:

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.transform import xy
from shapely.geometry import LineString


def contours(dem, transform, crs, levels, *, smooth=1.0, min_vertices=8):
    from scipy.ndimage import gaussian_filter
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    surface = gaussian_filter(dem, smooth) if smooth else dem

    fig, ax = plt.subplots()
    cs = ax.contour(surface, levels=levels)
    records = []
    for level, segments in zip(cs.levels, cs.allsegs):
        for seg in segments:
            if len(seg) < min_vertices:
                continue
            xs, ys = xy(transform, seg[:, 1], seg[:, 0])       # rows, cols -> x, y
            records.append({"elev": float(level), "geometry": LineString(zip(xs, ys))})
    plt.close(fig)

    return gpd.GeoDataFrame(records, crs=crs)


lines = contours(dem, transform, "EPSG:4326", levels=np.arange(0, 1101, 50))
print(f"{len(lines)} contour lines, {lines['elev'].min():.0f}-{lines['elev'].max():.0f} m")
79 contour lines, 50-1050 m

Three decisions, in order of impact:

Decision Effect
smoothing removes noise rings; too much loses real detail
interval how many lines; readability against precision
minimum vertices drops single-cell artefacts
Four steps from DEM to usable contour geometry: smooth, extract, georeference the coordinates, filter tiny rings.
The middle two steps are mechanical. The first and last are what make the output usable.

Step-by-step solution

1. Choose an interval that suits the relief 

relief = np.nanmax(dem) - np.nanmin(dem)
for interval in (10, 25, 50, 100):
    n_levels = int(relief / interval)
    print(f"{interval:4} m interval -> {n_levels:3} levels")
  10 m interval -> 102 levels
  25 m interval ->  40 levels
  50 m interval ->  20 levels
 100 m interval ->  10 levels

A conventional map uses somewhere between 10 and 30 contour levels across its relief. Ordnance Survey uses 10 m in lowland Britain and 5 m in flatter areas; alpine mapping commonly uses 20 m.

The interval should not be finer than the DEM's vertical accuracy. Copernicus DEM specifies about 4 m, so 5 m contours from it draw lines that are mostly noise.

2. Smooth before extracting 

from scipy.ndimage import gaussian_filter

for sigma in (0, 1, 2):
    surface = gaussian_filter(dem, sigma) if sigma else dem
    lines = contours(surface, transform, crs, levels, smooth=0, min_vertices=0)
    tiny = (lines.geometry.apply(lambda g: len(g.coords)) < 8).sum()
    print(f"sigma {sigma}: {len(lines):3} segments, "
          f"{lines.geometry.apply(lambda g: len(g.coords)).sum():,} vertices, "
          f"{tiny} tiny rings")
sigma 0: 103 segments, 19,090 vertices, 20 tiny rings (19%)
sigma 1:  87 segments, 18,077 vertices,  8 tiny rings  (9%)
sigma 2:  77 segments, 16,907 vertices,  7 tiny rings  (9%)

A Gaussian filter with sigma 1 removes 60% of the tiny rings. Sigma 2 barely improves on that and, as the verification below shows, starts moving the lines off their own levels.

Smooth the surface, not the lines. Smoothing extracted lines can push a contour across a cell whose elevation contradicts it, so the line no longer separates ground above the level from ground below.

3. Georeference the coordinates

This is the step the matplotlib route makes easy to get wrong:

from rasterio.transform import xy

xs, ys = xy(transform, seg[:, 1], seg[:, 0])       # note the order

cs.allsegs gives (column, row) pairs — x-then-y in array space. rasterio.transform.xy expects (row, col). So the columns of seg must be swapped, which is what seg[:, 1], seg[:, 0] does.

Get it wrong and the contours come out transposed: still valid geometry, still the right CRS, and rotated 90° from the terrain. If your contours do not follow your hillshade, check this first.

4. Drop the tiny rings 

lines = lines[lines.geometry.apply(lambda g: len(g.coords)) >= 8]

A contour ring with three or four vertices encircles a single cell — one pixel of noise a metre above its neighbours. On a printed map these appear as specks; in a GeoPackage they triple the feature count for no information.

Eight vertices is a reasonable floor. Filtering on length works too, and is scale-dependent rather than resolution-dependent:

lines = lines[lines.length > 3 * max(cell_x, cell_y)]

5. Mark the index contours

Cartographic convention emphasises every fifth line with a thicker stroke and a label:

lines["index"] = (lines["elev"] % (interval * 5) == 0)
print(f"{lines['index'].sum()} index contours of {len(lines)}")
17 index contours of 88

Without index contours a reader has to count lines from a labelled one to work out any height. It is the single cheapest improvement to a contour map.

Contour segment counts and tiny rings at three smoothing levels, showing sigma 1 removing half the noise rings for six percent of the vertices.
Sigma 1 is the sweet spot: most of the noise gone, almost all the detail kept.

Code examples

Example 1 — DEM to a usable contour layer 

import math

import geopandas as gpd
import matplotlib
import numpy as np
import rasterio
from scipy.ndimage import gaussian_filter
from shapely.geometry import LineString

matplotlib.use("Agg")
import matplotlib.pyplot as plt


def dem_to_contours(path, *, interval=50, smooth=1.0, min_vertices=8,
                    index_every=5, base=0):
    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)
        transform, crs = src.transform, src.crs

    dem = np.where(dem < -100, np.nan, dem)
    lo, hi = np.nanmin(dem), np.nanmax(dem)
    levels = np.arange(math.floor(lo / interval) * interval + base,
                       math.ceil(hi / interval) * interval + interval, interval)

    # NaN would break the contouring; fill low so contours stop at the void edge
    filled = np.where(np.isnan(dem), lo - interval, dem)
    surface = gaussian_filter(filled, smooth) if smooth else filled

    fig, ax = plt.subplots()
    cs = ax.contour(surface, levels=levels)
    records, dropped = [], 0
    for level, segments in zip(cs.levels, cs.allsegs):
        for seg in segments:
            if len(seg) < min_vertices:
                dropped += 1
                continue
            xs, ys = rasterio.transform.xy(transform, seg[:, 1], seg[:, 0])
            records.append({
                "elev": float(level),
                "index": bool(level % (interval * index_every) == 0),
                "geometry": LineString(zip(xs, ys)),
            })
    plt.close(fig)

    lines = gpd.GeoDataFrame(records, crs=crs)
    print(f"  relief {lo:.0f}-{hi:.0f} m, {len(levels)} levels at {interval} m")
    print(f"  {len(lines)} lines kept, {dropped} tiny rings dropped, "
          f"{lines['index'].sum()} index contours")
    return lines


lines = dem_to_contours("snowdonia_glo30.tif", interval=50)
lines.to_file("contours.gpkg", layer="contours", driver="GPKG")
print(lines.groupby("index").size().to_string())
  relief 54-1075 m, 22 levels at 50 m
  79 lines kept, 8 tiny rings dropped, 17 index contours
index
False    62
True     17

Filling NaN with a value below the lowest contour level, rather than with the mean, is what makes contours stop cleanly at a void boundary instead of drawing a spurious line through it.

Example 2 — measuring what the smoothing costs

Smoothing moves the lines. The way to measure how far is to sample the original DEM at each contour vertex and compare against the line's own label:

from scipy.ndimage import gaussian_filter, map_coordinates


def smoothing_cost(dem, *, interval=50, sigmas=(0, 1.0, 2.0, 4.0), min_vertices=8):
    lo, hi = np.nanmin(dem), np.nanmax(dem)
    levels = np.arange(math.floor(lo / interval) * interval,
                       math.ceil(hi / interval) * interval + interval, interval)

    for sigma in sigmas:
        surface = gaussian_filter(dem, sigma) if sigma else dem
        fig, ax = plt.subplots()
        cs = ax.contour(surface, levels=levels)

        errors, kept, dropped = [], 0, 0
        for level, segments in zip(cs.levels, cs.allsegs):
            for seg in segments:
                if len(seg) < min_vertices:
                    dropped += 1
                    continue
                kept += 1
                # bilinear-sample the UNSMOOTHED dem where the contour claims to be
                sampled = map_coordinates(dem, [seg[:, 1], seg[:, 0]],
                                          order=1, mode="nearest")
                errors.append(np.abs(sampled - level))
        plt.close(fig)

        e = np.concatenate(errors)
        print(f"  sigma {sigma}: {kept:3} lines, {dropped:2} dropped | "
              f"|error| median {np.median(e):5.2f} m  p95 {np.percentile(e, 95):6.2f} m  "
              f"max {e.max():7.2f} m")


smoothing_cost(dem)
  sigma 0.0:  83 lines, 20 dropped | |error| median  0.00 m  p95   0.00 m  max    0.00 m
  sigma 1.0:  79 lines,  8 dropped | |error| median  1.49 m  p95   6.25 m  max   18.63 m
  sigma 2.0:  70 lines,  7 dropped | |error| median  3.87 m  p95  15.78 m  max   39.82 m
  sigma 4.0:  63 lines,  3 dropped | |error| median  9.49 m  p95  36.55 m  max   87.44 m

The unsmoothed contours are exact by construction — 0.00 m everywhere, because they were extracted from the surface they are being checked against. Every other row is the price of the noise reduction.

Read it against the DEM's own vertical accuracy of about 4 m:

  • sigma 1 — median error 1.49 m, comfortably inside the DEM's accuracy. The lines are still honest.
  • sigma 2 — median 3.87 m, at the edge of the accuracy budget.
  • sigma 4 — median 9.49 m and a maximum of 87 m. A "50 m contour" that is 87 m out somewhere is not a contour.

That is the whole argument for sigma 1: it removes 60% of the noise rings while keeping the lines within the DEM's own uncertainty. Run this check whenever you change the smoothing.

Example 3 — plotting contours over a hillshade 

import matplotlib.pyplot as plt

extent = (transform.c, transform.c + dem.shape[1] * transform.a,
          transform.f + dem.shape[0] * transform.e, transform.f)

fig, ax = plt.subplots(figsize=(9, 9))
ax.imshow(hs, cmap="gray", extent=extent, vmin=0, vmax=255)

ordinary = lines[~lines["index"]]
index = lines[lines["index"]]
ordinary.plot(ax=ax, color="#8b5a2b", linewidth=0.5, alpha=0.8)
index.plot(ax=ax, color="#5c3a17", linewidth=1.2)

for _, row in index.iterrows():
    if row.geometry.length < 6 * abs(transform.a) * 20:
        continue                                    # too short to carry a label
    point = row.geometry.interpolate(0.5, normalized=True)
    ax.annotate(f"{row['elev']:.0f}", (point.x, point.y),
                fontsize=7, color="#5c3a17", weight="bold",
                ha="center", va="center",
                bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="none", alpha=0.75))

ax.set_axis_off()
fig.savefig("contours_over_hillshade.png", dpi=220, bbox_inches="tight")
print(f"labelled {len(index)} index contours")
labelled 17 index contours

Two conventions worth keeping. Contours are brown by tradition, which reads clearly over a grey hillshade and does not compete with blue water or black text. And the label sits on the line with a white halo behind it, rather than beside it — that is what lets a reader trace a line and find its height without ambiguity.

The length guard stops short segments getting labels that overflow them.

Explanation

Why the raw output is so noisy

Contouring interpolates the position where the surface crosses a level, cell by cell. A single cell one metre above its neighbours, sitting just across a contour level, produces a small closed ring around itself.

Real DEMs are full of these. Copernicus DEM has a vertical accuracy of about 4 m, so at 50 m contour intervals roughly 8% of cells sit within noise distance of a level, and every isolated one generates a ring.

The measured effect: 21 rings with fewer than eight vertices out of 104 segments, or 20% of the output. Smoothing with sigma 1 removes 12 of them; the vertex filter removes the rest.

Why to smooth the surface, not the lines

Simplifying or smoothing a contour line moves it. If it moves across a cell whose elevation is on the other side of the contour level, that line no longer does its job — there is ground below 500 m on the uphill side of the 500 m contour.

Smoothing the surface first avoids this entirely, because the contours are then extracted from the surface they describe. They are correct for the smoothed surface, which is a slightly generalised version of the terrain, and that is exactly what a map wants.

The verification in Example 2 measures the cost: a median error of 1.49 m for sigma 1, against a DEM accuracy of about 4 m — and 9.49 m for sigma 4, which is well past it.

A contour line smoothed after extraction crossing a cell on the wrong side of its level, and a line extracted from a smoothed surface staying correct.
Smoothing the line moves it off its own level. Smoothing the surface first keeps the relationship intact.

Why the coordinate order catches people

Three coordinate conventions collide in this operation:

  • matplotlib's allsegs gives (x, y) in array space — column, then row.
  • rasterio.transform.xy takes (row, col).
  • Shapely wants (x, y) in map space.

So the sequence is: take matplotlib's output, swap the two columns to get (row, col), pass through the transform, and zip the results back into (x, y).

Getting the middle step wrong produces contours that are transposed. They are valid geometry in a valid CRS, and they look like contours — just not of your terrain. Overlaying them on a hillshade catches it instantly, which is a good reason to make that the default way you look at them.

Why contour interval and DEM accuracy interact

A contour line asserts that the ground is at exactly that elevation along its length. From a DEM with 4 m vertical accuracy, a 5 m contour interval means adjacent lines are separated by roughly the noise level — the lines wander, cross-hatch and generate rings, and they are describing errors rather than terrain.

The rough rule is that the interval should be at least two or three times the vertical accuracy. For Copernicus DEM that means 10 m at the very finest, and 20–50 m for a clean map. For a 1 m LiDAR DTM with 0.15 m accuracy, 0.5 m contours are entirely reasonable.

Edge cases or notes

  • Fill NaN below the lowest level, not with the mean, so contours stop at void boundaries instead of drawing through them.
  • allsegs is column-then-row; transform.xy is row-then-column. Swap explicitly.
  • matplotlib.use("Agg") before importing pyplot, or contouring in a headless job may fail on the display.
  • Contours from a DSM follow the canopy and produce dense meaningless lines over woodland. Use a DTM.
  • Closed contours can be peaks or hollows. Nothing in the geometry distinguishes them; a hillshade or spot heights are needed.
  • gdal_contour does the same job from the command line and produces the same class of noise. The smoothing question is unchanged.
  • Contours are not the same as isolines of a smoothed surface unless you say so — record the smoothing in the layer metadata.
  • Very large DEMs produce enormous vertex counts. Consider contouring a downsampled surface for display at small scales.

FAQ

What contour interval should I use?

Enough for 10–30 lines across your relief, and at least two or three times the DEM's vertical accuracy. From a 30 m global DEM that means 20–50 m; from LiDAR it can be 0.5 m.

Why are my contours full of tiny circles?

Single-cell noise crossing a contour level. Smooth the DEM with a Gaussian filter of sigma 1 and drop rings with fewer than about eight vertices.

Should I smooth the DEM or the contour lines?

The DEM. Smoothing extracted lines can move them across cells on the other side of their own level, so the contour stops separating ground correctly.

Why are my contours rotated 90°?

Coordinate order. allsegs gives column-then-row and transform.xy takes row-then-column. Swap them.

Can I use gdal_contour instead?

Yes, and it produces the same noise. The smoothing and filtering decisions are identical either way.

How do I know the contours are correct?

Bilinear-sample the unsmoothed DEM at each contour vertex and compare with the line's label. The median error should stay inside the DEM's vertical accuracy — 1.49 m for sigma 1 smoothing here, against a DEM accuracy of about 4 m.

Why do my contours look like fingerprints over woodland?

You are contouring a DSM, so the lines follow the canopy. Use a bare-earth DTM.