Python LiDAR

Contour Generation from LiDAR DTMs

Contours are still one of the most requested LiDAR deliverables. Engineers design on them, planners read them, and many regulations and legacy CAD workflows expect them. A 1 m LiDAR DTM contains everything needed, but contouring it naively produces lines that zigzag around every cell, break into thousands of tiny closed loops on flat ground, and imply an accuracy the data does not have. This topic in the Ground Filtering and DTM/DSM Generation section covers the full workflow: picking an interval that the vertical accuracy supports, smoothing the surface so the lines represent landform rather than noise, generating contours with GDAL from Python, cleaning and attributing them, and delivering them in a form GIS and CAD users can load directly.

Contours from raw and smoothed DTMs Two panels of the same hillside. Contours from the raw 1 metre DTM are jagged with small spurious closed loops on the flatter ground at the bottom. Contours from a lightly smoothed DTM are clean, evenly spaced curves that follow the landform, with no spurious loops. raw 1 m DTM smoothed DTM zigzags and spurious loops on flat ground lines follow the landform

# Prerequisites

  • A bare-earth DTM GeoTIFF from class 2 returns, ideally hydro-flattened and void-filled; see generating a DTM GeoTIFF with writers.gdal and filling NoData voids.
  • GDAL 3.x Python bindings (from osgeo import gdal, ogr), plus rasterio, NumPy, SciPy, GeoPandas and Shapely.
  • The DTM’s vertical accuracy, from vertical accuracy assessment, to set a defensible interval.
  • A projected CRS in metres (or feet if the client works in feet), with the vertical datum known.
  • The client’s conventions: interval, index contour spacing, attribute names, and whether contours must be closed at tile edges.

# Core Workflow Architecture

  1. Choose the interval. A common rule of thumb is that the contour interval should be at least about twice the NVA at 95 % confidence, so that a contour’s position is meaningful. For QL2 data with NVA around 0.2 m, 0.5 m contours are defensible; 0.25 m contours are not.
  2. Mosaic and prepare. Contour a seamless mosaic, not individual tiles, to avoid broken lines at tile edges. Fill voids and apply hydro-flattening first.
  3. Smooth. Apply a Gaussian or median filter sized to a few cells to remove micro-relief that would produce zigzags, without shifting the landform.
  4. Contour. Run gdal.ContourGenerateEx with a fixed interval and base, writing an elevation attribute.
  5. Clean. Drop closed loops below a minimum length, simplify slightly, and tag index contours (every fifth line, typically).
  6. Deliver. Write to GeoPackage or Shapefile with elevation, index flag and CRS; optionally DXF for CAD.

# Full Implementation

python
"""DTM to clean, attributed contours in a GeoPackage."""
from __future__ import annotations

import logging
from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
from osgeo import gdal, ogr, osr
from scipy import ndimage as ndi

gdal.UseExceptions()
log = logging.getLogger("contours")


def smooth_dtm(src: Path, dst: Path, sigma_cells: float = 1.5) -> None:
    with rasterio.open(src) as ds:
        z = ds.read(1, masked=True)
        profile = ds.profile
    filled = z.filled(np.nan)
    valid = ~np.isnan(filled)
    # Normalized Gaussian: smooth values and weights separately so NoData does not bleed in.
    num = ndi.gaussian_filter(np.where(valid, filled, 0.0), sigma_cells)
    den = ndi.gaussian_filter(valid.astype(float), sigma_cells)
    out = np.where(valid, num / np.maximum(den, 1e-6), profile.get("nodata", -9999))
    profile.update(dtype="float32")
    with rasterio.open(dst, "w", **profile) as o:
        o.write(out.astype("float32"), 1)


