Contours Come Out Jagged or Full of Tiny Rings

Problem statement

The contour layer generated, and it is not a map:

• Hundreds of tiny closed rings, each two or three vertices, scattered across gentle ground.
• Staircase lines that follow cell edges instead of the terrain.
• A file with ten times the vertices the map needs.
• Contours over woodland that look like fingerprints, dense and meaningless.
• Lines that cross each other, which contours cannot do.

Every one of these is real information in the DEM — noise, quantisation, canopy — faithfully converted into geometry. The contouring algorithm is not wrong. It is doing exactly what you asked with a surface that was not ready.

Quick answer

Measure the problem before treating it. Here is the Copernicus DEM of the Snowdon massif, reprojected to British National Grid at 30 m and contoured at 50 m with nothing done to it:

raw_lines = extract(dem, levels_for(dem, 50), transform, crs)
vertices = raw_lines["vertices"]
print(f"{len(raw_lines)} lines, {vertices.sum():,} vertices")
print(f"tiny rings (<8 verts): {(vertices < 8).sum()} ({(vertices < 8).mean():.0%})")
print(f"median vertices per line: {vertices.median():.0f}")

117 lines, 27,015 vertices
tiny rings (<8 verts): 18 (15%)
median vertices per line: 67

Smooth to an error budget, drop rings, find staircases

1. Smooth the surface, and measure what it costs

from scipy.ndimage import gaussian_filter, map_coordinates


def contour_and_check(surface, original, interval):
    """Contour the surface; sample the ORIGINAL where each line claims to be."""
    fig, ax = plt.subplots()
    cs = ax.contour(surface, levels=levels_for(original, interval))
    kept, tiny, errors = 0, 0, []
    for level, segments in zip(cs.levels, cs.allsegs):
        for seg in segments:
            if len(seg) < 8:
                tiny += 1
                continue
            kept += 1
            sampled = map_coordinates(original, [seg[:, 1], seg[:, 0]], order=1, mode="nearest")
            errors.append(np.abs(sampled - level))
    plt.close(fig)
    return kept, tiny, np.concatenate(errors)


for sigma in (0, 1.0, 2.0, 4.0):
    surface = gaussian_filter(dem, sigma) if sigma else dem
    kept, tiny, errors = contour_and_check(surface, dem, interval=50)
    print(f"sigma {sigma}: {kept:3} lines, {tiny:2} tiny  |  "
          f"|error| median {np.median(errors):5.2f} m  max {errors.max():6.2f} m")

sigma 0:  99 lines, 20 tiny  |  |error| median  0.00 m  max   0.00 m
sigma 1.0:  94 lines,  8 tiny  |  |error| median  1.19 m  max  15.05 m
sigma 2.0:  83 lines,  5 tiny  |  |error| median  3.24 m  max  39.38 m
sigma 4.0:  73 lines,  5 tiny  |  |error| median  8.22 m  max  86.09 m

The error is measured by sampling the unsmoothed DEM where each contour claims to be. Sigma 1 removes most of the tiny rings for a median error well inside Copernicus DEM's own 4 m vertical accuracy, so the lines are still honest. Sigma 4 removes nearly all of them and moves the lines by more than the DEM's accuracy — a 50 m contour that is tens of metres out somewhere is no longer a contour.

Smoothing is not free. Measure it, and stop before the error reaches the DEM's accuracy.

2. Drop the tiny rings that survive

lines = extract(gaussian_filter(dem, 1.0), levels_for(dem, 50), transform, crs)
lines = lines[lines["vertices"] >= 8]

A three-vertex closed ring encircles one cell. It carries no information a reader can use and it inflates the feature count.

Filtering on length is the scale-aware alternative — the layer is in metres, so the threshold can be too:

min_length = 4 * 30                            # four 30 m cells
print(f"kept {(lines.length > min_length).sum()} of {len(lines)} lines longer than {min_length} m")

kept 94 of 94 lines longer than 120 m

3. Diagnose staircase lines

Staircases mean the contour is following cell boundaries rather than interpolating across them, which happens when the elevation is quantised. Copernicus ships float32; plenty of DEMs arrive as whole metres:

for name, z in [("float32 as published", dem), ("rounded to int16", np.round(dem).astype("int16"))]:
    unique = np.unique(z)
    print(f"{name:21} {len(unique):>7,} distinct values, smallest step {np.diff(unique).min():.4f} m")

float32 as published   98,205 distinct values, smallest step 0.0000 m
rounded to int16        1,010 distinct values, smallest step 1.0000 m

