Calculate Raster Statistics in PyQGIS

Summarizing the values in a raster — its minimum, maximum, mean, and standard deviation — is the starting point for choosing a color ramp, normalizing data, or quality-checking a DEM. PyQGIS exposes this directly through layer.dataProvider().bandStatistics(...), which returns a QgsRasterBandStatistics object, and through the native:zonalstatisticsfb algorithm for per-polygon summaries. Both are essential tools in Raster Analysis Workflows.

This page covers single-band statistics, building a histogram of value frequencies, and computing zonal statistics that summarize a raster within each polygon of a vector layer.

Prerequisites

• QGIS 3.34 LTR (bundled Python 3.12) with Processing available.
• A single- or multi-band raster (a GeoTIFF DEM, NDVI grid, or similar).
• For zonal statistics, a polygon vector layer overlapping the raster.
• The QGIS Python Console (Plugins > Python Console).

Band indices in the data-provider API are 1-based: band 1 is the first band, not band 0.

Read Band Statistics

Load the raster and ask its data provider for statistics on a band. The bandStatistics method accepts a flag for which statistics to compute and returns an object with the results.

from qgis.core import QgsRasterLayer, QgsRasterBandStats

raster = QgsRasterLayer("/data/dem.tif", "dem")
if not raster.isValid():
    raise RuntimeError("Raster failed to load")

provider = raster.dataProvider()
stats = provider.bandStatistics(1, QgsRasterBandStats.All)

print(f"Min:    {stats.minimumValue:.3f}")
print(f"Max:    {stats.maximumValue:.3f}")
print(f"Mean:   {stats.mean:.3f}")
print(f"StdDev: {stats.stdDev:.3f}")
print(f"Range:  {stats.range:.3f}")
print(f"Sum:    {stats.sum:.3f}")

Breakdown: bandStatistics(1, QgsRasterBandStats.All) computes the full set for band 1. QgsRasterBandStats.All forces a complete scan; passing a narrower flag such as QgsRasterBandStats.Min | QgsRasterBandStats.Max is faster when you only need extremes. The returned object exposes minimumValue, maximumValue, mean, stdDev, range, and sum as plain attributes — no method calls needed.

For large rasters where an exact scan is slow, request statistics over a sample extent and resolution to get a fast approximation:

extent = raster.extent()
sample_size = 250000
approx = provider.bandStatistics(1, QgsRasterBandStats.All, extent, sample_size)
print(f"Approx mean over {sample_size} px: {approx.mean:.3f}")

Breakdown: Supplying an extent and a sampleSize caps how many pixels are read, trading a little accuracy for speed on multi-gigabyte rasters. With no sample size the provider reads every pixel, which is exact but slow.

The min and max you get here are what you would feed into a color ramp — see Apply a Color Ramp to a Raster in PyQGIS for styling with these values.

Loop Over Multiple Bands

For multi-band imagery (e.g. RGB or multispectral), iterate the band count and collect per-band statistics:

from qgis.core import QgsRasterLayer, QgsRasterBandStats

raster = QgsRasterLayer("/data/landsat.tif", "scene")
provider = raster.dataProvider()

for band in range(1, raster.bandCount() + 1):
    s = provider.bandStatistics(band, QgsRasterBandStats.All)
    name = raster.bandName(band)
    print(f"{name}: mean={s.mean:.2f} min={s.minimumValue:.2f} max={s.maximumValue:.2f}")

Breakdown: bandCount() gives the number of bands and range(1, bandCount() + 1) produces 1-based indices. bandName(band) returns a human-readable label. This pattern is the basis for per-band normalization or building a summary report across a scene.

Build a Histogram

A histogram shows how many pixels fall into each value bin — useful for spotting bimodal distributions, saturation, or nodata spikes. The provider's histogram method returns the bin counts:

from qgis.core import QgsRasterLayer, QgsRasterBandStats

raster = QgsRasterLayer("/data/dem.tif", "dem")
provider = raster.dataProvider()
stats = provider.bandStatistics(1, QgsRasterBandStats.All)

bin_count = 20
hist = provider.histogram(
    1,                       # band
    bin_count,               # number of bins
    stats.minimumValue,      # range minimum
    stats.maximumValue,      # range maximum
)

width = (stats.maximumValue - stats.minimumValue) / bin_count
for i, count in enumerate(hist.histogramVector):
    lo = stats.minimumValue + i * width
    print(f"[{lo:8.2f}, {lo + width:8.2f}): {count}")

Breakdown: histogram(band, binCount, min, max) returns a QgsRasterHistogram; its histogramVector is a list of pixel counts per bin. Bounding the range with the actual min/max from bandStatistics keeps the bins meaningful. The loop reconstructs each bin's interval from the bin width so the output is readable.

Compute Zonal Statistics

Zonal statistics summarize raster values inside each polygon of a vector layer — average elevation per catchment, mean NDVI per field. Use native:zonalstatisticsfb, which writes the summary columns into a new vector layer (the "fb" variant returns a fresh feature-based output rather than editing in place):

import processing

