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. Here the raster is WorldPop's 2020 population estimate for the UK on a 1 km grid, and the zones are Leeds's 33 wards (the setup is under The data used on this page). The question sounds trivial and the naive answer is wrong in at least three ways:

import geopandas as gpd
import rasterio
from rasterio.mask import mask

wards = gpd.read_file("wards.gpkg")
with rasterio.open("population_1km.tif") as src:
    for ward in wards.to_crs(src.crs).head(2).itertuples():
        arr, _ = mask(src, [ward.geometry], crop=True)
        print(f"{ward.ward_name:24} {arr.sum():>14,.0f}")

Ardsley & Robin Hood         -6,674,615
Armley                       -2,671,290

Negative population. The NoData sentinel is -99999, mask fills every cell outside the ward with it, and it is being summed as if it were data.

Fix that and the next problem appears: the totals do not add up, 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: a zone smaller than a cell can contain no cell centre at all, and comes back empty.

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:

import pandas as pd
from rasterstats import zonal_stats

with rasterio.open("population_1km.tif") as src:
    zones = wards.to_crs(src.crs)
stats = zonal_stats(zones, "population_1km.tif", stats=["sum", "mean", "count"], all_touched=False)
wards = wards.join(pd.DataFrame(stats))
print(wards[["ward_name", "sum", "mean", "count"]].head(3).round(1).to_string(index=False))

            ward_name     sum   mean  count
 Ardsley & Robin Hood 25317.8  684.3     37
               Armley 28683.4 2206.4     13
Bramley & Stanningley 18575.5 1857.6     10

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

One CRS, the meaning of the values, the right statistic

1. Get both inputs into the same CRS

with rasterio.open("population_1km.tif") as src:
    raster_crs = src.crs
wards = gpd.read_file("wards.gpkg")
print(f"raster {raster_crs}, wards {wards.crs}")
if wards.crs != raster_crs:
    wards = wards.to_crs(raster_crs)

raster EPSG:4326, wards EPSG:27700

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

WorldPop's cells hold people per cell, so sum is right here. 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       # a projected CRS in metres
    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:

centres = zonal_stats(wards, "population_1km.tif", stats=["sum"])
touched = zonal_stats(wards, "population_1km.tif", stats=["sum"], all_touched=True)
a, b = sum(s["sum"] for s in centres), sum(s["sum"] for s in touched)
print(f"cell centres  {a:>12,.0f}")
print(f"all touched   {b:>12,.0f}  ({b / a - 1:+.0%})")

cell centres       788,164
all touched      1,388,589  (+76%)

On a 1 km grid over wards of a few square kilometres, a large share of every ward's cells sit on a boundary — which is why all_touched=True overstates the city's population by so much. Use it 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

Leeds's 488 LSOAs — census zones of roughly 1,500 people each — are often smaller than a 1 km cell:

lsoas = gpd.read_file("lsoas.gpkg").to_crs(raster_crs)
lsoa_stats = zonal_stats(lsoas, "population_1km.tif", stats=["sum", "count"])
empty = [i for i, s in enumerate(lsoa_stats) if not s["count"]]

print(f"{len(empty)} of {len(lsoas)} LSOAs contain no cell centre; their sum is {lsoa_stats[empty[0]]['sum']}")
area = lsoas.to_crs(27700).area / 1e6
is_empty = lsoas.index.isin(empty)
print(f"median area: empty LSOAs {area[is_empty].median():.2f} km², the rest {area[~is_empty].median():.2f} km²")

144 of 488 LSOAs contain no cell centre; their sum is None
median area: empty LSOAs 0.29 km², the rest 0.57 km²

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
patched = zonal_stats(lsoas.iloc[empty], RASTER, stats=["sum"], all_touched=True)

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

Here (a) would be wrong. LSOAs are drawn to hold similar numbers of people, so the small ones are small because they are dense — the empty ones are half the area of the rest — and filling them with zero would erase some of the most populated neighbourhoods in the city. Never silently apply (a) to a mean or an elevation either. "No cells" and "zero" are different statements, and the null, empty, missing and invalid distinction applies to values just as much as to geometry. Example 3 shows the better answer for small zones.

5. Join the results back and sanity-check

Check the total against something you trust. The 2021 Census counted every ward:

census = pd.read_csv("census_2021_wards.csv")
check = wards[["ward_code"]].assign(worldpop=[s["sum"] for s in centres]).merge(census, on="ward_code")
print(f"WorldPop 2020, cell centres  {check['worldpop'].sum():>10,.0f}")
print(f"Census 2021                  {check['population'].sum():>10,.0f}")
print(f"largest ward difference      {(check['worldpop'] / check['population'] - 1).abs().max():>10.0%}")

WorldPop 2020, cell centres     788,164
Census 2021                     811,952
largest ward difference             23%