An int16 DEM has 1 m quantisation. Contouring it at 1 m intervals means every contour sits exactly on a quantisation step, so it traces the boundary between two integer values — a staircase.

Two fixes: cast to float and smooth slightly, or use an interval several times the quantum.

4. Recognise the DSM texture

Epping Forest, 2 km square, from the Environment Agency's lidar surface model (DSM) and terrain model (DTM) at 5 m, contoured at 5 m:

for name in ("epping_dsm.tif", "epping_dtm.tif"):
    with rasterio.open(name) as src:
        z = src.read(1, masked=True).filled(np.nan).astype("float64")
        found = extract(z, levels_for(z, 5), src.transform, crs)
    short = found[found.length < 200]
    print(f"{name}: {len(found):,} lines, {len(short) / len(found):.0%} shorter than 200 m, "
          f"{found['vertices'].sum():,} vertices")

epping_dsm.tif: 16,878 lines, 99% shorter than 200 m, 218,942 vertices
epping_dtm.tif: 118 lines, 69% shorter than 200 m, 19,485 vertices

The surface model gives over a hundred times as many lines, almost all of them short: each tree crown is a hill of its own. No smoothing parameter fixes it — a 20 m tree is a real 20 m feature of the surface, and the contour algorithm is right to draw it. Use a bare-earth DTM. See DEM, DSM and DTM.

5. Keep contours off the edge of missing data

Contouring needs a value in every cell, so gaps are usually filled with something. The full reprojected grid has empty corners where the rotated lon/lat window did not reach. Try the two common fills, and the alternative of not filling at all:

from rasterio.features import shapes
from shapely.geometry import shape
from shapely.ops import unary_union

void = unary_union([shape(g) for g, v in shapes(np.isnan(dem_full).astype("uint8"),
                                                 transform=bng_transform) if v == 1])
along_void = void.buffer(45)                       # within 1.5 cells of the gap
full_levels = levels_for(dem_full, 50)
for label, surface in [
        ("filled with the mean", np.where(np.isnan(dem_full), np.nanmean(dem_full), dem_full)),
        ("filled below every level", np.where(np.isnan(dem_full), full_levels.min() - 50, dem_full)),
        ("masked, not filled", np.ma.masked_invalid(dem_full))]:
    found = extract(surface, full_levels, bng_transform, crs)
    print(f"{label:25} {found.length.sum() / 1000:6.1f} km of contour, "
          f"{found.intersection(along_void).length.sum() / 1000:5.1f} km of it along the gap")

filled with the mean       861.2 km of contour, 153.9 km of it along the gap
filled below every level  1020.8 km of contour, 313.5 km of it along the gap
masked, not filled         713.2 km of contour,   5.9 km of it along the gap

Both fills draw contours around the hole. The mean puts a flat plateau at 404 m, so the levels either side of it trace the gap's outline. Filling below every level is worse, not better: every level between the fill and the ground at the gap's edge now crosses there, so a whole bundle of contours runs along the boundary. A masked array is the fix — matplotlib's contouring skips masked cells, and the lines stop at the edge of the data:

surface = np.ma.masked_invalid(dem)

A noisy-contour diagnostic and a clean-up pass

The data used on this page

The Copernicus DEM (GLO-30) over the Snowdon massif, read from its public bucket and reprojected to British National Grid at 30 m, so that lengths are in metres; an inner window with no gaps is used everywhere except step 5. And the Environment Agency's lidar DSM and DTM for 2 km² of Epping Forest at 5 m. extract() is the matplotlib-and-rasterio contouring from the contours how-to, returning a GeoDataFrame with a vertex count per line.

import numpy as np
import rasterio
import requests
import geopandas as gpd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from rasterio.merge import merge
from rasterio.transform import from_origin, xy
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.windows import Window
from shapely.geometry import LineString

BBOX = (-4.14, 53.02, -3.98, 53.11)                    # the Snowdon massif, WGS84
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")
mosaic, src_transform = merge([COP.format("W005"), COP.format("W004")], bounds=BBOX)
crs = rasterio.crs.CRS.from_epsg(27700)
bng_transform, width, height = calculate_default_transform(
    "EPSG:4326", crs, mosaic.shape[2], mosaic.shape[1], *BBOX, resolution=30)
dem_full = np.full((height, width), np.nan, dtype="float32")
reproject(mosaic[0], dem_full, src_transform=src_transform, src_crs="EPSG:4326",
          dst_transform=bng_transform, dst_crs=crs, dst_nodata=np.nan, resampling=Resampling.bilinear)

