Short answer

To turn a digital elevation model into contour lines with GDAL, run gdal_contour -a elev -i 10 dem.tif contours.gpkg. The -i flag sets the contour interval in the DEM's vertical unit (usually metres), and -a elev writes each line's elevation into an attribute called elev. The quality of the result is decided before you run the command: by the DEM resolution, by whether you smoothed it, and by choosing an interval that fits both the terrain relief and the map scale. Raw lidar DEMs without smoothing produce "spaghetti" contours that no cartographer will accept.

How gdal_contour works

gdal_contour walks the raster and traces the iso-lines where the surface crosses each elevation level. With -i, levels are evenly spaced starting from an offset you can set with -off. So -i 10 -off 0 produces 0, 10, 20, 30… If you only want specific levels — say bathymetric breaks at 5, 10, 20, 50, 100 m — list them explicitly with repeated -fl (fixed-level) flags instead of -i:

gdal_contour -a depth -fl 5 -fl 10 -fl 20 -fl 50 -fl 100 \
  bathy.tif breaks.gpkg

By default it produces line geometry. Since GDAL 2.4 you can produce filled bands (polygons between levels) with -p, writing minimum and maximum elevation attributes via -amin and -amax — useful for hypsometric tinting:

gdal_contour -p -amin elev_min -amax elev_max -i 50 \
  dem.tif bands.gpkg

Prefer GeoPackage (.gpkg) output over Shapefile. Shapefiles cap field names at 10 characters, cannot exceed 2 GB, and store no real CRS metadata; a dense contour set on a 1 m DEM blows past those limits easily.

Choosing the contour interval

The interval is a cartographic decision, not a default. Two factors drive it:

  • Map scale. A practical heuristic for topographic maps is interval ≈ scale-denominator / 2000 metres. That gives ~10 m at 1:25,000 and ~20 m at 1:50,000 — values you will recognise from national mapping series. At 1:10,000, 5 m intervals are typical.
  • Relief. In flat terrain (floodplains, plateaus) a 10 m interval can leave whole sheets with no lines; drop to 1–2 m. In high-relief alpine terrain, a 5 m interval produces an unreadable wall of lines; widen to 20–25 m and rely on index contours for orientation.

Add index contours — heavier, labelled lines every fifth interval — to make the map readable. gdal_contour does not flag them, so derive the flag afterward. In a GeoPackage you can add a field and update it with OGR SQL:

ogr2ogr -f GPKG contours_idx.gpkg contours.gpkg \
  -sql "SELECT *, CASE WHEN CAST(elev AS INTEGER) % 50 = 0 \
        THEN 1 ELSE 0 END AS is_index FROM contour"

Here index contours fall every 50 m; symbolise is_index = 1 with a thicker stroke and labels.

A complete worked workflow

Take a 1 m lidar DEM of a study catchment, target a 1:5,000 engineering plan with 1 m contours and 5 m index lines.

  1. Inspect the DEM. Confirm units and NoData:
gdalinfo -stats catchment_dem.tif

Verify the vertical unit is metres and the NoData value is set; contours traced across an unset NoData value can spike or wrap around voids.

  1. Smooth lightly. A 1 m lidar surface carries roof edges, vehicles, and noise that become jagged micro-contours. Apply a gentle low-pass filter before contouring. With GDAL's nearblack/gdal_calc toolset a simple option is a focal mean; many teams use a small Gaussian or the QGIS r.neighbors / Smooth algorithms. The goal is to remove sub-metre noise without flattening real breaks of slope.

  2. Generate contours into GeoPackage with the elevation attribute and edge handling:

gdal_contour -a elev -i 1 -off 0 \
  -f GPKG smoothed_dem.tif contours.gpkg
  1. Drop tiny loops. Lidar noise produces dozens of closed contours a few cells across. Filter them out by length:
ogr2ogr -f GPKG contours_clean.gpkg contours.gpkg \
  -dialect SQLITE \
  -sql "SELECT * FROM contour WHERE ST_Length(geom) > 15"

The 15 m threshold is illustrative; tune it to your DEM resolution and noise level so you remove artifacts without deleting real small knolls.

  1. Smooth the line geometry for cartography. gdal_contour outputs vertex-per-cell-crossing lines that look angular at large scale. Apply a line-smoothing step (QGIS Smooth with a few iterations, or a Chaikin/spline simplification) — applied to the lines, not the DEM, this is purely visual and must be documented so it is never mistaken for added accuracy.

  2. Flag and label index contours as shown above (elev % 5 = 0), then style: thin grey for intermediate lines, heavier brown for index lines with elevation labels following the slope.

Common pitfalls and why they happen

  • Spaghetti contours. Caused by contouring a noisy high-resolution DEM directly. Each tiny elevation wiggle crosses a level and spawns a line. Smooth the DEM first; do not just simplify the lines, because over-simplification then shifts contours away from the true surface.
  • Contours that "climb" through buildings. Lidar digital surface models (DSM) include rooftops and trees. For terrain contours you need a digital terrain model (DTM); contouring a DSM produces lines around every building.
  • Wrong values everywhere. The DEM's vertical unit was feet, not metres, but you assumed metres — your "10 m" interval is really ~3 m of relief. Always check units in gdalinfo.
  • Lines stop at tile edges. Contouring tiles separately leaves gaps at seams. Build a single mosaic or VRT and contour the whole extent at once.
  • Self-intersecting or invalid geometry after smoothing. Aggressive spline smoothing can push lines through each other. Re-check with ST_IsValid and reduce smoothing strength.

Validation

Confirm contour count and value range:

ogrinfo -so contours_clean.gpkg contour
ogrinfo -sql "SELECT MIN(elev), MAX(elev), COUNT(*) FROM contour" \
  contours_clean.gpkg

The elevation range must sit inside the DEM's min/max from gdalinfo. Overlay contours on a hillshade of the same DEM: lines should run perpendicular to the steepest descent, bunch tightly on steep slopes, spread on flats, and form closed loops around peaks and depressions. Where two adjacent contours touch or cross, the DEM or the smoothing is at fault — contours of different elevations can never intersect.

Bathyl perspective

Contours are an interpretation of a surface, and the smoothing choices that make them legible can quietly imply precision the DEM does not have. We keep the unsmoothed source contours archived next to the cartographic version and record the interval, smoothing method, and any line filtering, so the map reads cleanly while the underlying accuracy stays auditable.

Related reading

Sources