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
matplotlib.use("Agg")
import matplotlib.pyplot as plt

cs = plt.contour(dem, levels=range(0, 1101, 50))
print(type(cs).__name__, cs.allsegs[10][0][:2])

QuadContourSet [[12.59194567  0.        ]
 [13.          0.58171698]]

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 (step 2 has the code):

raw DEM   104 segments   19,227 vertices   21 rings with <8 vertices (20%)

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

78 contour lines, 100-1050 m

Three decisions, in order of impact:

An interval for the relief, smoothing, georeferencing

1. Choose an interval that suits the relief

levels = np.arange(0, 1101, 50)
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's 1:25,000 Explorer maps use 5 m contours in lowland Britain and 10 m in the mountains; its 1:50,000 Landranger maps use 10 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: 104 segments, 19,227 vertices, 21 tiny rings
sigma 1:  87 segments, 18,153 vertices, 9 tiny rings
sigma 2:  77 segments, 16,995 vertices, 7 tiny rings

A Gaussian filter with sigma 1 removes more than half of the tiny rings, 12 of 21. 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 = contours(dem, transform, crs, levels, smooth=1.0, min_vertices=0)
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:

cell = max(abs(transform.a), abs(transform.e))      # in the layer's own units — degrees here
print(f"{(lines.length > 3 * cell).sum()} of {len(lines)} lines longer than three cells")

78 of 78 lines longer than three cells

5. Mark the index contours

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

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

17 index contours of 78

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.

A contour layer over a hillshade

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.

import numpy as np
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"}
    crs = src.crs
with rasterio.open("snowdonia_glo30.tif", "w", **profile) as dst:
    dst.write(mosaic)
dem = mosaic[0].astype("float64")
print(f"snowdonia_glo30.tif {dem.shape[1]} × {dem.shape[0]}, {dem.min():.0f} to {dem.max():.0f} m, {crs}")

snowdonia_glo30.tif 240 × 252, 54 to 1075 m, EPSG:4326

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)

    # the Gaussian filter needs numbers in the gaps; contouring must not see them
    gaps = np.isnan(dem)
    filled = np.where(gaps, np.nanmean(dem), dem)
    smoothed = gaussian_filter(filled, smooth) if smooth else filled
    surface = np.ma.masked_array(smoothed, mask=gaps)

    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
  78 lines kept, 8 tiny rings dropped, 17 index contours
index
False    61
True     17

The gaps are filled only so the Gaussian filter has numbers to work with, then masked again before contouring. Contouring a filled gap — with the mean, or with a value below every level — draws contours around its edge; a masked array makes them stop there. Contours come out jagged measures the difference.

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:  83 lines, 20 dropped | |error| median  0.00 m  p95   0.00 m  max    0.00 m
  sigma 1.0:  78 lines,  8 dropped | |error| median  1.49 m  p95   6.27 m  max   18.63 m
  sigma 2.0:  70 lines,  6 dropped | |error| median  3.86 m  p95  15.76 m  max   39.82 m
  sigma 4.0:  64 lines,  3 dropped | |error| median  9.49 m  p95  36.48 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 more than half 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 math

def hillshade(z, cell_x=27.9, cell_y=30.7, azimuth=315, altitude=45):
    gy, gx = np.gradient(z, cell_y, cell_x)
    slope, aspect = np.arctan(np.hypot(gx, gy)), np.arctan2(-gy, -gx)
    az, alt = math.radians(360 - azimuth + 90), math.radians(altitude)
    return np.clip(math.sin(alt) * np.cos(slope) + math.cos(alt) * np.sin(slope) * np.cos(az - aspect), 0, 1) * 255

hs = hillshade(dem)
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.

Noisy rings and smoothing the surface, not the lines

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, and at 50 m intervals 15% of the cells in this window sit within 4 m of a contour level — close enough that noise decides which side they fall on. Where such a cell is isolated, it 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.

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.

NaN masking, allsegs order and gdal_contour

• Mask NaN rather than filling it. Any fill value — the mean, or one below every level — makes contours trace the edge of the gap.
• 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.

Internal links

• Contours come out jagged or full of tiny rings — diagnosing and fixing noisy output
• Digital elevation models explained — vertical accuracy and its effect on the interval
• How to make a hillshade from a DEM in Python — the standard companion layer
• How to simplify geometry in Python with GeoPandas and Shapely — why simplifying contours is different
• The raster data model explained — the transform used to georeference
• Raster and vector do not line up in Python — when the contours land in the wrong place
• How to read and write GeoPackage files in Python — writing the contour layer
• How to add map legends and labels in GeoPandas — labelling index contours
• How to generate depth contours from bathymetry — graduated intervals for a bimodal surface

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.