rows, cols = np.where(np.isfinite(dem_full).all(axis=1))[0], np.where(np.isfinite(dem_full).all(axis=0))[0]
inner = Window(20, 20, width - 40, height - 40)
dem = dem_full[20:-20, 20:-20].astype("float64")          # an inner window with no gaps
transform = rasterio.windows.transform(inner, bng_transform)
assert np.isfinite(dem).all()

def lidar(coverage_slug, coverage_id, path):
    r = requests.get(f"https://environment.data.gov.uk/spatialdata/{coverage_slug}/wcs", timeout=300,
                     params=[("service", "WCS"), ("version", "2.0.1"), ("request", "GetCoverage"),
                             ("CoverageId", coverage_id), ("format", "image/tiff"),
                             ("subset", "E(541000,543000)"), ("subset", "N(197000,199000)"),
                             ("SCALEFACTOR", "0.2")])
    r.raise_for_status()
    with open(path, "wb") as f:
        f.write(r.content)
lidar("lidar-composite-digital-surface-model-first-return-dsm-1m",
      "df4e3ec3-315e-48aa-aaaf-b5ae74d7b2bb__Lidar_Composite_Elevation_FZ_DSM_1m", "epping_dsm.tif")
lidar("lidar-composite-digital-terrain-model-dtm-1m",
      "13787b9a-26a4-4775-8523-806d13af58fc__Lidar_Composite_Elevation_DTM_1m", "epping_dtm.tif")


def levels_for(z, interval):
    lo, hi = np.nanmin(z), np.nanmax(z)
    return np.arange(np.floor(lo / interval) * interval, np.ceil(hi / interval) * interval + interval, interval)


def extract(surface, levels, transform, crs, min_vertices=0):
    """Contour a surface with matplotlib and return georeferenced lines."""
    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) < max(min_vertices, 2):
                continue
            xs, ys = xy(transform, seg[:, 1], seg[:, 0])
            records.append({"elev": float(level), "vertices": len(seg), "geometry": LineString(zip(xs, ys))})
    plt.close(fig)
    return gpd.GeoDataFrame(records, columns=["elev", "vertices", "geometry"], geometry="geometry", crs=crs)


print(f"Snowdon DEM in British National Grid: {dem.shape[1]} × {dem.shape[0]} at 30 m, "
      f"{dem.min():.0f}–{dem.max():.0f} m; {np.isnan(dem_full).sum():,} empty corner cells in the full grid")

Snowdon DEM in British National Grid: 328 × 304 at 30 m, 52–1071 m; 7,228 empty corner cells in the full grid

Example 1 — a diagnostic for a noisy contour layer

import pandas as pd


def diagnose_contours(lines, dem, *, cell_size=30.0, interval=None, dem_accuracy=4.0):
    problems = []
    vertices = lines["vertices"]
    lengths = lines.length

    tiny = (vertices < 8).mean()
    if tiny > 0.05:
        problems.append(f"{tiny:.0%} of lines have under 8 vertices — "
                        f"single-cell noise; smooth the DEM")

    short = (lengths < 4 * cell_size).mean()
    if short > 0.2:
        problems.append(f"{short:.0%} of lines are under 4 cells long — "
                        f"either canopy (use a DTM) or too fine an interval")

    if interval is not None:
        if interval < 2 * dem_accuracy:
            problems.append(f"interval {interval} m is under 2x the DEM accuracy "
                            f"({dem_accuracy} m) — the lines describe noise")
        n_levels = lines["elev"].nunique()
        if n_levels > 40:
            problems.append(f"{n_levels} levels — over about 30 is unreadable on a map")

    if np.issubdtype(dem.dtype, np.integer):
        problems.append(f"DEM dtype is {dem.dtype} — quantised to 1 m; "
                        f"expect staircase lines at fine intervals")

    print(f"  {len(lines)} lines, {vertices.sum():,} vertices, "
          f"median {vertices.median():.0f} per line")
    print(f"  length: median {lengths.median():,.0f} m, "
          f"shortest {lengths.min():,.0f} m, longest {lengths.max():,.0f} m")
    for problem in problems:
        print(f"  ✗ {problem}")
    if not problems:
        print("  ✓ contour layer looks clean")
    return problems


