Generating Smooth Contours Without Staircase Artifacts

The contour set comes back looking like a pixel-art staircase: long runs along cell edges joined by right-angle turns, every line following the raster grid rather than the terrain. It is unmistakably wrong and the instinct — run a line-smoothing filter over the geometry — makes the drawing look correct while moving the contours off the surface they claim to describe.

The staircase is a symptom of the DEM, not of the contour generator. This page identifies which of the three causes produced it and fixes each at the source.

Three ways a DEM produces stepped contours

Elevation quantisation. The DEM is stored as an integer type, or as a float that was rounded, so elevations take discrete values. A contour at 118.0 m then follows the boundary between cells valued 118 and 119, which is a cell edge — so the line is literally a staircase because the surface is.

Nearest-neighbour resampling. A DEM reprojected or resampled with nearest-neighbour carries each output cell the value of one input cell. Every cell in a block of four then shares a value, and the contour traces the block boundary.

A grid coarser than the contour interval implies. With a 50 cm interval on a 1 m grid over gentle terrain, consecutive cells often differ by less than the interval, so a contour crosses whole runs of identical cells and follows their edges.

All three produce the same appearance and have the same cure: the surface must vary smoothly at the scale the contour interval samples.

A quantised surface and the contour it produces Two elevation profiles across the same slope with the contour each produces. In the first, elevations are quantised to whole metres, so the profile is a staircase of flat treads and the contour at a given level runs along a cell edge for the whole width of a tread before stepping. In the second, elevations vary continuously, so the profile is a smooth ramp and the contour crosses each cell at an interpolated position, producing a smooth line. A note states that the contour generator interpolates within cells in both cases, and that it is the surface rather than the interpolation that differs. quantised elevations contour level the contour follows cell edges continuous elevations contour level the contour crosses cells at interpolated points The generator interpolates within cells in both cases; only the surface differs.

Figure 1 — The staircase is in the data. A contour generator can only interpolate between the values it is given, and between two identical values there is nothing to interpolate.

Minimal reproducible solution

Diagnose first, because the three causes have different fixes and one of them is not a smoothing problem at all.

import numpy as np
import rasterio


