Terrain & Interpolation Analysis in PyQGIS

Terrain analysis is the part of geoprocessing where a wrong answer looks entirely convincing. A slope raster computed from a DEM in degrees of latitude produces gradients that are wrong by a factor of about a hundred thousand and still renders as a plausible-looking picture of hills. An interpolated surface built from thirty boreholes covers the whole study area with numbers that carry no warning about the twenty kilometres between samples.

This guide sits inside Spatial Data Processing & Automation and covers the two halves of the subject: deriving products from an elevation raster you already have, and building a raster surface from scattered point measurements. Both are short pieces of PyQGIS wrapped around decisions — units, cell size, search radius — that determine whether the output means anything.

Where terrain products come fromAn elevation raster feeds directly into slope, aspect, hillshade and contour derivation. Scattered point measurements must first be interpolated into a continuous surface, after which they feed the same derivative chain. Zonal statistics summarise any of these rasters back into vector features.Every product below depends on the same grid being rightelevation rasterDEM · DSM · DTMpoint samplesboreholes · gauges · spot heightsinterpolateIDW · TIN · krigingslope · aspecthillshadecontoursprofileszonalstatisticsback to vectorthe one thing that breaks all of ithorizontal units in degrees while vertical units are metres

The DEM is a grid of numbers, and the numbers have units

An elevation raster is a regular grid where each cell holds a height. Three properties of that grid determine everything downstream: the cell size, the vertical unit, and the horizontal unit of the CRS.

Slope is a ratio of a vertical change to a horizontal distance. If the raster is in EPSG:4326 the horizontal distance between cell centres is measured in degrees while the height is in metres, and every gradient computed from it is nonsense — around a hundred thousand times too shallow at mid-latitudes. QGIS does not warn about this, because a raster in a geographic CRS is perfectly legitimate for many other purposes.

from qgis.core import QgsRasterLayer

dem = QgsRasterLayer("/data/dem.tif", "dem")
crs = dem.crs()
print(crs.authid(), crs.isGeographic(), dem.rasterUnitsPerPixelX())

Breakdown: isGeographic() returning True is the flag that every terrain script should check before doing anything else. rasterUnitsPerPixelX() reports the cell size in the CRS's own units, so a value like 0.0002777 confirms degrees at a glance where 25.0 confirms metres. The fix is to reproject the DEM into a projected CRS before deriving anything, resampling bilinearly because elevation is continuous.

Cell size sets the ceiling on what any derivative can say. A slope raster from a 30 m DEM describes 30 m averages; it cannot tell you about a 5 m embankment, no matter how the result is styled. Resampling a coarse DEM to a fine grid produces smooth-looking output with no additional information in it, and is one of the more common ways an analysis oversells itself.

Slope, aspect, hillshade and the z-factor

The three standard derivatives are all one call each, through either the native algorithms or their GDAL equivalents.

import processing

processing.run("native:slope", {
    "INPUT": "/data/dem_27700.tif",
    "Z_FACTOR": 1.0,
    "OUTPUT": "/data/output/slope.tif",
})

Breakdown: Z_FACTOR scales the vertical values before the gradient is computed, and it exists for exactly one reason: mismatched units. A DEM in feet on a CRS in metres needs 0.3048; a DEM and CRS both in metres needs 1.0. It is not a style control, and using it to exaggerate relief silently corrupts every number the raster reports. Exaggeration belongs in the hillshade, which is a picture, not in the slope, which is a measurement.

Native slope returns degrees. gdal:slope offers a percent option through AS_PERCENT, which matters because engineering thresholds are usually written as percentages while ecological ones are usually degrees — and 30 percent and 30 degrees are very different slopes.

Hillshade is the one derivative that is purely cartographic, and it is the place where exaggeration is legitimate:

processing.run("native:hillshade", {
    "INPUT": "/data/dem_27700.tif",
    "Z_FACTOR": 2.0,
    "AZIMUTH": 315,
    "V_ANGLE": 45,
    "OUTPUT": "/data/output/hillshade.tif",
})

