Smoothing a LiDAR DTM Before Contouring

TL;DR: Apply a NoData-aware Gaussian filter with sigma of 1–2 cells to a 1 m DTM, measure the change against the original (the 95th percentile of |Δz| should stay well below half the contour interval), and increase sigma only if contours still zigzag or form tiny loops. Use a median filter instead where spikes remain, and mask breaklines such as road edges and river banks if they must stay sharp.

# Context and Motivation

This guide is part of Contour Generation from LiDAR DTMs. A 1 m LiDAR DTM records relief at a scale no contour map intends to show: plough furrows, tyre ruts, residual low vegetation, interpolation texture. Contouring it directly draws every one of those, so lines zigzag and flat fields fill with small closed loops. Smoothing removes relief below a chosen scale before contouring. Done well, it changes the surface by a few centimetres and makes the contours readable; done badly, it rounds ridges, fills ditches and moves lines by more than the data’s own accuracy. The difference is measurement.

# Prerequisites and Assumptions

• A void-filled DTM GeoTIFF with NoData set.
• Python with rasterio, NumPy and SciPy.
• The planned contour interval, which sets how much smoothing change is acceptable.
• Optionally, breakline polygons or a mask of features that must not be smoothed.

# Step-by-Step Implementation

# Step 1 — Smooth with NoData awareness

Filter values and a validity mask separately and divide, so NoData neither spreads nor pulls edge values toward the fill value.

# Step 2 — Measure the change

Compute Δz = smoothed − original on valid cells and summarize: median, 95th and 99th percentiles of |Δz|.

# Step 3 — Compare with the interval

A reasonable ceiling: 95th percentile of |Δz| under a quarter of the contour interval. Beyond that, smoothing is moving contours visibly.

# Step 4 — Sweep sigma

Try 0.5, 1, 1.5, 2, 3 cells; contour each and count closed loops shorter than a threshold. Choose the smallest sigma at which loops drop off.

# Step 5 — Protect breaklines if required

Blend the original back in along masked features: out = where(mask, original, smoothed), with a feathered mask to avoid steps.

# Complete Working Example

"""Sweep Gaussian smoothing of a DTM and report surface change and contour clutter."""
from __future__ import annotations

import numpy as np
import rasterio
from scipy import ndimage as ndi
from skimage import measure

DTM = "mosaic/county_north_dtm_1m.tif"
INTERVAL = 0.5


def smooth(z: np.ndarray, valid: np.ndarray, sigma: float) -> np.ndarray:
    num = ndi.gaussian_filter(np.where(valid, z, 0.0), sigma)
    den = ndi.gaussian_filter(valid.astype(float), sigma)
    return np.where(valid, num / np.maximum(den, 1e-6), np.nan)


def small_loops(z: np.ndarray, interval: float, max_len_cells: float = 40) -> int:
    levels = np.arange(np.nanmin(z) // interval * interval + interval, np.nanmax(z), interval)
    zz = np.where(np.isnan(z), np.nanmin(z) - 1000, z)
    count = 0
    for lv in levels:
        for c in measure.find_contours(zz, lv):
            closed = np.allclose(c[0], c[-1])
            length = np.sum(np.hypot(*np.diff(c, axis=0).T))
            count += int(closed and length < max_len_cells)
    return count


with rasterio.open(DTM) as ds:
    window = rasterio.windows.Window(0, 0, 2000, 2000)       # a 2 km test block
    z = ds.read(1, window=window, masked=True).filled(np.nan).astype(float)
valid = ~np.isnan(z)

print(f"{'sigma':>6} {'p50|dz|':>9} {'p95|dz|':>9} {'p99|dz|':>9} {'small loops':>12}")
for sigma in (0.0, 0.5, 1.0, 1.5, 2.0, 3.0):
    s = z if sigma == 0 else smooth(z, valid, sigma)
    dz = np.abs(s - z)[valid]
    print(f"{sigma:>6} {np.percentile(dz, 50):>9.3f} {np.percentile(dz, 95):>9.3f} "
          f"{np.percentile(dz, 99):>9.3f} {small_loops(s, INTERVAL):>12}")

Illustrative output for rolling farmland at 1 m:

 sigma   p50|dz|   p95|dz|   p99|dz|  small loops
   0.0     0.000     0.000     0.000         4118
   0.5     0.008     0.031     0.058         1307
   1.0     0.015     0.052     0.094          402
   1.5     0.021     0.071     0.131          118
   2.0     0.027     0.089     0.170           61
   3.0     0.038     0.121     0.244           39

At 0.5 m contours, sigma 1.5 keeps the 95th-percentile change at 7 cm — well under a quarter of the interval — while removing 97 percent of the tiny loops.

# Key Parameter Table

# Verification

• Change statistics within the acceptance rule for the chosen interval.
• Visual comparison. Overlay contours from original and smoothed DTMs; lines should shift slightly and lose zigzags, not move across features.
• Features preserved. Profile across a known ditch and ridge; depth and height should change by no more than a few centimetres.

# Gotchas and Edge Cases

Smoothing across NoData. A plain gaussian_filter on a raster with −9999 fill values drags edge cells down by hundreds of metres. Always use the normalized form or fill voids first.

Smoothing hydro-flattened water. Smoothing blurs the sharp shoreline of flattened lakes. Mask water polygons out of the smoothing and paste flattened values back afterwards.

Units of sigma. Sigma is in cells. On a 0.5 m DTM, sigma 1.5 cells is 0.75 m — half the physical smoothing of the same sigma on a 1 m DTM. Express the choice in metres when comparing projects.

Over-smoothing to hide classification errors. Heavy smoothing can make bumps from misclassified vegetation disappear from contours, but they remain in the DTM. Fix classification rather than masking it.

# Frequently Asked Questions

How much should I smooth a DTM before contouring?

Enough to remove micro-relief below the scale the contours represent, but not so much that the surface moves by more than about a quarter of the contour interval. For a 1 metre DTM and 0.5 metre contours, a Gaussian sigma of 1 to 2 cells is typical.

Gaussian or median filter?

Gaussian for general micro-relief such as furrows and interpolation texture; median where isolated spikes or pits remain. Median filters preserve edges better but can look blocky on smooth slopes.

Does smoothing change the accuracy of the DTM?

It changes the surface slightly, so it should be applied to a contouring copy, not the delivered DTM. Measure the change and keep it well below the vertical accuracy.

How do I avoid smoothing across NoData?

Filter the values with NoData set to zero and a separate validity mask, then divide the two results. This normalized filter keeps edge values correct and leaves voids untouched.

# Related

• Contour Generation from LiDAR DTMs — the full workflow
• Generating Contours from a DTM with gdal_contour — contouring the smoothed surface
• Exporting Contours to GeoPackage — cleaning the result
• Filling NoData Voids in DTM Rasters — preparing the input
• Removing Pits and Spikes from a DSM — median filtering for surfaces