The city totals agree to within a few per cent — a modelled estimate for 2020 against a count a year later — and individual wards differ far more, because on a 1 km grid most of a ward's population sits in boundary cells assigned whole to one side or the other. The sum of zones should be at or below the raster's own total over the same area; 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.

Zonal statistics with rasterio alone, and categories

The data used on this page

• population_1km.tif — WorldPop's 2020 unconstrained population estimate for the United Kingdom, aggregated to 30 arc-seconds (about 1 km), people per cell, float32 with NoData -99999, cut to a window around Leeds.
• wards.gpkg and lsoas.gpkg — Leeds's 33 electoral wards (ONS, December 2024) and 488 LSOAs (ONS, 2021), in British National Grid.
• census_2021_wards.csv — usual residents per ward from the 2021 Census, via Nomis.
• elevation_30m.tif — the Copernicus 30 m DEM over Leeds, for example 1.
• landcover_10m_bng.tif — ESA WorldCover 2021 over Leeds, reprojected to British National Grid at 10 m with nearest-neighbour resampling, for example 2.

import io
from pathlib import Path

import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio
import requests
from rasterio.warp import Resampling, calculate_default_transform, reproject
from rasterio.windows import from_bounds

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
wards = gpd.read_file(f"{ONS}/Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
                      "?where=LAD24NM%3D%27Leeds%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson")
wards = (wards.set_crs(27700, allow_override=True).rename(columns={"WD24CD": "ward_code", "WD24NM": "ward_name"})
              .sort_values("ward_code").reset_index(drop=True))
wards.to_file("wards.gpkg")

features, offset = [], 0
while True:
    page = requests.get(f"{ONS}/Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query",
                        params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "LSOA21CD", "outSR": 27700,
                                "f": "geojson", "resultOffset": offset}, timeout=120).json()
    features += page["features"]
    if not page.get("properties", {}).get("exceededTransferLimit"):
        break
    offset += len(page["features"])
gpd.GeoDataFrame.from_features(features, crs=27700).to_file("lsoas.gpkg")

pop = pd.read_csv(io.StringIO(requests.get("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv", timeout=120,
                  params={"date": "latest", "geography": ",".join(wards["ward_code"]), "c2021_restype_3": 0,
                          "measures": 20100, "select": "geography_code,obs_value"}).text))
pop.columns = ["ward_code", "population"]
pop.to_csv("census_2021_wards.csv", index=False)

west, south, east, north = wards.to_crs(4326).total_bounds
pad = 0.05
Path("worldpop.tif").write_bytes(requests.get(
    "https://data.worldpop.org/GIS/Population/Global_2000_2020_1km/2020/GBR/gbr_ppp_2020_1km_Aggregated.tif",
    timeout=300).content)


def cut(src_path, dst_path):
    with rasterio.open(src_path) as src:
        window = from_bounds(west - pad, south - pad, east + pad, north + pad,
                             transform=src.transform).round_offsets().round_lengths()
        profile = src.profile | {"driver": "GTiff", "width": int(window.width), "height": int(window.height),
                                 "transform": src.window_transform(window), "compress": "deflate"}
        with rasterio.open(dst_path, "w", **profile) as dst:
            dst.write(src.read(window=window))


cut("worldpop.tif", "population_1km.tif")
cut("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_W002_00_DEM/"
    "Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif", "elevation_30m.tif")
cut("https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/ESA_WorldCover_10m_2021_v200_N51W003_Map.tif",
    "landcover_wgs84.tif")
with rasterio.open("landcover_wgs84.tif") as src:
    t, w, h = calculate_default_transform(src.crs, "EPSG:27700", src.width, src.height, *src.bounds, resolution=10)
    with rasterio.open("landcover_10m_bng.tif", "w", **(src.profile | {"crs": "EPSG:27700", "transform": t,
                                                                         "width": w, "height": h})) as dst:
        reproject(rasterio.band(src, 1), rasterio.band(dst, 1), resampling=Resampling.nearest)
print(f"{len(wards)} wards, {len(features)} LSOAs, census total {pop['population'].sum():,}")

33 wards, 488 LSOAs, census total 811,952

Example 1: zonal statistics with rasterio alone

No extra dependency, and it makes every decision visible:

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 and the cells outside the zone
            rows.append({
                key_col: getattr(feat, key_col),
                "count": int(vals.size),
                "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 pd.DataFrame(rows)


elevation = zonal("elevation_30m.tif", gpd.read_file("wards.gpkg"), "ward_name")
print(elevation.sort_values("mean", ascending=False).head(3).round(0).to_string(index=False))

        ward_name  count  mean   min   max   p90
Guiseley & Rawdon  24155 191.0 109.0 341.0 266.0
   Otley & Yeadon  27456 145.0  42.0 280.0 223.0
        Alwoodley  23190 139.0  79.0 185.0 157.0

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:

CLASSES = {10: "tree", 20: "shrub", 30: "grass", 40: "crop", 50: "built", 60: "bare", 80: "water", 90: "wetland"}


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])            # m², so the CRS must be projected
        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_10m_bng.tif", gpd.read_file("wards.gpkg"), "ward_name")