result = processing.run("native:zonalstatisticsfb", {
    "INPUT": "/data/parcels.gpkg",
    "RASTER": "/data/dem.tif",
    "RASTER_BAND": 1,
    "COLUMN_PREFIX": "elev_",
    "STATISTICS": [0, 1, 2],   # 0=Count, 1=Sum, 2=Mean
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

for feat in result.getFeatures():
    print(feat["id"], feat["elev_mean"])

Breakdown: STATISTICS is a list of integer codes — 0 count, 1 sum, 2 mean, 3 median, 4 stddev, 5 min, 6 max, and so on. Each requested statistic becomes a field named COLUMN_PREFIX plus the statistic (e.g. elev_mean). The RASTER_BAND is 1-based, matching the provider API. Because the input and raster must overlap in the same CRS, reproject first if they differ — see Coordinate Reference Systems.

To pull the global mean directly into your zonal calc, you can pre-read it with bandStatistics and compare per-zone means against it for anomaly detection. Once you have the raster trimmed to a region with Clip a Raster by a Mask Layer in PyQGIS, the same statistics calls report on just that area.

Exclude Nodata from the Result

Statistics are only trustworthy if nodata pixels are excluded. The data-provider statistics respect the band's declared nodata value automatically, but if a raster carries an undeclared fill (a stray -9999 or 0), you must register it first so it does not skew the mean:

from qgis.core import QgsRasterLayer, QgsRasterBandStats, QgsRasterRange

raster = QgsRasterLayer("/data/dem.tif", "dem")
provider = raster.dataProvider()

# Declare the fill value as nodata for band 1, then recompute
provider.setUserNoDataValue(1, [QgsRasterRange(-9999, -9999)])
provider.setUseSourceNoDataValue(1, True)

stats = provider.bandStatistics(1, QgsRasterBandStats.All)
print(f"Mean excluding fill: {stats.mean:.3f}")

Breakdown: setUserNoDataValue registers -9999 as nodata for band 1 using a QgsRasterRange, and setUseSourceNoDataValue ensures any nodata already in the file is honored too. After this, bandStatistics ignores those pixels, so the mean and standard deviation reflect real data only. Import QgsRasterRange from qgis.core alongside the other classes.

Write Statistics to a Report

For documentation or QA it is useful to dump the figures to a small text or CSV report rather than just printing them. Reuse the per-band loop and write rows:

import csv
from qgis.core import QgsRasterLayer, QgsRasterBandStats

raster = QgsRasterLayer("/data/landsat.tif", "scene")
provider = raster.dataProvider()

with open("/data/output/band_stats.csv", "w", newline="") as fh:
    writer = csv.writer(fh)
    writer.writerow(["band", "name", "min", "max", "mean", "stddev"])
    for band in range(1, raster.bandCount() + 1):
        s = provider.bandStatistics(band, QgsRasterBandStats.All)
        writer.writerow([
            band, raster.bandName(band),
            f"{s.minimumValue:.4f}", f"{s.maximumValue:.4f}",
            f"{s.mean:.4f}", f"{s.stdDev:.4f}",
        ])

print("Wrote band_stats.csv")

Breakdown: The standard-library csv module writes a header then one row per band, formatting each statistic to four decimals. This produces a portable summary you can attach to a data delivery or diff between processing runs to catch regressions.

QGIS Version Compatibility

The code targets QGIS 3.34 LTR (Python 3.12).

The older native:zonalstatistics algorithm edited the input layer in place and is deprecated in favor of native:zonalstatisticsfb. Prefer the fb variant on all current versions; it leaves the source untouched and returns a clean output layer.

Troubleshooting

• All statistics return 0 or nan. The band has nodata where you expect data, or you used a 0-based band index. Use band 1 for the first band and confirm the nodata value with provider.sourceNoDataValue(1).
• Statistics seem wrong on a huge raster. You may have passed a sampleSize, giving an approximation. Drop the sample arguments for an exact full-scan result.
• Zonal output has empty stat columns. The polygons do not overlap the raster, or their CRS differs. Reproject the vector layer to the raster's CRS and confirm extents intersect.
• AttributeError on stats object. Use attributes (stats.mean) not method calls (stats.mean()). QgsRasterBandStatistics exposes plain fields.
• Histogram bins look skewed. You did not bound the range; pass stats.minimumValue and stats.maximumValue so bins span the real data range.

Conclusion

PyQGIS gives you raster statistics at two scales. bandStatistics reads whole-band min, max, mean, and standard deviation in one call, and histogram reveals the value distribution behind those numbers — both exact or sampled for speed. For per-feature summaries, native:zonalstatisticsfb writes count, mean, and other statistics into each polygon. Remembering that band indices are 1-based and that vector and raster must share a CRS is what keeps these results trustworthy across DEMs, imagery, and analysis pipelines.

Frequently Asked Questions

Are raster band indices 0-based or 1-based? 1-based. Band 1 is the first band in both bandStatistics(1, ...) and the RASTER_BAND parameter of zonal statistics. Using 0 raises an error or returns empty results.

How do I compute statistics faster on a large raster? Pass an extent and a sampleSize to bandStatistics. This reads a capped number of pixels and returns an approximation, which is far quicker than an exact full scan and usually accurate enough for ramp limits.

What is the difference between zonalstatistics and zonalstatisticsfb? The legacy native:zonalstatistics modifies the input vector layer in place. native:zonalstatisticsfb returns a new feature-based output and leaves the source untouched — the recommended choice on current QGIS.

Why are my zonal statistics columns empty? The polygons and the raster do not overlap, or they are in different coordinate reference systems. Reproject the vector layer to the raster's CRS and verify their extents intersect before running.

Related

• Raster Analysis Workflows
• Clip a Raster by a Mask Layer in PyQGIS
• Apply a Color Ramp to a Raster in PyQGIS