def contour(dtm: Path, gpkg: Path, interval: float, base: float = 0.0) -> None:
    ds = gdal.Open(str(dtm))
    band = ds.GetRasterBand(1)
    srs = osr.SpatialReference(wkt=ds.GetProjection())
    drv = ogr.GetDriverByName("GPKG")
    if gpkg.exists():
        drv.DeleteDataSource(str(gpkg))
    out = drv.CreateDataSource(str(gpkg))
    layer = out.CreateLayer("contours_raw", srs, ogr.wkbLineString)
    layer.CreateField(ogr.FieldDefn("id", ogr.OFTInteger))
    layer.CreateField(ogr.FieldDefn("elev", ogr.OFTReal))
    nodata = band.GetNoDataValue()
    opts = [f"LEVEL_INTERVAL={interval}", f"LEVEL_BASE={base}", "ID_FIELD=0", "ELEV_FIELD=1"]
    if nodata is not None:
        opts.append(f"NODATA={nodata}")
    gdal.ContourGenerateEx(band, layer, options=opts)
    out = None


def clean(gpkg: Path, interval: float, index_every: int = 5,
          min_loop_m: float = 40.0, simplify_m: float = 0.25) -> gpd.GeoDataFrame:
    gdf = gpd.read_file(gpkg, layer="contours_raw")
    gdf = gdf.explode(index_parts=False).reset_index(drop=True)
    closed = gdf.geometry.is_ring
    gdf = gdf[~(closed & (gdf.length < min_loop_m))].copy()
    gdf["geometry"] = gdf.geometry.simplify(simplify_m, preserve_topology=True)
    gdf["elev"] = gdf["elev"].round(3)
    steps = np.round(gdf["elev"] / interval).astype(int)
    gdf["index"] = (steps % index_every == 0).astype(int)
    gdf = gdf[["elev", "index", "geometry"]]
    gdf.to_file(gpkg, layer="contours", driver="GPKG")
    log.info("%d contour lines, %d index", len(gdf), int(gdf["index"].sum()))
    return gdf


def run(dtm: Path, out_dir: Path, interval: float = 0.5) -> gpd.GeoDataFrame:
    out_dir.mkdir(parents=True, exist_ok=True)
    smooth = out_dir / "dtm_smooth.tif"
    gpkg = out_dir / "contours.gpkg"
    smooth_dtm(dtm, smooth)
    contour(smooth, gpkg, interval)
    return clean(gpkg, interval)


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    run(Path("mosaic/county_north_dtm_1m.tif"), Path("out/contours"), interval=0.5)

# Code Breakdown

Normalized Gaussian smoothing. A plain Gaussian filter on a raster with NoData either spreads NoData or pulls edge cells toward the fill value. Smoothing the values and a validity mask separately and dividing keeps edges correct and leaves NoData cells untouched.

sigma_cells: 1.5. On a 1 m DTM this removes variation at the scale of a couple of metres — ploughing furrows, small debris, residual vegetation bumps — while moving contour positions by only centimetres on real slopes. The dedicated guide on smoothing a LiDAR DTM before contouring quantifies the trade-off.

ContourGenerateEx with named options. LEVEL_INTERVAL and LEVEL_BASE define the elevations; ID_FIELD and ELEV_FIELD give the indexes of the fields created on the layer. Passing NODATA stops contours being drawn around the NoData boundary.

Explode, then filter loops. GDAL may emit multi-part lines. Exploding makes each part a row, so tiny closed rings — spurious loops on flat ground — can be dropped by length.

Topology-preserving simplification. A tolerance a quarter of a cell removes the stair-step vertices contouring produces on a grid without letting adjacent contours cross.

Index contours by integer step. Computing the step number from the elevation and interval avoids floating-point surprises (elev % 2.5 on doubles is unreliable).

# Parameter Reference Table

Parameter Type Default here Typical range Effect
interval float, m 0.5 0.25–5 Vertical spacing; ≥ about 2 × NVA95
base float, m 0.0 — Offset of the contour series
sigma_cells float 1.5 0.5–4 Smoothing strength before contouring
min_loop_m float, m 40 10–200 Shortest closed contour kept
simplify_m float, m 0.25 0.1–1.0 Vertex reduction tolerance
index_every int 5 4–10 Every nth contour flagged as index
Interval against accuracy A vertical scale showing a DTM's 95 percent vertical accuracy band of plus or minus 0.2 metres around a true surface. Contours at 0.25 metre spacing fall within each other's uncertainty bands and cannot be distinguished reliably. Contours at 0.5 metre spacing are separated by more than the band, so each line's position is meaningful. 0.25 m interval 0.5 m interval uncertainty bands overlap: lines not meaningful bands separated: each line is meaningful shaded: ±0.2 m (95 %) around each contour level

