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 of 400,000 vertices for an area that needs 20,000.
• 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:

import numpy as np

vertices = lines.geometry.apply(lambda g: len(g.coords))
print(f"{len(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}")

103 lines, 19,090 vertices
tiny rings (<8 verts): 20 (19%)
median vertices per line: 34

Step-by-step solution #

1. Smooth the surface, and measure what it costs #

from scipy.ndimage import gaussian_filter, map_coordinates

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

sigma 0.0:  83 lines, 20 tiny  |  |error| median  0.00 m  max   0.00 m
sigma 1.0:  79 lines,  8 tiny  |  |error| median  1.49 m  max  18.63 m
sigma 2.0:  70 lines,  7 tiny  |  |error| median  3.87 m  max  39.82 m
sigma 4.0:  63 lines,  3 tiny  |  |error| median  9.49 m  max  87.44 m

The error is measured by sampling the unsmoothed DEM where each contour claims to be. Sigma 1 removes 60% of the tiny rings for a median error of 1.49 m — inside Copernicus DEM's own 4 m accuracy, so the lines are still honest.

Sigma 4 removes almost all of them and puts the median error at 9.49 m, with a maximum of 87 m. A 50 m contour that is 87 m out 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 = lines[lines.geometry.apply(lambda g: len(g.coords)) >= 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:

min_length = 4 * max(cell_x, cell_y)          # four cells around
lines = lines[lines.length > min_length]
print(f"kept {len(lines)} lines longer than {min_length:.0f} 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:

unique = np.unique(dem[np.isfinite(dem)])
steps = np.diff(unique)
print(f"dtype {dem.dtype}, {len(unique):,} distinct values")
print(f"smallest step {steps.min():.4f} m, most common {np.bincount((steps * 1000).astype(int)).argmax() / 1000:.3f} m")

dtype int16, 1,021 distinct values
smallest step 1.0000 m, most common 1.000 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 #

short = lines[lines.length < 200]
print(f"{len(short)} lines shorter than 200 m ({len(short) / len(lines):.0%})")
print(f"they hold {short.geometry.apply(lambda g: len(g.coords)).sum():,} vertices")

If most of your lines are short and they cluster in woodland and around buildings, you are contouring a canopy. 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. Fix lines that stop in the middle of the map #

print(f"NaN cells: {np.isnan(dem).sum():,}")
print(f"fill value used: {fill}")
print(f"lowest contour level: {levels.min()}")

If NaN was filled with the mean elevation, contours at levels near the mean will draw around the void boundary, producing a line that stops abruptly or encircles a hole.

Fill below the lowest level instead:

filled = np.where(np.isnan(dem), levels.min() - interval, dem)

Now no contour level intersects the fill value, so the lines simply stop at the void edge.

Code examples #

Example 1 — a diagnostic for a noisy contour layer #

import numpy as np
import pandas as pd


def diagnose_contours(lines, dem, *, cell_size=30.0, interval=None, dem_accuracy=4.0):
    problems = []
    vertices = lines.geometry.apply(lambda g: len(g.coords))
    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")

    finite = dem[np.isfinite(dem)]
    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=28.0, interval=50)

  103 lines, 19,090 vertices, median 34 per line
  length: median 1,847 m, shortest 62 m, longest 24,193 m
  ✗ 19% of lines have under 8 vertices — single-cell noise; smooth the DEM

One finding, and it names the fix. Run it on the DSM version of the same area and the second and third checks fire too.

Example 2 — choosing the smoothing by an error budget #

from scipy.ndimage import gaussian_filter, map_coordinates


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
    lo, hi = np.nanmin(dem), np.nanmax(dem)
    levels = np.arange(np.floor(lo / interval) * interval,
                       np.ceil(hi / interval) * interval + interval, interval)

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

    rows = []
    for sigma in candidates:
        surface = gaussian_filter(dem, sigma) if sigma else dem
        fig, ax = plt.subplots()
        cs = ax.contour(surface, levels=levels)
        errors, kept, tiny = [], 0, 0
        for level, segments in zip(cs.levels, cs.allsegs):
            for seg in segments:
                if len(seg) < 8:
                    tiny += 1
                    continue
                kept += 1
                errors.append(np.abs(
                    map_coordinates(dem, [seg[:, 1], seg[:, 0]], order=1, mode="nearest")
                    - level))
        plt.close(fig)
        e = np.concatenate(errors)
        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     83    20        0.00     0.00
   0.5     84     9        0.42     1.72
   1.0     79     8        1.49     6.25
   1.5     76     4        2.63    10.86
   2.0     70     7        3.87    15.78
   3.0     65     3        6.58    25.79
   4.0     63     3        9.49    36.55

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)

    lo, hi = np.nanmin(dem), np.nanmax(dem)
    levels = np.arange(np.floor(lo / interval) * interval,
                       np.ceil(hi / interval) * interval + interval, interval)

    filled = np.where(np.isnan(dem), levels.min() - interval, dem)   # below every level
    surface = gaussian_filter(filled, sigma) if sigma else filled

    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from shapely.geometry import LineString
    import geopandas as gpd
    import rasterio

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

    lines = gpd.GeoDataFrame(records, crs=crs)
    if crs and not crs.is_geographic:
        keep = lines.length > min_length_cells * cell_size
        dropped_short = int((~keep).sum())
        lines = lines[keep]

    lines["index"] = lines["elev"] % (interval * 5) == 0
    print(f"  sigma {sigma}, interval {interval} m")
    print(f"  {len(lines)} lines, {dropped_tiny} tiny dropped, "
          f"{dropped_short} short dropped")
    print(f"  {lines.geometry.apply(lambda g: len(g.coords)).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
  79 lines, 8 tiny dropped, 0 short dropped
  18,077 vertices total

From 103 lines and 19,090 vertices to 79 lines and 18,077 — a 23% reduction in features and 5% in vertices, with the median line position 1.49 m from the true surface. The vertex saving is modest because the tiny rings held few vertices each; the feature-count saving is what makes the layer usable.

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

Explanation #

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, roughly 8% of cells sit within noise distance of some level. 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 20 rings from 103 segments. 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 the fill value matters #

Contouring cannot handle NaN, so it must be filled. Filling with the mean puts a large flat region at a value that some contour level will pass through — and the algorithm dutifully draws that level around the void boundary.

The result is a contour that appears to describe terrain and is actually tracing the edge of missing data. Filling below the lowest level guarantees no contour intersects it, so the lines simply terminate at the void.

The same logic applies to sea: if sea is encoded as 0 and your lowest contour is 0, you will get a contour around every coastline that is really the NoData boundary.

Edge cases or notes #

• 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.
• Fill NaN below the lowest level, never with the mean.
• 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
• Topology in GIS explained: slivers, gaps and shared boundaries — why crossing lines are invalid
• 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 removes 60% of the noise rings for a 1.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 stop in the middle of the map? #

NaN was filled with a value that a contour level passes through — usually the mean. Fill below the lowest level instead.

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.