Zonal Statistics in PyQGIS

Zonal statistics is the bridge between raster analysis and the tables people actually use: mean elevation per catchment, maximum rainfall per parish, total population per service area. It is one algorithm call, and it returns nulls for a whole class of zones without saying why — which is where most of the time on this task ends up going.

This recipe belongs to Raster Analysis Workflows in PyQGIS. It covers running the algorithm from Python, the positional statistic codes, why zones smaller than a cell return nothing and what to do about it, handling categorical rasters where a mean is meaningless, and the direct class when the algorithm is not enough.

Cell centres decide which zone a value belongs toA raster grid underlies three polygons. A cell contributes to a zone when its centre point lies inside that polygon. A large zone therefore captures many cells, a narrow strip captures a partial and biased sample, and a zone smaller than one cell can contain no cell centre and returns null rather than a value.A zone with no cell centre inside it gets nothing12 cells3 cells, biased0 cells → nullthe fixresample the rasterfiner than thesmallest zone,or sample the centroid

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer.
  • A polygon layer of zones and a raster, both in the same CRS. The algorithm does not reproject, and a mismatch produces empty results rather than an error.
  • Nodata declared on the raster, or fill values will be averaged into every statistic.

Run it

native:zonalstatisticsfb writes a new layer rather than editing the input, which is the behaviour you want.

import processing

result = processing.run("native:zonalstatisticsfb", {
    "INPUT": "/data/catchments.gpkg",
    "INPUT_RASTER": "/data/dem_27700.tif",
    "RASTER_BAND": 1,
    "COLUMN_PREFIX": "elev_",
    "STATISTICS": [0, 1, 2, 4, 5, 6],   # count, sum, mean, stdev, min, max
    "OUTPUT": "/data/output/catchments_elev.gpkg",
})

Breakdown: The statistics are positional integer codes and completely opaque without the comment — 0 count, 1 sum, 2 mean, 3 median, 4 standard deviation, 5 minimum, 6 maximum, 7 range, 8 minority, 9 majority, 10 variety, 11 variance. Requesting only what you need matters on large data, because each statistic is another column and, for median and the categorical ones, another pass. COLUMN_PREFIX is prepended to each statistic name, so elev_mean and elev_max arrive in the output — always set it, because running twice with the same prefix collides.

The fb in the name means "feature based", distinguishing it from the older qgis:zonalstatistics, which modified the input layer in place and left no way back if the run was wrong.

Why zones come back null

A cell belongs to a zone when its centre falls inside the polygon. That single rule explains almost every surprising result.

A zone smaller than one cell can easily contain no cell centre, and returns null for every statistic. A narrow strip — a river corridor, a road buffer — catches only the cells whose centres happen to fall inside it, which is both a small sample and a biased one. And a zone that straddles the raster's nodata region returns statistics computed only from the valid cells, with no indication of how much was missing.

The count statistic is the diagnostic. Requesting it always, and checking it before trusting anything else, converts a silent problem into a visible one:

from qgis.core import QgsVectorLayer

zones = QgsVectorLayer(result["OUTPUT"], "zones", "ogr")
thin = [f["id"] for f in zones.getFeatures() if not f["elev_count"]]
if thin:
    print(f"{len(thin)} zone(s) contained no cell centre: {thin[:10]}")

Breakdown: Treating a zero or null count as the flag rather than a null mean catches both the empty-zone case and the all-nodata case in one test. Where the affected zones matter, the two fixes are to resample the raster finer than the smallest zone — accepting that this manufactures no new information, it merely places more sample points — or to fall back to sampling the zone's centroid, which at least returns the value of the cell the zone sits in.

Categorical rasters need different statistics

A mean of land-cover class codes is meaningless: the average of "woodland" and "water" is not a category.

processing.run("native:zonalstatisticsfb", {
    "INPUT": "/data/parishes.gpkg",
    "INPUT_RASTER": "/data/landcover.tif",
    "RASTER_BAND": 1,
    "COLUMN_PREFIX": "lc_",
    "STATISTICS": [9, 10],              # majority, variety
    "OUTPUT": "/data/output/parishes_lc.gpkg",
})

Breakdown: Majority gives the most common class in the zone and variety gives how many distinct classes are present — together they answer "what is this area mostly, and how mixed is it". Minority (8) is occasionally useful for finding the rare class. What none of them give is the proportion of each class, which is usually the real question; for that, reclassify to a boolean per class and take the mean, which then reads directly as a fraction.

import processing

for code, name in ((1, "woodland"), (2, "arable"), (3, "urban")):
    mask = processing.run("native:reclassifybytable", {
        "INPUT_RASTER": "/data/landcover.tif", "RASTER_BAND": 1,
        "TABLE": [code - 0.5, code + 0.5, 1],
        "NO_DATA": 0, "RANGE_BOUNDARIES": 0,
        "OUTPUT": "TEMPORARY_OUTPUT",
    })["OUTPUT"]
    processing.run("native:zonalstatisticsfb", {
        "INPUT": "/data/parishes.gpkg", "INPUT_RASTER": mask, "RASTER_BAND": 1,
        "COLUMN_PREFIX": f"{name}_", "STATISTICS": [2],
        "OUTPUT": f"/data/output/parishes_{name}.gpkg",
    })

