How to Calculate Zonal Statistics in Python #

Problem statement #

You have a population raster and a layer of wards, and you need the total population in each ward. The question sounds trivial and the naive answer is wrong in at least three ways:

for ward in wards.itertuples():
    arr, _ = mask(src, [ward.geometry], crop=True)
    print(ward.name, arr.sum())

Ancoats     -3746109
Ardwick     -2914003

Negative population. The NoData sentinel is -9999 and it is being summed as if it were data.

Fix that and the next problem appears: the totals do not add up to the national figure, because cells straddling a ward boundary are counted whole for one ward and not at all for the other. Fix that and the third problem appears: on a 30 m raster, a ward smaller than a few cells returns nan, because no cell centre falls inside it.

Zonal statistics is the operation of summarising raster values within vector zones. Getting the summary right is easy. Getting the cells right is the actual work.

Quick answer #

Use rasterstats when you can:

from rasterstats import zonal_stats
import geopandas as gpd

wards = gpd.read_file("wards.gpkg").to_crs("EPSG:27700")

stats = zonal_stats(
    wards, "population_30m.tif",
    stats=["sum", "mean", "count", "min", "max"],
    nodata=-9999,
    all_touched=False,
)
wards = wards.join(gpd.pd.DataFrame(stats))
print(wards[["ward_name", "sum", "mean", "count"]].head())

        ward_name      sum    mean  count
0  Ancoats           14203    9.81   1448
1  Ardwick           18770   11.42   1644
2  Piccadilly         9114    7.03   1297

If rasterstats is not available, the pure-rasterio version is in Example 1 and is about fifteen lines.

Step-by-step solution #

1. Get both inputs into the same CRS #

import rasterio, geopandas as gpd

with rasterio.open("population_30m.tif") as src:
    raster_crs = src.crs

wards = gpd.read_file("wards.gpkg")
if wards.crs != raster_crs:
    wards = wards.to_crs(raster_crs)

Reproject the vector, not the raster — resampling a population raster changes the totals, which is precisely what you are trying to measure. See raster and vector do not line up when the results come back empty.

2. Decide what the raster values mean, because it changes the statistic #

The common error is summing a density raster. If cells hold people per square kilometre, the sum over a ward is not a population — it is people-per-km² multiplied by a cell count, which is a number with no meaning. Convert to counts first:

with rasterio.open("density_per_km2.tif") as src:
    cell_km2 = (src.res[0] * src.res[1]) / 1_000_000
    counts = src.read(1, masked=True) * cell_km2      # now summable

3. Choose the cell-inclusion rule #

By default a cell belongs to a zone if the cell centre falls inside the polygon. Each cell therefore belongs to at most one zone, so sums across zones do not double-count. That property is worth a lot, and it is why all_touched=False is the right default for totals.

all_touched=True includes every cell the polygon touches. Cells on a shared boundary then belong to both neighbouring zones, and the sum of all zones exceeds the raster total.

totals = zonal_stats(wards, RASTER, stats=["sum"], nodata=-9999)
touched = zonal_stats(wards, RASTER, stats=["sum"], nodata=-9999, all_touched=True)

print(sum(s["sum"] for s in totals))     # 1_842_003
print(sum(s["sum"] for s in touched))    # 2_104_887  ← 14% over

Use all_touched=True only when zones are small relative to the cell size and missing them entirely is the worse error — and never when the totals must reconcile.

4. Handle zones that contain no cells #

stats = zonal_stats(wards, RASTER, stats=["sum", "count"], nodata=-9999)
empty = [i for i, s in enumerate(stats) if not s["count"]]
print(f"{len(empty)} zones with no cells")

A zone smaller than a cell, or one that falls entirely on NoData, returns count = 0 and sum = None. Three ways to respond, and the right one depends on the question:

# (a) treat as zero — correct for a count raster over genuinely empty land
df["sum"] = df["sum"].fillna(0)

# (b) fall back to all_touched for the small ones only
small = wards.loc[empty]
patched = zonal_stats(small, RASTER, stats=["sum"], nodata=-9999, all_touched=True)

# (c) keep as missing — correct when "no data here" is not "zero here"

Never silently apply (a) to a mean or an elevation. "No cells" and "zero" are different statements, and the null, empty, missing and invalid distinction applies to values just as much as to geometry.

5. Join the results back and sanity-check #

import pandas as pd

df = pd.DataFrame(stats)
wards = wards.reset_index(drop=True).join(df)

raster_total = None
with rasterio.open(RASTER) as src:
    raster_total = src.read(1, masked=True).sum()

print(f"raster total      {raster_total:,.0f}")
print(f"sum of zones      {wards['sum'].sum():,.0f}")
print(f"cells in zones    {wards['count'].sum():,}")
print(f"zones with 0 cells{(wards['count'] == 0).sum():>6}")

raster total      1,904,551
sum of zones      1,842,003
cells in zones    198,447
zones with 0 cells     3

The zone total should be at or below the raster total — the difference is cells outside every zone. Above it means double-counting, which means all_touched=True or overlapping zones. Both are worth knowing about before the number goes in a report.