Breakdown: An azimuth of 315 — light from the north-west — is the cartographic convention, and departing from it triggers relief inversion, where readers see valleys as ridges. The vertical angle of 45 degrees is a reasonable default; lower angles lengthen shadows and pick out subtle features at the cost of blowing out steep ground. Combined with a thematic layer set to multiply blending, a hillshade turns a flat choropleth into a map that reads as terrain.

How slope and aspect are computed per cellA three by three neighbourhood is centred on each cell. Differences between the west and east columns give the gradient in the x direction; differences between the north and south rows give the gradient in the y direction. Slope is the magnitude of that gradient and aspect is its direction, which is why both are undefined on perfectly flat ground.Slope is a magnitude; aspect is the direction of the same vector10410811310110611299103109the 3×3 window, in metresdz/dx, dz/dyfrom opposite columns and rowsslope = |gradient|aspect = its bearingon flat groundthe gradient is zero,so aspect has nomeaningful value

Contours: turning a surface back into lines

Contours are the oldest terrain product and still the most readable on paper.

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

Breakdown: The interval should be chosen against the terrain and the output scale, not against a round number: ten metres on a flat floodplain produces three lines, and on a mountainside produces an unreadable mat. FIELD_NAME writes the elevation into an attribute, which is what makes index contours possible later — a rule-based renderer can then draw every fifth line heavier and label only those.

Contours inherit every artefact in the DEM. Stepping in a coarse or integer-valued raster shows up as terracing, and the honest fix is a mild smoothing of the DEM before contouring rather than smoothing the lines afterwards, which produces contours that no longer match the surface they claim to describe.

Interpolating a surface from points

When there is no raster, one has to be built from samples. QGIS ships inverse distance weighting and TIN interpolation natively, and both take a slightly awkward layer-specification string.

layer = QgsProject.instance().mapLayersByName("boreholes")[0]
spec = f"{layer.source()}::~::0::~::{layer.fields().indexOf('depth')}::~::0"

processing.run("qgis:idwinterpolation", {
    "INTERPOLATION_DATA": spec,
    "DISTANCE_COEFFICIENT": 2.0,
    "EXTENT": layer.extent(),
    "PIXEL_SIZE": 25,
    "OUTPUT": "/data/output/depth_idw.tif",
})

Breakdown: The ::~:: separator is a QGIS-internal encoding of source, geometry type, attribute index and use-z flag; building it by hand is unpleasant but it is the documented interface. DISTANCE_COEFFICIENT is the power in the weighting: 2 is the usual default, higher values make the surface hug each sample and produce visible bull's-eyes, lower values smooth towards the mean. PIXEL_SIZE should be no finer than about a quarter of the typical sample spacing — anything finer manufactures detail the data cannot support.

TIN interpolation is the other native option and behaves quite differently: it triangulates between samples and interpolates linearly inside each triangle, so the surface passes exactly through every measured point and stops abruptly at the convex hull. That makes it honest about extrapolation in a way IDW is not, and a poor choice when the underlying phenomenon is smooth. The IDW guide and the TIN guide work through the trade-off with real parameters.

Choosing an interpolation method

The methods available in and around QGIS differ in what they assume, and the assumption is the thing to match against your data rather than the name.

IDW assumes only that nearby samples are more relevant than distant ones. It never produces a value outside the range of the observations, which makes it safe and slightly dull, and it always produces some value, which makes it dangerous far from the samples. It is the right default for dense, evenly spread measurements — rainfall gauges across a region, soil samples on a grid.

TIN assumes the surface is piecewise planar between samples. It reproduces every observation exactly and refuses to extrapolate past the convex hull, which is the honest behaviour for spot heights and survey points where the measurements are trusted and the gaps are genuinely unknown. Its weakness is visible: triangular facets and creases along the triangle edges, which look wrong on anything meant to be smooth.

Kriging, available through the SAGA and GRASS providers rather than natively, assumes the spatial correlation structure can be estimated from the data and produces a variance surface alongside the estimate. That second output is the reason to reach for it — a map of how uncertain each cell is answers the question the other two methods leave implicit.