diagnose_contours(raw_lines, dem, cell_size=30.0, interval=50)
with rasterio.open("epping_dsm.tif") as src:
    canopy = src.read(1, masked=True).filled(np.nan)
    diagnose_contours(extract(canopy, levels_for(canopy, 5), src.transform, crs), canopy,
                      cell_size=5.0, interval=5, dem_accuracy=0.15)

  117 lines, 27,015 vertices, median 67 per line
  length: median 1,539 m, shortest 9 m, longest 44,226 m
  ✗ 15% of lines have under 8 vertices — single-cell noise; smooth the DEM
  16878 lines, 218,942 vertices, median 5 per line
  length: median 16 m, shortest 0 m, longest 25,674 m
  ✗ 67% of lines have under 8 vertices — single-cell noise; smooth the DEM
  ✗ 59% of lines are under 4 cells long — either canopy (use a DTM) or too fine an interval

One finding on the mountain DEM, and it names the fix. On the Epping surface model the canopy check fires as well.

Example 2 — choosing the smoothing by an error budget

def choose_smoothing(dem, *, interval=50, budget=None, dem_accuracy=4.0,
                     candidates=(0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0)):
    """Pick the largest sigma whose median contour error stays inside the budget."""
    budget = budget if budget is not None else dem_accuracy / 2
    rows = []
    for sigma in candidates:
        surface = gaussian_filter(dem, sigma) if sigma else dem
        kept, tiny, e = contour_and_check(surface, dem, interval)
        rows.append({"sigma": sigma, "lines": kept, "tiny": tiny,
                     "median_err": round(float(np.median(e)), 2),
                     "p95_err": round(float(np.percentile(e, 95)), 2)})
    frame = pd.DataFrame(rows)
    within = frame[frame["median_err"] <= budget]
    best = within["sigma"].max() if len(within) else 0.0
    print(frame.to_string(index=False))
    print(f"\nbudget {budget} m -> use sigma {best}")
    return float(best)


choose_smoothing(dem, interval=50)

 sigma  lines  tiny  median_err  p95_err
   0.0     99    20        0.00     0.00
   0.5     97    10        0.31     1.41
   1.0     94     8        1.19     5.35
   1.5     86     8        2.17     9.64
   2.0     83     5        3.24    14.13
   3.0     78     5        5.62    23.98
   4.0     73     5        8.22    33.64

budget 2.0 m -> use sigma 1.0

An explicit budget turns "how much should I smooth" from taste into arithmetic. Half the DEM's vertical accuracy is a defensible default: the contour's own displacement then contributes less than the DEM's uncertainty already does.

Example 3 — a full clean-up pass

def clean_contours(dem, transform, crs, *, interval=50, sigma=None,
                   min_vertices=8, min_length_cells=4, cell_size=30.0):
    sigma = sigma if sigma is not None else choose_smoothing(dem, interval=interval)
    levels = levels_for(dem, interval)
    gaps = np.isnan(dem)
    filled = np.where(gaps, np.nanmean(dem), dem)               # only so the filter has numbers
    smoothed = gaussian_filter(filled, sigma) if sigma else filled
    surface = np.ma.masked_array(smoothed, mask=gaps)           # then hide the gaps again

    lines = extract(surface, levels, transform, crs)
    dropped_tiny = int((lines["vertices"] < min_vertices).sum())
    lines = lines[lines["vertices"] >= min_vertices]
    keep = lines.length > min_length_cells * cell_size
    dropped_short = int((~keep).sum())
    lines = lines[keep].copy()
    lines["index"] = lines["elev"] % (interval * 5) == 0

    print(f"  sigma {sigma}, interval {interval} m")
    print(f"  {len(lines)} lines, {dropped_tiny} tiny dropped, {dropped_short} short dropped")
    print(f"  {lines['vertices'].sum():,} vertices total")
    return lines


lines = clean_contours(dem, transform, crs, interval=50, sigma=1.0)
lines.to_file("contours_clean.gpkg", layer="contours", driver="GPKG")

  sigma 1.0, interval 50 m
  94 lines, 6 tiny dropped, 0 short dropped
  25,682 vertices total

The feature count falls much more than the vertex count, because the tiny rings held few vertices each; the feature-count saving is what makes the layer usable. The median line position stays within step 1's measured error of the true surface.

The order matters: fill, smooth, mask the gaps again, extract, then filter. Filtering before extraction is impossible, and smoothing after extraction moves lines off their own levels.

Single-cell rings, integer steps and filled gaps

Why single-cell noise makes rings

Contouring finds where the surface crosses a level, cell by cell. A cell whose elevation sits just above a level, surrounded by cells just below it, is a closed crossing — so the algorithm draws a small ring around it.

With a 4 m vertical accuracy and 50 m contour intervals, about one cell in seven sits within 4 m of some level — 15% in the Snowdon window of the contours how-to. Most of them are next to other cells on the same side and produce nothing; the isolated ones produce rings.