Code examples #

Example 1: zonal statistics with rasterio alone #

No extra dependency, and it makes every decision visible:

import numpy as np, rasterio, geopandas as gpd
from rasterio.mask import mask

def zonal(raster_path, gdf, key_col, all_touched=False):
    rows = []
    with rasterio.open(raster_path) as src:
        gdf = gdf.to_crs(src.crs)
        for feat in gdf.itertuples():
            try:
                arr, _ = mask(src, [feat.geometry], crop=True,
                              filled=False, all_touched=all_touched)
            except ValueError:                      # no overlap at all
                rows.append({key_col: getattr(feat, key_col), "count": 0})
                continue
            vals = arr[0].compressed()              # drops NoData
            rows.append({
                key_col: getattr(feat, key_col),
                "count": int(vals.size),
                "sum":   float(vals.sum()) if vals.size else None,
                "mean":  float(vals.mean()) if vals.size else None,
                "min":   float(vals.min()) if vals.size else None,
                "max":   float(vals.max()) if vals.size else None,
                "p90":   float(np.percentile(vals, 90)) if vals.size else None,
            })
    return gpd.pd.DataFrame(rows)

wards = gpd.read_file("wards.gpkg")
print(zonal("elevation_25m.tif", wards, "ward_name").head())

filled=False returns a masked array, and .compressed() gives only the valid cells as a flat array. That one line is the whole NoData fix: the sentinel never reaches the statistic.

This version is slower than rasterstats for thousands of small zones, because it opens a window per zone. For a few hundred zones the difference is irrelevant, and the explicitness is worth more.

Example 2: categorical zonal statistics — area per class #

For a land-cover raster the useful output is not a mean but a breakdown:

import numpy as np, rasterio, geopandas as gpd, pandas as pd
from rasterio.mask import mask

CLASSES = {1: "urban", 2: "arable", 3: "grassland", 4: "woodland", 5: "water"}

def class_areas(raster_path, gdf, key_col):
    rows = []
    with rasterio.open(raster_path) as src:
        gdf = gdf.to_crs(src.crs)
        cell_area = abs(src.res[0] * src.res[1])
        for feat in gdf.itertuples():
            arr, _ = mask(src, [feat.geometry], crop=True, filled=False)
            vals = arr[0].compressed()
            codes, counts = np.unique(vals, return_counts=True)
            row = {key_col: getattr(feat, key_col), "cells": int(vals.size)}
            for code, n in zip(codes, counts):
                row[f"{CLASSES.get(int(code), code)}_ha"] = n * cell_area / 10_000
            rows.append(row)
    return pd.DataFrame(rows).fillna(0)

areas = class_areas("landcover_25m.tif", gpd.read_file("wards.gpkg"), "ward_name")
print(areas.head())

     ward_name  cells  urban_ha  arable_ha  grassland_ha  woodland_ha  water_ha
0    Ancoats     1448     78.25      3.12          6.94         2.31      0.88
1    Ardwick     1644     84.50      1.50         11.06         4.13      0.44

Two things make this correct. Areas come from the cell count times cell area, which requires a projected CRS — do this in degrees and res is in degrees and the hectares are fiction. And np.unique on the compressed array means NoData never becomes a class.

Example 3: area-weighted statistics for small zones #

When zones are comparable in size to cells, whole-cell assignment is too crude. Weight each cell by the fraction of it inside the zone:

import numpy as np, rasterio, geopandas as gpd
from rasterio.features import rasterize
from shapely.geometry import box

def area_weighted_mean(raster_path, geom, subsample=8):
    """Weight cells by how much of each falls inside the polygon."""
    with rasterio.open(raster_path) as src:
        from rasterio.windows import from_bounds
        win = from_bounds(*geom.bounds, transform=src.transform).round_lengths().round_offsets()
        data = src.read(1, window=win, masked=True)
        transform = src.window_transform(win)

    # rasterise the polygon at N× resolution, then average down to get coverage
    fine = rasterize(
        [(geom, 1)],
        out_shape=(data.shape[0] * subsample, data.shape[1] * subsample),
        transform=transform * transform.scale(1 / subsample, 1 / subsample),
        fill=0, dtype="uint8",
    )
    coverage = fine.reshape(
        data.shape[0], subsample, data.shape[1], subsample
    ).mean(axis=(1, 3))                       # 0.0 … 1.0 per cell

    weights = np.where(data.mask, 0, coverage)
    if weights.sum() == 0:
        return None
    return float((data.filled(0) * weights).sum() / weights.sum())

wards = gpd.read_file("wards.gpkg").to_crs("EPSG:27700")
print(area_weighted_mean("elevation_25m.tif", wards.geometry.iloc[0]))

The trick is rasterising the polygon eight times finer than the data grid and averaging each 8 × 8 block, which gives each cell a coverage fraction between 0 and 1 at 1/64-cell precision. A boundary cell 30% inside the zone contributes 30% of its weight instead of 100% or nothing.