# Validation and Integrity Checks

  • No crossings. Contours of different elevations must never intersect. A spatial self-join of the output checks it; any hit points to a simplification tolerance that is too large.
  • Elevations on the series. Every elev value should be base + k × interval to rounding.
  • Continuity at mosaic joins. If you contoured tiles separately, check for dangling line ends along tile edges; the fix is to contour the mosaic.
  • Spot comparison. Sample the original (unsmoothed) DTM along a few contours; values should scatter tightly around the contour elevation, with a standard deviation well below the interval.
python
from shapely.strtree import STRtree
g = gpd.read_file("out/contours/contours.gpkg", layer="contours")
tree = STRtree(g.geometry.values)
pairs = tree.query(g.geometry.values, predicate="crosses")
bad = [(i, j) for i, j in zip(*pairs) if i < j and g.elev.iat[i] != g.elev.iat[j]]
assert not bad, f"{len(bad)} crossing contour pairs"

# Index Contours, Labels and Cartography

Contour sets are read by people, and a few cartographic conventions make them far easier to read. Index contours — every fifth line at a 0.5 m interval, so every 2.5 m — are drawn heavier and carry elevation labels; intermediate contours are thin and unlabelled. Labels sit along the line, oriented uphill so that the top of the text faces higher ground, and are placed on straight stretches away from dense clusters. On very flat ground, supplementary contours at half the interval, drawn dashed, can show subtle relief where standard contours are hundreds of metres apart.

Index contours carry the labels Five roughly parallel contour lines on a slope. The middle one is an index contour, drawn heavier, with the label 215.0 placed along it and oriented so the text's top faces uphill. The other four are thin intermediate contours at 0.5 metre spacing without labels. An arrow indicates the uphill direction. 215.0 uphill heavy: index contour, labelled · thin: intermediate, 0.5 m apart

None of this needs to be baked into the geometry. Deliver lines with an elev value and an index flag, and let the client’s GIS or CAD style them. What does need care in the data is line continuity: labelling engines place labels on long, unbroken lines, so contours split at every tile edge or simplified into fragments label poorly.

# Contours for CAD and Engineering Users

Many engineering users will load contours into CAD rather than GIS. Three details avoid friction there. First, CAD tools expect 3D polylines or 2D lines with elevation; writing each line’s Z coordinate equal to its elevation (a 3D line string) serves both. Second, units must match the drawing — US projects often work in US survey feet, so contour the DTM in the delivery units or convert elevations and interval together, never one without the other. Third, a DXF export with the elevation in a layer name or attribute follows most CAD conventions; GDAL’s DXF driver writes it from the same GeoDataFrame used for the GeoPackage.

Engineers also tend to ask for finer intervals than the data supports. The honest response is to deliver the defensible interval and, if needed, a separately labelled set of supplementary contours with an explicit note on accuracy, rather than a dense set that looks precise and is not.

# Performance Tuning

Contouring is fast; a 20 km × 20 km mosaic at 1 m takes a minute or two. The expensive parts are smoothing very large rasters in memory and writing hundreds of thousands of features.

  • Contour a VRT mosaic. gdal.BuildVRT over tiles gives GDAL a seamless virtual raster without copying data.
  • Smooth in blocks with overlap. For rasters too large for memory, smooth windows with a margin of 4 × sigma and write the interiors.
  • Coarsen for small-scale maps. 2 m or 5 m contours for a regional map can come from a DTM resampled to 2–5 m, which is much faster and smoother.

# Common Errors and Troubleshooting

Contours zigzag along cell edges. No smoothing, or simplification set to zero. Apply a small Gaussian and simplify at a quarter of a cell.