The measured result on this DEM was 18 lines with fewer than eight vertices among 117. Smoothing removes them by removing the isolated highs and lows that create them.

Why smoothing the lines is worse than smoothing the surface

It is tempting to run simplify() or a spline over noisy contour lines. It produces smoother geometry and breaks the contours' defining property.

A contour separates ground above its level from ground below. Move the line and it can cross a cell that contradicts it — and two adjacent contours moved independently can cross each other, which is geometrically impossible for real contours and immediately visible on a map.

Smoothing the surface keeps the guarantee: the lines are exact contours of a slightly generalised surface. Example 2 measures how far that generalised surface has moved from the original.

Why an integer DEM staircases

An int16 DEM stores whole metres. The surface is therefore a set of flat plateaux with 1 m steps between them, not a continuous surface.

Contour at 1 m intervals and every level coincides with a step, so the algorithm has no gradient to interpolate along — the crossing is exactly at the cell boundary, and the line follows the cell edges.

Casting to float does not add information but does let a small Gaussian filter create gradients within the plateaux, which the contouring can then interpolate along. Alternatively use an interval of 5 m or more, so each contour crosses several quantisation steps.

Why any fill value draws the gap

Contouring finds where the surface crosses a level. Fill a gap with a single value and the gap becomes a flat plateau at that height, with a cliff all round its edge wherever the real ground is higher or lower. Every contour level between the plateau and the ground at the edge crosses that cliff, so the algorithm draws each of them around the gap.

The mean puts the plateau in the middle of the relief, so the levels near the mean trace the gap. Filling below the lowest level is the advice you will often read, and step 5 shows why it is wrong: the plateau is now below every level, so every level lower than the ground at the edge runs along the boundary — twice the length of the mean fill on this DEM.

A masked array has no plateau and no cliff. The contouring treats masked cells as absent and ends each line where the data ends. The same logic applies to sea encoded as 0: contour it as data and the 0 m and every low level trace the NoData boundary, not the coast.

Crossing contours, vertex counts and woodland texture

• Contours cannot cross. If yours do, they were smoothed after extraction.
• Measure the smoothing error by sampling the unsmoothed DEM at the contour vertices. Half the DEM's vertical accuracy is a defensible budget.
• Do not fill gaps before contouring — with the mean or with a value below every level. Contour a masked array.
• A closed contour may be a peak or a hollow. The geometry does not say; add spot heights or a hillshade.
• Integer DEMs staircase at intervals near the quantum. Cast to float and smooth, or widen the interval.
• gdal_contour produces the same noise — the fix is the same, applied to the input.
• Vertex counts scale badly. A 1 m interval on a 30 m DEM over a large area can produce millions of vertices and a GeoPackage nobody can draw.
• Woodland fingerprints mean a DSM. No parameter fixes it.

Internal links

• How to generate contour lines from a DEM in Python — the working pipeline
• Digital elevation models explained — vertical accuracy, and DSM canopy
• Slope values are wrong or absurdly steep — the same DEM problems in a different product
• How to simplify geometry in Python with GeoPandas and Shapely — and why not to use it on contours
• How to make a hillshade from a DEM in Python — for spotting DEM artefacts before contouring
• Topological Coverage: What "Clean" Means for Adjacent Polygons
• How to reduce GIS file size in Python without wrecking the data — after the vertex count explodes
• The raster data model explained — dtype and NoData

FAQ

Why are there hundreds of tiny circles in my contours?

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

How much should I smooth?

Enough that the median contour position error stays below about half the DEM's vertical accuracy. On a Copernicus DEM that is sigma 1, which here removes 12 of the 20 tiny rings for a 1.19 m median error.49 m median error.

Can I smooth the contour lines instead?

No. Moving a line can put it on the wrong side of cells it is supposed to separate, and independently smoothed neighbours can cross — which real contours never do.

Why do my contours look like staircases?

The DEM is quantised, usually int16 with 1 m steps, and your interval is near the quantum. Cast to float and smooth slightly, or widen the interval.

Why do contours run around the edge of a gap?

The gap was filled with a single value — the mean, or a value below every level — and the contouring drew the cliff at its edge. Contour a masked array (np.ma.masked_invalid) instead, and the lines stop where the data stops.

Why does woodland produce a fingerprint pattern?

You are contouring a DSM, so the lines follow the canopy. Use a bare-earth DTM; no smoothing parameter will fix it.

My contour file is enormous. What do I do?

Widen the interval — vertex count scales roughly with the number of levels — then drop short lines and tiny rings.