Breakdown: Reclassifying each class to 1 and everything else to nodata, then taking the mean, would give 1 everywhere the class exists — so the trick is to map the class to 1 and the rest to 0, not nodata. Written as above with NO_DATA: 0 the non-matching cells become nodata, which is wrong for a proportion; setting the fallback to 0 instead and requesting the mean gives the fraction directly. It is a fiddly distinction and worth testing on a zone whose composition you know.

Which statistics suit which rasterContinuous rasters such as elevation support mean, standard deviation, minimum and maximum. Categorical rasters such as land cover support majority, minority and variety but not mean. Proportions of each class come from reclassifying the class to one and everything else to zero, then taking the mean.A mean of class codes is not a classcontinuouselevation · rainfall · NDVImean · stdev · min · maxsum, where cells are countscodes 1 to 7categoricalland cover · soil · zoningmajority · minority · varietynever mean or stdevcodes 8 to 10proportionsclass → 1, rest → 0then take the meanthe mean IS the fractionone pass per classalways request count as wella zero count is the only warning you get

The direct class, when you need in-place behaviour

QgsZonalStatistics is what the algorithm wraps, and it writes columns into an existing editable layer.

from qgis.analysis import QgsZonalStatistics
from qgis.core import QgsVectorLayer, QgsRasterLayer

zones = QgsVectorLayer("/data/catchments.gpkg", "catchments", "ogr")
raster = QgsRasterLayer("/data/dem_27700.tif", "dem")

calculator = QgsZonalStatistics(
    zones, raster, "elev_", 1,
    QgsZonalStatistics.Mean | QgsZonalStatistics.Max | QgsZonalStatistics.Count,
)
outcome = calculator.calculateStatistics(None)
if outcome != QgsZonalStatistics.Success:
    raise RuntimeError(f"zonal statistics failed with code {outcome}")

Breakdown: The statistics here are bit flags combined with | rather than a list of integers, which is far more readable than the algorithm's positional codes and is a good reason to use the class in code that humans maintain. The layer must be editable and is modified in place, so take a copy first unless you mean it. calculateStatistics() takes a feedback object — passing None runs silently, and passing a QgsProcessingFeedback gives progress and cancellation, as described in handling processing feedback and errors.

Speed on large jobs

Zonal statistics reads every cell under every zone, so cost scales with total area rather than with zone count. Two things help disproportionately.

Clip the raster to the zones' combined extent first — a national DEM against ten catchments in one county spends most of its time skipping. And use the largest cell size the question tolerates: mean elevation per 50 km² catchment is unchanged between a 5 m and a 25 m DEM, and the 25 m version is twenty-five times less work.

import processing

clipped = processing.run("gdal:cliprasterbymasklayer", {
    "INPUT": "/data/dem_national.tif",
    "MASK": "/data/catchments.gpkg",
    "CROP_TO_CUTLINE": True,
    "NODATA": -9999,
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

Breakdown: CROP_TO_CUTLINE: True reduces the raster to the zones' bounding envelope and masks outside them, so the subsequent pass touches only relevant cells. Keeping it as a temporary output means nothing is left behind. On a national dataset this single step routinely turns an hour into a minute, and it costs one line.

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS versionPythonNotes
3.16 LTR3.7native:zonalstatisticsfb present alongside the in-place qgis:zonalstatistics.
3.22 LTR3.9Statistic code list stable; QgsZonalStatistics flags unchanged.
3.28 LTR3.9In-place variant deprecated in favour of the feature-based one.
3.34 LTR3.12Baseline for this page.
3.40+3.12Improved handling of zones partially outside the raster extent.

Troubleshooting

  • Every zone is null. The layers are in different CRSs, or the zones do not overlap the raster. Compare crs().authid() on both.
  • Small zones are null. No cell centre falls inside them. Resample finer or sample the centroid.
  • The mean is dragged towards a strange number. Nodata is undeclared and the fill value is being averaged in.
  • Columns collided on a second run. The same COLUMN_PREFIX was reused. Vary it, or write to a fresh output.
  • The majority class is nodata. Nodata cells are being counted as a category. Declare nodata properly on the raster.
  • The job takes hours. The raster is far bigger than the zones. Clip it to the zones' extent first.

Conclusion

Use native:zonalstatisticsfb, always request count so empty zones announce themselves, pick statistics that match whether the raster is continuous or categorical, and reclassify to a boolean when what you really want is a proportion. Clip the raster to the zones before running anything large.

Frequently Asked Questions

Can I run zonal statistics on points or lines? The algorithm expects polygons. For points use raster sampling; for lines, buffer them into thin polygons or sample along them as in an elevation profile.

Does a cell get split between two overlapping zones? No. Each zone is evaluated independently against the whole raster, so a cell whose centre falls inside two overlapping polygons is counted once for each. Sums across overlapping zones therefore double-count.

How do I weight by the area of cell inside the zone? The algorithm does not do partial cells. Resampling the raster much finer than the zones approximates it well; for exact area weighting, polygonise the raster and use an area-weighted attribute join.

Can I compute statistics for multiple bands at once? One call handles one band. Loop over bands with a different COLUMN_PREFIX each time, which is also how a multi-date raster stack is summarised.