for power in (1.0, 2.0, 4.0):
    processing.run("qgis:idwinterpolation", {
        "INTERPOLATION_DATA": spec,
        "DISTANCE_COEFFICIENT": power,
        "EXTENT": layer.extent(),
        "PIXEL_SIZE": 25,
        "OUTPUT": f"/data/output/idw_p{int(power)}.tif",
    })

Breakdown: Running the same interpolation at three powers and comparing takes a minute and settles an argument that otherwise runs on intuition. At a power of 1 the surface is broad and smooth; at 4 each sample sits in its own small plateau with steep walls between. The right value is the one where features you can independently justify survive and features you cannot disappear — the same test that applies to a heatmap radius.

Validation deserves more than a look. Hold back a random tenth of the samples, interpolate from the rest, and sample the result at the held-back points: the mean absolute difference between predicted and observed is a number you can put in a report, and it usually deflates confidence in a satisfying way.

Getting numbers back out

A raster is not a deliverable for most audiences. Zonal statistics summarises it into the vector features people actually work with.

processing.run("native:zonalstatisticsfb", {
    "INPUT": "/data/catchments.gpkg",
    "INPUT_RASTER": "/data/output/slope.tif",
    "RASTER_BAND": 1,
    "COLUMN_PREFIX": "slope_",
    "STATISTICS": [2, 3, 6],          # mean, median, max
    "OUTPUT": "/data/output/catchments_slope.gpkg",
})

Breakdown: native:zonalstatisticsfb writes a new layer, unlike the older in-place qgis:zonalstatistics which modified the input and caused a good deal of grief. The STATISTICS codes are positional and worth writing with a comment, because [2, 3, 6] is unreadable six months later. Small polygons relative to the cell size are the usual source of nulls here: a zone smaller than one cell may contain no cell centre at all, and the result is null rather than the value of the cell it sits in.

Nodata, edges and the halo of wrong values

Every neighbourhood operation has an edge problem. A 3×3 window centred on a cell at the raster boundary has no neighbours on one side, and a window overlapping a nodata gap has no neighbours on that side either. What the algorithm does about it determines whether the output has a rim of nonsense.

QGIS's native derivatives set the output to nodata wherever the window is incomplete, which is the conservative and correct choice — one cell of nodata around every gap and along every edge. GDAL's versions offer COMPUTE_EDGES, which extrapolates rather than dropping the cell. That flag is convenient for a hillshade, where an approximate value at the border is better than a visible black line, and inadvisable for slope, where it invents measurements.

import processing

processing.run("gdal:hillshade", {
    "INPUT": "/data/dem_27700.tif",
    "BAND": 1,
    "Z_FACTOR": 2.0,
    "COMPUTE_EDGES": True,
    "OUTPUT": "/data/output/hillshade_edges.tif",
})

Breakdown: Turning edge computation on for the picture and leaving it off for the measurements is the arrangement that keeps both honest. If the DEM is a tile from a larger dataset, the better fix is to derive from a buffered extract and clip afterwards, so the edge artefacts fall outside the area you keep — see clipping a raster by a mask layer.

Gaps inside the DEM deserve a decision rather than a default. Filling them by interpolation before deriving anything gives a continuous product with quietly invented terrain in the holes; leaving them produces an honest raster with a nodata halo one cell wide around each. Which is right depends entirely on whether a downstream consumer will treat a value as a measurement, and it is worth recording the choice next to the code.

Working at scale

Terrain algorithms are memory-bound and read every pixel, so the difference between a workable script and an overnight job is usually structural rather than a matter of tuning.

Three habits do most of the work. Process by tile rather than by scene: a national DEM handled as a mosaic of tiles with a small overlap runs in bounded memory and parallelises across processes, and the overlap is what stops the tile seams appearing in the derivatives. Write intermediates as compressed GeoTIFFs with tiling enabled, since a slope raster of Float32 zeros compresses to almost nothing and the read cost dominates. And skip materialising anything that is only consumed once — a virtual raster describing the calculation avoids an entire read-write cycle.

processing.run("gdal:buildvirtualraster", {
    "INPUT": ["/data/tiles/n01.tif", "/data/tiles/n02.tif"],
    "RESOLUTION": 0,
    "SEPARATE": False,
    "OUTPUT": "/data/output/mosaic.vrt",
})