def diagnose_quantisation(dem_path: str) -> str:
    """Name the cause from the DEM's values, not from the contour geometry."""
    with rasterio.open(dem_path) as ds:
        band = ds.read(1, masked=True)
        dtype = ds.dtypes[0]
        cell = abs(ds.transform.a)

    vals = band.compressed()
    if vals.size < 1000:
        return "too little valid data to diagnose"

    if dtype.startswith(("int", "uint")):
        return f"integer dtype {dtype} — elevations are quantised to 1 unit"

    # Distinct values relative to sample size: heavy repetition means either
    # rounding or nearest-neighbour resampling.
    sample = vals[:: max(1, vals.size // 200_000)]
    distinct_ratio = np.unique(sample).size / sample.size
    if distinct_ratio < 0.02:
        return (f"only {distinct_ratio:.1%} distinct values — rounded elevations "
                "or nearest-neighbour resampling")

    # Local relief per cell against the contour interval matters too.
    gy, gx = np.gradient(band.filled(np.nan))
    step = float(np.nanmedian(np.hypot(gx, gy)))
    return f"continuous; median elevation change per cell is {step:.3f} m at {cell:.2f} m cells"

The distinct-value ratio is the diagnostic that separates rounding from nearest-neighbour resampling in practice. A genuinely continuous float DEM has almost as many distinct values as cells; anything below a couple of percent means values are being repeated at scale, and repeated values are what a staircase is made of.

The fix is at the surface, in this order.

import numpy as np
import rasterio
from scipy.ndimage import gaussian_filter


def prepare_for_contours(src_path: str, dst_path: str,
                         sigma_cells: float = 1.0) -> None:
    """Write a float32 DEM, smoothed enough that the interval is resolvable.

    Casting to float32 removes integer quantisation; the smoothing then
    removes the residual terracing that rounding left behind.
    """
    with rasterio.open(src_path) as src:
        band = src.read(1, masked=True)
        profile = src.profile

    arr = band.filled(np.nan).astype(np.float32)
    valid = np.isfinite(arr)

    num = gaussian_filter(np.where(valid, arr, 0.0), sigma_cells)
    den = gaussian_filter(valid.astype(np.float32), sigma_cells)
    out = np.where(den > 1e-6, num / np.maximum(den, 1e-6), np.nan)
    out[~valid] = np.nan

    profile.update(dtype="float32", nodata=np.nan, compress="deflate", predictor=3)
    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(out, 1)
        dst.update_tags(CONTOUR_PREP_SIGMA_CELLS=str(sigma_cells))

Recording the sigma in the tags is what makes a later revision reproducible, and it is the same value the contour and hillshade workflow uses for every other derived product — so the contours and the slope raster describe one surface.

Edge-case matrix

Cause Diagnostic Fix
Integer DEM dtype is int or uint Cast to float32, then smooth
Rounded floats very few distinct values Smooth; the information is gone
Nearest-neighbour resampling values repeat in blocks Resample again from the source, bilinear
Interval finer than the relief continuous, small per-cell change Widen the interval
Interval finer than the accuracy continuous Widen; see the interval rule
Smoothing applied to the lines geometry departs from the DEM Undo; smooth the surface instead
Voids filled with a constant flat plateaus, ringed by contours Fill properly, or contour unfilled
Contours from a hillshade nonsense Contour the DEM, not a rendering

The nearest-neighbour row is the one worth going back for. If the DEM was resampled that way, smoothing recovers a plausible surface and not the original one — the information was discarded at the resample. Going back to the source raster and resampling bilinearly costs one command and produces a genuinely better DEM.

Verify the fix worked

Two things must hold: the lines must be smooth, and they must still lie on the surface.

import numpy as np
import fiona
import rasterio


def contour_turn_angles(gpkg: str, layer: str | None = None) -> np.ndarray:
    """Interior angles at contour vertices, in degrees."""
    angles = []
    with fiona.open(gpkg, layer=layer) as src:
        for feat in src:
            pts = np.asarray(feat["geometry"]["coordinates"], dtype=float)[:, :2]
            if len(pts) < 3:
                continue
            v1, v2 = pts[1:-1] - pts[:-2], pts[2:] - pts[1:-1]
            cos = np.einsum("ij,ij->i", v1, v2) / (
                np.linalg.norm(v1, axis=1) * np.linalg.norm(v2, axis=1) + 1e-12)
            angles.append(np.degrees(np.arccos(np.clip(cos, -1, 1))))
    return np.concatenate(angles) if angles else np.array([])


def assert_smooth_and_on_surface(gpkg: str, dem_path: str,
                                 max_right_angle_frac: float = 0.02,
                                 tol_m: float = 0.02) -> None:
    """Smooth lines that still sit at their stated elevation."""
    ang = contour_turn_angles(gpkg)
    right = float(np.mean(np.abs(ang - 90.0) < 5.0))
    assert right <= max_right_angle_frac, (
        f"{right:.1%} of vertices turn through ~90° — the staircase remains")

    worst = 0.0
    with rasterio.open(dem_path) as dem, fiona.open(gpkg) as src:
        for feat in src:
            elev = float(feat["properties"]["elev"])
            pts = [tuple(p[:2]) for p in feat["geometry"]["coordinates"][::10]]
            for v in dem.sample(pts, indexes=1):
                if np.isfinite(v[0]):
                    worst = max(worst, abs(float(v[0]) - elev))
    assert worst <= tol_m, (
        f"contours depart from the DEM by up to {worst:.3f} m — the geometry "
        "was smoothed rather than the surface")

Running both assertions together is the point. Line smoothing passes the first and fails the second; leaving the DEM quantised passes the second and fails the first. Only fixing the surface passes both.

Three approaches against the two tests A comparison table of three approaches against two acceptance tests. Leaving the quantised DEM alone fails the smoothness test and passes the on-surface test. Smoothing the contour geometry after generation passes the smoothness test and fails the on-surface test, because the lines have been moved away from the elevations they claim. Smoothing the DEM before generating passes both. A note observes that each single test can be satisfied by an approach that is wrong, which is why both are needed. smooth geometry? on the surface? leave the DEM quantised smooth the contour lines smooth the DEM first failspasses passesfails passespasses Each test alone can be satisfied by an approach that is wrong.

Figure 2 — Why both assertions are needed. The tempting fix passes exactly the test that was being looked at, which is what makes it tempting.

When to escalate

  • The DEM is float, has plenty of distinct values, and contours still step. The interval is finer than the terrain’s relief per cell, so consecutive cells straddle several intervals in some places and none in others. Widen the interval rather than smoothing further; smoothing at that point starts erasing real features.
  • Smoothing enough to fix the steps erases a kerb or a bund. The interval is too fine for the data. This is the accuracy question rather than a rendering one, and it is resolved by the two-to-three-times-RMSE rule in the contour workflow.
  • Contours are smooth in the open and stepped over one region. That region was filled or resampled differently — usually a void patched with a constant. Check the fill mask before touching the smoothing, since the plateau is data, not quantisation.

Finally, resist the temptation to fix this at the drawing stage even when it is quicker. A contour set that has been geometrically smoothed carries no record of the fact, so the next person to sample the DEM along those lines — for a profile, a volume, or a cross-section — gets values that disagree with the contour labels and has no way to know why. Smoothing the surface leaves the two products consistent by construction, and the sigma recorded in the raster tags says exactly what was done.

← Deriving Contours and Hillshade with GDAL

Choosing the smoothing sigma against what it erases A chart of two quantities against the Gaussian smoothing sigma in cells. The fraction of contour vertices turning through approximately ninety degrees falls steeply as sigma rises, reaching near zero by about one cell. The smallest terrain feature that survives grows steadily with sigma, passing a typical kerb height at about one and a half cells and a bund at about three. A shaded band between one and one and a half cells is marked as the working range where the staircase is gone and real features are retained. working range 0 1.0 1.5 2.5 4.0 Gaussian sigma, in cells right-angle vertices smallest surviving feature Below the band the staircase remains; above it a kerb has been smoothed away and the DEM no longer shows it.

Figure 3 — The sigma is a survey decision. One cell removes the artefact; three removes features the client is paying to see, and nothing in the contour geometry distinguishes the two.