Create Contours from a DEM in PyQGIS

Contours are the most readable terrain product ever invented and the easiest to get wrong. The algorithm call is one line; the decisions around it — what interval, how to handle a DEM that steps in whole metres, which lines to label — are what separate a map somebody can navigate with from a mat of tangled spaghetti.

This recipe belongs to Terrain & Interpolation Analysis in PyQGIS. It covers producing lines and filled bands, choosing an interval against the terrain rather than for roundness, dealing with terracing in coarse or integer DEMs, and setting up index contours with labels that follow the curve.

Prerequisites

QGIS 3.34 LTR (bundled Python 3.12) or newer.
A DEM with nodata declared. Undeclared fill values produce contour rings around every gap at absurd elevations.
A projected CRS is not strictly required for contouring — the vertical values are read directly — but it is required for anything you do with the result, and for the smoothing step below.

Generate the lines

gdal:contour is the workhorse and takes a fixed interval.

import processing

processing.run("gdal:contour", {
    "INPUT": "/data/dem_27700.tif",
    "BAND": 1,
    "INTERVAL": 10.0,
    "FIELD_NAME": "elev",
    "BASE": 0.0,
    "CREATE_3D": False,
    "IGNORE_NODATA": False,
    "NODATA": None,
    "OFFSET": 0.0,
    "OUTPUT": "/data/output/contours_10m.gpkg",
})

Breakdown: FIELD_NAME writes the elevation of each line into an attribute, and skipping it produces a layer of unlabelled, unstyleable lines — always set it. BASE anchors the sequence, so a base of 0 with an interval of 10 yields lines at 0, 10, 20; a base of 5 yields 5, 15, 25, which is occasionally what a local convention requires. IGNORE_NODATA set to False is the safe default: with it True, GDAL treats gaps as interpolable and draws contours across them.

CREATE_3D writes the elevation into the geometry's Z coordinate as well as the attribute. It is worth turning on when the contours will feed a TIN interpolation or a 3D view, and worth leaving off otherwise, because a Z-aware layer is fussier to edit and some downstream algorithms drop the dimension silently.

Choosing the interval

There is no correct interval, only intervals that suit the terrain and the output scale.

The usable rule is a limit on line density: contours closer than about half a millimetre on the printed page merge into a solid tone. At 1:25 000, half a millimetre is 12.5 metres on the ground, so on a slope of 30 degrees — about 0.58 metres of rise per metre of run — the closest acceptable interval is around 7 metres. Round to 10 and the steepest ground stays readable.

def suggest_interval(map_scale, max_slope_degrees, min_mm=0.5):
    import math
    ground_mm = map_scale * min_mm / 1000.0
    return ground_mm * math.tan(math.radians(max_slope_degrees))

print(suggest_interval(25000, 30))     # ≈ 7.2 metres

Breakdown: The function converts the minimum printable separation into a ground distance, then multiplies by the gradient to get the vertical interval that produces exactly that separation on the steepest ground. It is a floor rather than a recommendation — round up to a conventional number, and remember that the steepest slope in the DEM is often a data artefact rather than real terrain, so use a high percentile from zonal statistics rather than the maximum.

Terracing, and how to remove it honestly

A DEM stored as integers, or resampled from a coarser source, has flat steps. Contouring it produces contours that hug the step edges as angular, blocky lines that follow the pixel grid rather than the terrain.