Thousands of tiny loops. Micro-relief on flat ground crosses contour levels. Smooth more, raise the interval, or drop short closed loops.

Contours around NoData holes. NoData was not passed to the contouring call, or voids were not filled. Fill voids first, or pass NODATA so GDAL skips them.

Lines stop at tile edges. Tiles were contoured one by one. Contour a mosaic instead.

Contours drawn across lakes. The DTM was not hydro-flattened; interpolated water surfaces produce meaningless lines. Flatten water first, as in hydro-flattening water bodies in a DTM.

# Frequently Asked Questions

What contour interval should I use for LiDAR?

One the vertical accuracy supports: roughly twice the 95 percent vertical accuracy or more. For typical QL2 data that means 0.5 metres or coarser; tighter intervals imply precision the data does not have.

Should I smooth the DTM before contouring?

Yes, lightly. A small Gaussian filter removes micro-relief that makes contours zigzag and form spurious loops, while moving real contours by only centimetres. Heavy smoothing flattens real landforms, so keep it modest and check the result.

How do I generate contours in Python?

Use GDAL’s ContourGenerateEx on the DTM band with LEVEL_INTERVAL, LEVEL_BASE and field options, writing into an OGR layer such as a GeoPackage. GeoPandas can then clean and attribute the lines.

Should I contour tiles or a mosaic?

A mosaic. Contouring tiles separately breaks every line at every tile edge and can produce slightly different line positions on either side of a seam. Build a virtual mosaic with gdal.BuildVRT and contour that once; if the area is too large, contour overlapping blocks and merge lines by elevation afterwards.

Why do contours cross after simplification?

Because the simplification tolerance was large compared with the spacing between adjacent contours on steep ground. Use topology-preserving simplification, keep the tolerance to a fraction of a cell, and check the output for crossings with a spatial join before delivery.

Can I produce contours in feet from a metre DTM?

Yes, but convert elevations first: multiply the DTM by the appropriate foot factor, or contour in metres with an interval equivalent to the desired feet and convert the elevation attribute afterwards. Be explicit about US survey feet versus international feet, and state the vertical datum in the metadata.

What are index contours?

Every fifth contour, typically, drawn heavier and labelled on maps to make elevations easy to read. Flag them with an attribute so cartographic styling can pick them out.

Hands-on guides

Exporting Contours to GeoPackage --- title: "Exporting Contours to GeoPackage" description: "Deliver LiDAR contours as a clean, well-documented GeoPackage: 3D line geometry... Smoothing a LiDAR DTM Before Contouring --- title: "Smoothing a LiDAR DTM Before Contouring" description: "Choose and tune the smoothing applied to a LiDAR DTM before contouring:... Generating Contours from a DTM with gdal_contour --- title: "Generating Contours from a DTM with gdalcontour" description: "Create contour lines and filled contour polygons from a LiDAR DTM...

More in Ground Filtering and DTM/DSM Generation with PDAL

Canopy Height Models with filters.hag_nn --- title: "Canopy Height Models with filters.hagnn" description: "Producing a canopy height model in the point domain:... CSF Cloth Simulation Ground Filtering --- title: "CSF Cloth Simulation Ground Filtering" description: "Classify ground with PDAL filters.csf: how the cloth... DTM Raster Generation with PDAL --- title: "DTM Raster Generation with PDAL" description: "Turning classified ground returns into a bare-earth Digital... Hillshade, Slope and Aspect from LiDAR DTMs --- title: "Hillshade, Slope and Aspect from LiDAR DTMs" description: "Deriving hillshade, slope, and aspect rasters... DSM Generation from LiDAR with PDAL --- title: "DSM Generation from LiDAR with PDAL" description: "Building a Digital Surface Model from LiDAR first returns... PMF Ground Classification in PDAL --- title: "PMF Ground Classification in PDAL" description: "Using the Progressive Morphological Filter (filters.pmf) in... SMRF Ground Classification in PDAL --- title: "SMRF Ground Classification in PDAL" description: "How the Simple Morphological Filter (filters.smrf)...