Breakdown: A VRT is a small XML file describing where the pixels live rather than a copy of them, so building one over a directory of tiles costs milliseconds and gives every downstream algorithm a single seamless input. RESOLUTION: 0 keeps the highest resolution among the inputs; the alternatives average or take the lowest, neither of which is usually what a DEM mosaic wants. For the batch loop that runs a derivative over each tile, see running an algorithm over a folder of files.

Recording what the surface claims

Every raster in this chapter is a claim about places nobody measured, and the difference between a defensible deliverable and a pretty picture is usually documentation rather than technique. Three things belong with the output.

The provenance: which DEM or sample set, which date, which CRS, which cell size. A GeoTIFF can carry this in its metadata, and QGIS reads and writes it through the layer's metadata object.

from qgis.core import QgsRasterLayer

surface = QgsRasterLayer("/data/output/rainfall_idw.tif", "rainfall")
metadata = surface.metadata()
metadata.setTitle("Annual rainfall, IDW from 84 gauges")
metadata.setAbstract(
    "Inverse distance weighting, power 2, 100 m cells, masked to 5 km "
    "of a gauge. Mean absolute error 41 mm on a 10 percent holdout."
)
metadata.setLanguage("en")
surface.setMetadata(metadata)
surface.saveDefaultMetadata()

Breakdown: saveDefaultMetadata() writes a sidecar .qmd next to the raster, so the description travels with the file rather than living only in the project. Putting the parameters and the error in the abstract is what makes the file self-describing six months later, when the script that made it has been edited twice. More on this in reading and writing layer metadata.

The extent of validity: the mask polygon, saved alongside the raster rather than discarded after clipping. And the error estimate from cross-validation, which turns "this is our rainfall surface" into "this is our rainfall surface, accurate to about 40 mm where gauges are within 5 km". The second sentence is the one that survives review.

Key takeaways

  • Check crs().isGeographic() before any terrain derivative; a DEM in degrees produces slopes that are wrong by orders of magnitude.
  • Z_FACTOR corrects a vertical-unit mismatch. It is not a style control, and exaggeration belongs only in a hillshade.
  • Cell size caps what a derivative can describe; resampling finer adds smoothness, not information.
  • Keep hillshade azimuth at 315 degrees unless you have a specific reason, or readers will see the relief inverted.
  • Contour interval is a cartographic decision — choose it against the terrain and output scale, not for roundness.
  • IDW smooths towards the mean and always produces a value; TIN passes through every sample and stops at the hull.
  • Interpolate no finer than about a quarter of the sample spacing.
  • Use native:zonalstatisticsfb rather than the in-place variant, and expect nulls where zones are smaller than a cell.

Frequently Asked Questions

Why is my slope raster almost entirely zero? The DEM is in a geographic CRS, so the horizontal distances are degrees and the computed gradients are minuscule. Reproject to a projected CRS with bilinear resampling and recompute.

Should I use the native or the GDAL terrain algorithms? They agree closely. GDAL's versions expose more options — percent slope, edge handling, the Zevenbergen-Thorne algorithm — and are generally faster on large files. Native ones are simpler and integrate more smoothly with the Processing framework.

What cell size should I interpolate to? No finer than roughly a quarter of the average distance between samples. Beyond that the output is smoothly rendered guesswork, and it invites readers to zoom in to a level of detail the data never had.

Can I interpolate without an extent? The algorithms need one. Using the point layer's own extent clips the surface to the data, which is usually right; extending it beyond the samples means extrapolating, and IDW extrapolates to the global mean without saying so.

How do I mask out areas with no nearby samples? Buffer the sample points by a defensible search distance, dissolve, and clip the interpolated raster to that polygon. It converts an implicit claim into an explicit one and takes two algorithm calls.

Does a hillshade need to be the same resolution as the DEM? It is produced at the DEM's resolution by default and that is normally correct. Rendering a coarse hillshade under fine data looks blurry; the answer is a better DEM rather than resampling the hillshade.