smoothed = processing.run("native:cellstatistics", {
    "INPUT": ["/data/dem_27700.tif"],
    "STATISTIC": 2,
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

filtered = processing.run("grass7:r.neighbors", {
    "input": "/data/dem_27700.tif",
    "size": 3,
    "method": 0,                       # average
    "output": "TEMPORARY_OUTPUT",
})["output"]

processing.run("gdal:contour", {
    "INPUT": filtered, "BAND": 1, "INTERVAL": 10.0,
    "FIELD_NAME": "elev", "OUTPUT": "/data/output/contours_smooth.gpkg",
})

Breakdown: Smoothing the raster before contouring is the honest fix, because the resulting lines still describe a real surface — just a slightly generalised one. Smoothing the lines afterwards with native:smoothgeometry looks similar and is not equivalent: the smoothed line no longer sits at the elevation its attribute claims, so a point taken from it and sampled against the DEM disagrees. Use line smoothing only for purely decorative output where nobody will measure anything.

A 3×3 average is a gentle filter. If terracing survives it, the DEM's vertical resolution is genuinely coarser than the contour interval you have chosen, and the correct response is a larger interval rather than heavier smoothing.

Cleaning up the output

Raw contour output usually needs two passes before it is presentable, and both are cheap.

Short fragments accumulate wherever the surface just grazes a level — a handful of vertices describing a bump two metres across, which at map scale is a dot. Filtering by length removes them without touching anything real.

import processing

processing.run("native:extractbyexpression", {
    "INPUT": "/data/output/contours_10m.gpkg",
    "EXPRESSION": "$length > 40",
    "OUTPUT": "/data/output/contours_clean.gpkg",
})

Breakdown: The threshold is a map-scale judgement rather than a terrain one: at 1:25 000, forty metres is about 1.6 mm on the page, which is roughly the shortest line worth drawing. Using $length rather than a geometry function keeps the expression readable and lets it be tuned in the interface before being written into the script.

The second pass is ordering. Contours come out of GDAL in the order the algorithm found them, which means the drawing order is arbitrary and index contours can end up underneath ordinary ones where they cross. Sorting by elevation on export fixes it, and costs nothing.

processing.run("native:orderbyexpression", {
    "INPUT": "/data/output/contours_clean.gpkg",
    "EXPRESSION": '"elev"',
    "ASCENDING": True,
    "NULLS_FIRST": False,
    "OUTPUT": "/data/output/contours_final.gpkg",
})

Breakdown: Ordering a layer by an attribute fixes draw order for renderers that respect feature order, which includes the ordinary single-symbol and rule-based renderers. It is not a substitute for a symbol-level rendering pass when two classes genuinely must be drawn in separate sweeps, but for contours — where all the lines are the same kind of thing at different heights — it is exactly the right tool.

Index contours and labels that follow the line

Every fifth contour drawn heavier and labelled is the convention that makes a contour map readable at a glance. With the elevation in an attribute it is a rule-based renderer plus a label filter.

from qgis.core import (
    QgsVectorLayer, QgsProject, QgsPalLayerSettings, QgsProperty,
    QgsLabeling, QgsVectorLayerSimpleLabeling,
)

contours = QgsVectorLayer("/data/output/contours_10m.gpkg", "contours", "ogr")
QgsProject.instance().addMapLayer(contours)

settings = QgsPalLayerSettings()
settings.fieldName = "elev"
settings.placement = QgsPalLayerSettings.Line
settings.lineSettings().setPlacementFlags(
    QgsLabeling.OnLine | QgsLabeling.MapOrientation
)
settings.isExpression = False
settings.obstacle = False

filter_settings = settings
filter_settings.dataDefinedProperties().setProperty(
    QgsPalLayerSettings.Show,
    QgsProperty.fromExpression('"elev" % 50 = 0'),
)

contours.setLabeling(QgsVectorLayerSimpleLabeling(filter_settings))
contours.setLabelsEnabled(True)
contours.triggerRepaint()

Breakdown: The Show data-defined property is the clean way to label only some features — it evaluates per feature and suppresses the rest, without needing a separate layer or a subset string that would also hide the lines. placement = Line with OnLine puts the number in a gap in the contour rather than beside it, which is the cartographic convention and depends on the renderer being told to break the line under the label. obstacle = False stops the contours themselves being treated as things labels must avoid; leave it True and most labels are dropped. More on this in curved labels along lines.

Filled bands instead of lines

Sometimes the deliverable is coloured elevation zones rather than lines.

processing.run("gdal:contour_polygon", {
    "INPUT": "/data/dem_27700.tif",
    "BAND": 1,
    "INTERVAL": 50.0,
    "FIELD_NAME_MIN": "elev_min",
    "FIELD_NAME_MAX": "elev_max",
    "OUTPUT": "/data/output/bands_50m.gpkg",
})

Breakdown: The polygon variant writes both bounds of each band into attributes, so a graduated renderer on elev_min produces a hypsometric tint with no manual classification. Bands are considerably heavier than lines — a 50 m interval over a mountainous DEM easily produces tens of thousands of polygons — so it is worth simplifying the geometries before publishing to a web service.

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS version	Python	Notes
3.16 LTR	3.7	`gdal:contour` present; `gdal:contour_polygon` requires GDAL 2.4+.
3.22 LTR	3.9	Label `Show` data-defined property stable for per-feature filtering.
3.28 LTR	3.9	`native:smoothgeometry` available if decorative line smoothing is wanted.
3.34 LTR	3.12	Baseline for this page.
3.40+	3.12	Contour output honours the input's coordinate precision more closely.

Troubleshooting

Contour rings around every gap. Nodata is undeclared, so -9999 is being contoured as terrain.
The output has no elevation attribute. FIELD_NAME was omitted. Re-run; there is no way to recover it afterwards.
Lines are blocky and follow the pixel grid. The DEM is integer-valued or resampled. Filter the raster, not the lines.
No labels appear. setLabelsEnabled(True) was not called, or the contours are acting as obstacles to their own labels.
The map is a mat of lines. The interval is too fine for the output scale. Compute the floor from the steepest ground and round up.
Smoothed lines disagree with the DEM. They were smoothed after contouring. Redo the smoothing on the raster.

Conclusion

Set FIELD_NAME, leave IGNORE_NODATA off, and pick the interval from the steepest ground and the output scale rather than from habit. Fix terracing in the raster before contouring, use a per-feature Show expression for index contours, and reach for the polygon variant only when the deliverable is genuinely a tinted band map.

Frequently Asked Questions

Can I generate contours at specific non-uniform levels? Yes — gdal:contour accepts a FIXED_LEVELS string of comma-separated values instead of an interval, which is how flood-level or bathymetric contours at conventional depths are produced.

Why are my contours in the wrong place after reprojecting? Contours are derived at the DEM's grid, so reprojecting the lines afterwards resamples geometry rather than re-deriving it. Reproject the DEM and contour the reprojected raster instead.

How do I get closed contours for a bounded area? Contours are open where they meet the raster edge. Clip the DEM to a slightly larger extent than the study area, contour, then clip the lines back — that way each line is cut cleanly rather than ending at an arbitrary tile boundary.

Is there a native contour algorithm? The GDAL one is the standard route and is what the interface uses. There is no separately maintained native equivalent, and there is no need for one.

Terrain & Interpolation Analysis in PyQGIS — the guide this recipe belongs to
Generate Slope, Aspect and Hillshade in PyQGIS
Build a TIN Interpolation in PyQGIS
Curved Labels Along Lines in PyQGIS
Simplify a Geometry in PyQGIS