Generating Contours from a DTM with gdal_contour

TL;DR: gdal_contour -a elev -i 0.5 -f GPKG -nln contours dtm.tif contours.gpkg writes 0.5 m contour lines with an elev attribute. Add -p -amin elev_min -amax elev_max for filled elevation-band polygons, -fl 100 105 110 for fixed levels, and point it at a .vrt mosaic so lines are continuous across tiles. gdal.ContourGenerateEx does the same from Python.

# Context and Motivation

This guide is part of Contour Generation from LiDAR DTMs. gdal_contour ships with every GDAL installation, needs no code, and handles the geometry correctly — marching squares over the raster, lines joined into continuous features, NoData respected. It is the right first tool for contours, and the right production tool when contours are one step in a shell-driven batch. The Python API exposes exactly the same algorithm when contours are part of a larger Python workflow.

The command’s options map onto three decisions: which elevations to draw, what to store with each feature, and whether you want lines or filled bands.

# Prerequisites and Assumptions

• GDAL 3.x command-line tools and, for the Python route, the osgeo bindings.
• A void-filled, ideally smoothed DTM; see smoothing a LiDAR DTM before contouring.
• A NoData value set on the raster, so contours are not drawn around empty areas.

# Step-by-Step Implementation

# Step 1 — Mosaic tiles into a VRT