This costs roughly subsample² times the memory of the window, so subsample=8 is a reasonable ceiling. It matters most where it is most often skipped: small zones, coarse rasters, and any total that has to reconcile across neighbours.

Explanation #

Zonal statistics is a spatial join followed by a group-by, and it inherits its difficulties from the join half rather than the statistics half.

The join is between two incompatible geometries. A raster cell is a fixed square from an arbitrary grid; a zone is an arbitrary polygon. They almost never nest. So every zonal operation must answer a question that has no exact answer: does this cell, 40% inside the zone, belong to it?

Three answers are in common use, and each is right for a different question:

• Centre-in-polygon (the default) is a partition: every cell belongs to exactly one zone, so totals across zones never double-count. It is the correct choice whenever the numbers must reconcile. Its cost is that small zones can receive no cells at all.
• All-touched is a cover: every zone gets every cell it overlaps, so nothing is missed, but boundary cells are counted more than once. Correct when the question is "what is in here" and wrong when the question is "how much of the total is here".
• Area weighting is the honest answer and the expensive one. Each cell contributes in proportion to its overlap, so totals reconcile and small zones get sensible values. It costs a finer rasterisation per zone.

The second structural point is that the statistic must match the quantity. This is not a GIS issue at all — it is the ordinary distinction between extensive quantities (population, area, rainfall volume) which are additive, and intensive ones (density, temperature, slope) which are not. A raster gives no clue which it holds; the value 9.81 is just a number. Summing an intensive raster is the most common error in this whole area, and it produces results with the right order of magnitude, which is what makes it survive review.

Third, NoData must be excluded before the arithmetic, not after. This sounds obvious and yet it is the single most frequent bug, because mask() returns filled arrays by default and -9999 is a perfectly ordinary float to NumPy. filled=False plus .compressed(), or rasterstats's nodata= argument, is the fix. The tell-tale symptom is a negative sum or a mean pulled implausibly low — the same class of error described in the raster data model.

Finally, note what zonal statistics does not do: it does not resample. The cell values used are exactly the source values. This is deliberate and correct — resampling before summarising would invent data and change totals. If your zones and raster do not align, the answer is a better inclusion rule, not a resampled raster.

Edge cases or notes #

• rasterstats reads nodata from the file if you do not pass it, but only when the tag is set. If src.nodata is None, pass the value explicitly.
• Overlapping zones double-count by design. Check with gdf.geometry.overlaps before totalling — see how to find and fix gaps and overlaps in a polygon coverage.
• Zones must be valid. Self-intersecting polygons give wrong or empty results; run the validity fix first.
• Percentiles and medians cannot be combined across tiles, so a tiled zonal job cannot compute them from partials. Sums, counts, mins and maxes can.
• zonal_stats accepts a path, a GeoDataFrame or GeoJSON-like dicts, but a GeoDataFrame must already be in the raster's CRS — it does not reproject for you.
• stats="*" computes everything, including expensive ones. Ask for what you need.
• Area from cell counts requires a projected CRS. In EPSG:4326 the cell area varies with latitude and the number is meaningless.
• Multi-band rasters need band= (rasterstats) or a per-band loop; the default is band 1.
• categorical=True in rasterstats returns a class-count dict per zone, which is Example 2 in one argument.
• gdal_polygonize + a SQL group-by does the same job in PostGIS when the raster already lives there.

Internal links #

• The raster data model explained — NoData, dtype and why sums go negative
• How to clip a raster to a polygon in Python — extracting rather than summarising
• How to extract raster values at point locations — the point equivalent
• Raster resampling explained — why you should not resample before summarising
• How to aggregate spatial data by region in GeoPandas — the all-vector version of the same question
• Raster and vector do not line up in Python — when every zone comes back empty
• How to rasterize a vector layer in Python — the operation behind area weighting
• How to count points in polygons with GeoPandas — the same join, different data model

FAQ #

Why is my zonal sum negative? #

NoData is being included. Pass nodata= to zonal_stats, or use filled=False and .compressed() with rasterio.mask.

Should I use all_touched=True? #

Only when zones are small relative to cells and missing them entirely is worse than double-counting. It breaks reconciliation: the sum of all zones will exceed the raster total.

Why do some zones return None? #

No cell centre fell inside them — the zone is smaller than a cell, or it sits entirely on NoData. Retry those zones with all_touched=True, or use area weighting.

Should I sum or average? #

Sum counts, average rates. If the raster holds people per cell, sum. If it holds people per square kilometre, convert to counts first or take the mean and treat it as a mean density.

Do I need rasterstats? #

No. The fifteen-line rasterio.mask version in Example 1 does the same job and makes every decision explicit. rasterstats is faster across thousands of small zones.

How do I get the area of each land-cover class per zone? #

Count cells per class and multiply by cell area, as in Example 2 — or pass categorical=True to zonal_stats. Both need a projected CRS for the area to mean anything.

Why do my zone totals not add up to the raster total? #

Cells outside every zone are excluded, which makes the zone total lower. A zone total that is higher means double-counting from overlapping zones or all_touched=True.