cols = ["ward_name", "built_ha", "grass_ha", "tree_ha", "crop_ha", "water_ha"]
print(areas.sort_values("built_ha", ascending=False)[cols].head(3).round(0).to_string(index=False))
print(f"all wards: {areas.filter(like='_ha').sum().sum():,.0f} ha classified, "
      f"against {gpd.read_file('wards.gpkg').area.sum() / 10_000:,.0f} ha of ward area")
print("Leeds by class (ha):", areas.filter(like="_ha").sum().sort_values(ascending=False).round(0).astype(int).to_dict())

                  ward_name  built_ha  grass_ha  tree_ha  crop_ha  water_ha
        Hunslet & Riverside     650.0      55.0    163.0      4.0      13.0
          Beeston & Holbeck     512.0      73.0    190.0      1.0       0.0
Burmantofts & Richmond Hill     501.0     143.0    218.0      8.0       3.0
all wards: 55,170 ha classified, against 55,170 ha of ward area
Leeds by class (ha): {'grass_ha': 17184, 'tree_ha': 15520, 'built_ha': 11510, 'crop_ha': 10359, 'water_ha': 465, 'bare_ha': 93, 'wetland_ha': 38}

Two things make this correct. Areas come from the cell count times cell area, which requires a projected CRS — WorldCover arrives in degrees, and in degrees res is in degrees and the hectares are fiction, which is why the setup reprojected it (with nearest-neighbour resampling, so no class codes were invented). And np.unique on the compressed array means NoData never becomes a class. The classified area matching the ward area is the check that nothing was lost.

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:

from rasterio.features import geometry_window, rasterize


def area_weighted_sum(raster_path, geom, subsample=10):
    """Sum of cell values, each weighted by the fraction of the cell inside the polygon."""
    with rasterio.open(raster_path) as src:

        win = geometry_window(src, [geom])                   # whole cells, rounded outward
        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))
    return float((data.filled(0) * np.where(data.mask, 0, coverage)).sum())


with rasterio.open("population_1km.tif") as src:
    lsoa_zones = gpd.read_file("lsoas.gpkg").to_crs(src.crs)
    ward_zones = gpd.read_file("wards.gpkg").to_crs(src.crs)
weighted_lsoa = [area_weighted_sum("population_1km.tif", g) for g in lsoa_zones.geometry]
weighted_ward = [area_weighted_sum("population_1km.tif", g) for g in ward_zones.geometry]
print(f"LSOAs with no people:  centres {sum(1 for s in lsoa_stats if not s['sum'])}, "
      f"area-weighted {sum(1 for v in weighted_lsoa if v == 0)}")
print(f"sum over the 488 LSOAs: centres {sum(s['sum'] or 0 for s in lsoa_stats):,.0f}, "
      f"area-weighted {sum(weighted_lsoa):,.0f}")
print(f"sum over the 33 wards:  centres {sum(s['sum'] for s in centres):,.0f}, "
      f"area-weighted {sum(weighted_ward):,.0f}")

LSOAs with no people:  centres 144, area-weighted 0
sum over the 488 LSOAs: centres 788,164, area-weighted 786,393
sum over the 33 wards:  centres 788,164, area-weighted 786,303

The trick is rasterising the polygon ten times finer than the data grid and averaging each 10 × 10 block, which gives each cell a coverage fraction between 0 and 1 at 1/100-cell precision. A boundary cell 30% inside the zone contributes 30% of its weight instead of 100% or nothing, so every LSOA gets its own share. Notice what does not change: the city total. Centre assignment conserves it too — it simply hands the 144 small LSOAs' people to whichever neighbour holds the cell centre — so a correct total is no evidence that the per-zone figures are right.

This costs roughly subsample² times the memory of the window. exactextract (pip install exactextract) does the same job exactly rather than by subsampling, and is the better tool when you do this at scale.

A spatial join followed by a group-by

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.

Nodata in rasterstats and overlapping zones

• 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
• How to extract a time series per polygon from a NetCDF grid — zonal statistics with a time axis
• How to summarise catchment attributes — zones delineated from a DEM
• How to convert SAR backscatter to decibels — why zonal means of radar data need linear units
• How to estimate wind energy with a Weibull distribution — summarising a wind resource properly
• How to calculate population-weighted exposure — weighting a raster by population
• How to calculate forest loss by year for a region — zonal sums with area-weighted pixels
• How to summarise land surface temperature by neighbourhood — district heat statistics

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.