gdalbuildvrt mosaic.vrt dtm_tiles/*.tif creates a virtual raster with no copying; contouring it gives continuous lines across tile boundaries.

# Step 2 — Generate lines

-i sets the interval, -off the base offset, -a the elevation attribute name. -nln names the output layer.

# Step 3 — Or generate bands

-p switches to polygons; -amin and -amax name the attributes holding each band’s lower and upper bound.

# Step 4 — Use fixed levels when needed

-fl takes an explicit list of elevations — useful for flood levels, design heights or irregular intervals.

# Step 5 — Check the output

Count features, list distinct elevations, and look for contours along NoData edges.

# Complete Working Example

Shell, for a batch:

#!/usr/bin/env bash
set -euo pipefail
gdalbuildvrt -overwrite mosaic/dtm.vrt dtm_tiles/*.tif

# 0.5 m lines with an elevation attribute
gdal_contour -a elev -i 0.5 -off 0 -f GPKG -nln contours \
  mosaic/dtm.vrt out/contours_050.gpkg

# 1 m elevation bands for hypsometric tinting
gdal_contour -p -amin elev_min -amax elev_max -i 1.0 -f GPKG -nln bands \
  mosaic/dtm.vrt out/bands_100.gpkg

# design flood levels only
gdal_contour -a elev -fl 212.4 213.0 213.8 -f GPKG -nln flood_levels \
  mosaic/dtm.vrt out/flood_levels.gpkg

ogrinfo -so out/contours_050.gpkg contours | grep "Feature Count"

Python, when contours are part of a larger script:

"""Contour lines and bands from a DTM VRT with GDAL's Python API."""
from osgeo import gdal, ogr

gdal.UseExceptions()


def contour(src: str, dst: str, layer_name: str, interval: float, polygons: bool = False) -> int:
    ds = gdal.Open(src)
    band = ds.GetRasterBand(1)
    out = ogr.GetDriverByName("GPKG").CreateDataSource(dst)
    geom = ogr.wkbMultiPolygon if polygons else ogr.wkbLineString
    lyr = out.CreateLayer(layer_name, ds.GetSpatialRef(), geom)
    lyr.CreateField(ogr.FieldDefn("id", ogr.OFTInteger))
    if polygons:
        lyr.CreateField(ogr.FieldDefn("elev_min", ogr.OFTReal))
        lyr.CreateField(ogr.FieldDefn("elev_max", ogr.OFTReal))
        opts = [f"LEVEL_INTERVAL={interval}", "ID_FIELD=0", "ELEV_FIELD_MIN=1",
                "ELEV_FIELD_MAX=2", "POLYGONIZE=YES"]
    else:
        lyr.CreateField(ogr.FieldDefn("elev", ogr.OFTReal))
        opts = [f"LEVEL_INTERVAL={interval}", "ID_FIELD=0", "ELEV_FIELD=1"]
    nd = band.GetNoDataValue()
    if nd is not None:
        opts.append(f"NODATA={nd}")
    gdal.ContourGenerateEx(band, lyr, options=opts)
    n = lyr.GetFeatureCount()
    out = None
    return n


if __name__ == "__main__":
    import glob
    gdal.BuildVRT("mosaic/dtm.vrt", sorted(glob.glob("dtm_tiles/*.tif")))
    print("lines:", contour("mosaic/dtm.vrt", "out/contours_050.gpkg", "contours", 0.5))
    print("bands:", contour("mosaic/dtm.vrt", "out/bands_100.gpkg", "bands", 1.0, polygons=True))

# Using Elevation Bands

Polygon output is less familiar than lines but often more useful for analysis. Each band is an area, so questions that are awkward with lines become one-liners: how much of a parcel lies below a flood elevation, what share of a catchment is above 300 m, where does a site need cut or fill relative to a design level. With GeoPandas, intersect the bands with parcels and sum areas grouped by elev_min.

Bands also make hypsometric tinting trivial in any GIS: style by elev_min with a colour ramp and the landscape’s shape reads at a glance, without a separate raster. For that purpose bands can be coarser than the delivered contour lines — 5 or 10 m bands under 0.5 m lines is a common pairing.

Two cautions apply. Bands inherit every artefact of the surface, so smooth the DTM first; micro-relief creates thousands of tiny island polygons. And band boundaries are only as accurate as the DTM, so the same interval rule applies as for lines: do not slice more finely than the vertical accuracy supports. For flood mapping, a fixed-level band at the design elevation (-fl with -p) is usually what is really wanted.

# Key Parameter Table

# Verification

• Distinct elevations are exactly the interval series: ogrinfo -sql "SELECT DISTINCT elev FROM contours" out/contours_050.gpkg.
• No lines along NoData. Overlay contours on the DTM’s NoData mask; lines should not trace its boundary.
• Continuity. Pick a contour crossing a former tile edge and confirm it is one feature.

# Gotchas and Edge Cases

NoData not set. Without a NoData value, GDAL treats fill values such as −9999 as elevations and draws a dense ring of contours around every void. Set NoData on the DTM or pass -snodata.

Floating-point levels. An interval of 0.1 produces levels like 212.30000000000001. Round the attribute when cleaning, or use -fl with exact values.

Huge polygon features. Polygon mode on a large, detailed DTM creates bands with millions of vertices. Smooth first, and consider coarser intervals for bands than for lines.

Units. gdal_contour uses the raster’s values as they are. A DTM in feet produces contours in feet; make sure the interval is in the same unit.

# Frequently Asked Questions

How do I create contour lines from a DTM with GDAL?

Run gdal_contour with an attribute name, an interval and an output format, for example gdal_contour -a elev -i 0.5 -f GPKG dtm.tif contours.gpkg. From Python, gdal.ContourGenerateEx does the same.

How do I get filled elevation bands instead of lines?

Add the -p flag with -amin and -amax attribute names. Each output polygon covers the area between two consecutive levels and carries both bounds.

How do I contour many DTM tiles at once?

Build a VRT mosaic with gdalbuildvrt and contour the VRT. The lines are continuous across tile boundaries, which they would not be if each tile were contoured separately.

Can I choose specific elevations?

Yes. The -fl option takes a list of fixed levels, useful for flood elevations or design heights.

# Related

• Contour Generation from LiDAR DTMs — the full workflow
• Smoothing a LiDAR DTM Before Contouring — preparing the surface
• Exporting Contours to GeoPackage — cleaning and delivery
• Building a Seamless DTM Mosaic from Tiles — the mosaic contours come from
• Generating Slope and Aspect Rasters with gdaldem — the other GDAL terrain tool