How to Calculate Population-Weighted Pollution Exposure in Python
On this page
Problem statement
A pollution map says how dirty each place is. An exposure figure says how dirty the air is where people are — and the two differ, because people are not spread evenly. A city's area-average NO2 counts its forests, lakes and airfields as heavily as its dense inner districts; weighting by population counts each resident once.
The calculation is a weighted mean of two rasters, and nearly all of the effort is in making them line up: a population grid at one resolution and CRS, a pollution field at another, and a boundary to cut both to. The number that comes out depends on the pollution field's resolution as much as on the air, which is worth knowing before quoting it.
Quick answer
Put the pollution field onto the population grid, then weight:
population-weighted mean = Σ (population × concentration) / Σ population
import io, zipfile, requests, numpy as np, pandas as pd, geopandas as gpd, osmnx as ox
import rasterio, rasterio.mask
ox.settings.http_user_agent = "spatialworkflow-examples/1.0"
berlin = ox.geocode_to_gdf("Berlin, Germany")
tile = ("https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_POP_GLOBE_R2023A/"
"GHS_POP_E2020_GLOBE_R2023A_4326_30ss/V1-0/tiles/"
"GHS_POP_E2020_GLOBE_R2023A_4326_30ss_V1_0_R4_C20.zip")
zipfile.ZipFile(io.BytesIO(requests.get(tile, timeout=300).content)).extractall("ghsl")
with rasterio.open("ghsl/GHS_POP_E2020_GLOBE_R2023A_4326_30ss_V1_0_R4_C20.tif") as src:
pop, pop_tf = rasterio.mask.mask(src, berlin.geometry, crop=True, nodata=0)
pop = np.where(pop[0] > 0, pop[0], 0.0)
rows, cols = np.nonzero(pop)
lon, lat = rasterio.transform.xy(pop_tf, rows, cols)
cells = pd.DataFrame({"lon": lon, "lat": lat, "pop": pop[rows, cols]})
print(f"GHSL 2020: {cells['pop'].sum():,.0f} people in {len(cells):,} cells of 30 arc-seconds")
# CAMS annual mean NO2 on a 0.1° lattice covering the city
lons = np.round(np.arange(13.05, 13.80, 0.1), 2)
lats = np.round(np.arange(52.35, 52.70, 0.1), 2)
pts = [(la, lo) for la in lats for lo in lons]
annual = {}
for i in range(0, len(pts), 20):
batch = pts[i:i + 20]
j = requests.get("https://air-quality-api.open-meteo.com/v1/air-quality", params={
"latitude": ",".join(str(p[0]) for p in batch),
"longitude": ",".join(str(p[1]) for p in batch),
"hourly": "nitrogen_dioxide", "start_date": "2023-01-01", "end_date": "2023-12-31",
"domains": "cams_europe"}, timeout=300).json()
for p, r in zip(batch, j):
annual[p] = np.nanmean(np.array(r["hourly"]["nitrogen_dioxide"], dtype=float))
# nearest lattice point for every population cell
lat_i = np.clip(np.round((cells.lat - lats[0]) / 0.1).astype(int), 0, len(lats) - 1)
lon_i = np.clip(np.round((cells.lon - lons[0]) / 0.1).astype(int), 0, len(lons) - 1)
cells["no2"] = [annual[(lats[a], lons[b])] for a, b in zip(lat_i, lon_i)]
pw = np.average(cells.no2, weights=cells["pop"])
print(f"area-weighted mean NO2: {cells.no2.mean():.1f} µg/m³")
print(f"population-weighted mean NO2: {pw:.1f} µg/m³")
print(f"people above the WHO annual guideline (10 µg/m³): "
f"{cells.loc[cells.no2 > 10, 'pop'].sum() / cells['pop'].sum():.0%}")GHSL 2020: 3,583,694 people in 1,609 cells of 30 arc-seconds
area-weighted mean NO2: 11.5 µg/m³
population-weighted mean NO2: 12.3 µg/m³
people above the WHO annual guideline (10 µg/m³): 86%
The population-weighted mean is higher than the area mean: people live where the emissions are.
A population grid, a shared boundary, resample, weight
1. Choose a population grid
GHSL (global, 1 km and 100 m, 1975–2030), WorldPop (global, 100 m), or a national census grid. Note the year and whether it counts residents (night-time) — they all do, which is itself an assumption about where exposure happens.
2. Clip both layers to the same boundary
Use the administrative boundary the result will be reported for. Cells cut by the boundary keep their full population in a simple mask; for small areas, weight by the covered fraction.
3. Resample the pollution field to the population grid
For each population cell, take the concentration of the pollution cell containing it (nearest for coarse model output, bilinear for smooth surfaces). Never resample population by interpolation — it is a count, and counts must be summed, not averaged.
4. Compute the weighted mean and the distribution
The weighted mean is one number; the share of people above each guideline or limit is often more useful, and it is sensitive to the tail that a mean hides.
5. Break it down by district
The same weights within each district give district exposures that can be compared with deprivation, age structure or health outcomes.
6. State the resolution of the pollution field
A 10 km model cannot see roads; a 10 m land-use regression can. The same people get different exposures from each, and the difference is not noise.
Berlin exposure by district
The data used on this page
The quick answer fetches everything: Berlin's boundary from OpenStreetMap, the GHSL 2020 population grid at 30 arc-seconds, and 2023 hourly NO2 from the CAMS European ensemble through the Open-Meteo air-quality API on a 0.1° lattice. The ten background-station means in Example 2 are Berlin's 2023 annual means from the EEA, as computed in How to interpolate air quality stations in Python. The setup points OSMnx at an Overpass mirror and reads the two OpenStreetMap answers this page uses as saved on 9 October 2026, so the district boundaries stay as they were when the numbers below were printed.
import io
import os
import osmnx as ox
import requests
import zipfile
# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
reply = requests.get("https://data.spatialworkflow.io/osm-cache-population-weighted-exposure-python-2026-10-09.zip", timeout=300)
reply.raise_for_status()
with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
archive.extractall("cache")
ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.requests_timeout = 600
ox.settings.max_query_area_size = 50e9 # one Overpass query for the whole cityExample 1 — by district
districts = ox.features_from_place("Berlin, Germany",
{"boundary": "administrative", "admin_level": "9"})
# osmnx ORs the tags: this returns every administrative boundary. Keep Berlin's 12 Bezirke.
districts = districts[(districts.admin_level == "9")
& districts.geom_type.isin(["Polygon", "MultiPolygon"])][["name", "geometry"]]
pts_gdf = gpd.GeoDataFrame(cells, geometry=gpd.points_from_xy(cells.lon, cells.lat), crs=4326)
joined = gpd.sjoin(pts_gdf, districts.to_crs(4326), predicate="within")
by = joined.groupby("name").apply(
lambda g: pd.Series({"people": g["pop"].sum(),
"pw_no2": np.average(g.no2, weights=g["pop"])}), include_groups=False)
print(f"{len(by)} districts")
print(by.sort_values("pw_no2", ascending=False).round({"people": -3, "pw_no2": 1}).to_string())12 districts
people pw_no2
name
Mitte 309000.0 13.9
Charlottenburg-Wilmersdorf 314000.0 13.2
Friedrichshain-Kreuzberg 199000.0 13.2
Pankow 406000.0 12.7
Tempelhof-Schöneberg 368000.0 12.4
Lichtenberg 241000.0 12.4
Reinickendorf 264000.0 12.2
Neukölln 300000.0 12.0
Steglitz-Zehlendorf 327000.0 11.9
Spandau 255000.0 11.3
Marzahn-Hellersdorf 275000.0 11.1
Treptow-Köpenick 327000.0 10.8
District differences at this resolution reflect where each district sits in the city's gradient — central districts higher, outer ones lower. Within-district differences between main roads and back streets are invisible to a 10 km model.
Example 2 — sensitivity to the pollution field
xy_cells = cells[["lon", "lat"]].to_numpy()
st_bg = pd.DataFrame({"lon": [13.349326, 13.348775, 13.430844, 13.418833, 13.529504,
13.483079, 13.368103, 13.225144, 13.647050, 13.296081],
"lat": [52.543041, 52.485814, 52.489451, 52.513606, 52.485296,
52.644169, 52.398406, 52.473192, 52.447697, 52.653269],
"no2": [17.5, 15.0, 15.6, 15.7, 12.1, 8.8, 8.3, 7.6, 7.5, 6.7]})
# 2023 annual means at Berlin's ten background stations, from the EEA guide
def idw(xk, zk, xq):
scale = np.array([np.cos(np.radians(52.5)), 1.0])
d = np.hypot(*((xq[:, None, :] - xk[None, :, :]) * scale).transpose(2, 0, 1))
w = 1 / np.maximum(d, 1e-4) ** 2
return (w * zk).sum(1) / w.sum(1)
surfaces = {"CAMS, 0.1° model": cells.no2.to_numpy(),
"IDW, 10 background stations": idw(st_bg[["lon", "lat"]].to_numpy(),
st_bg.no2.to_numpy(), xy_cells),
"one number, mean of stations": np.full(len(cells), st_bg.no2.mean())}
for name, v in surfaces.items():
print(f"{name:<31} population-weighted {np.average(v, weights=cells['pop']):5.1f}, "
f"people above 15 µg/m³: {cells['pop'][v > 15].sum() / cells['pop'].sum():4.0%}")CAMS, 0.1° model population-weighted 12.3, people above 15 µg/m³: 0%
IDW, 10 background stations population-weighted 12.9, people above 15 µg/m³: 17%
one number, mean of stations population-weighted 11.5, people above 15 µg/m³: 0%
The weighted means are close — 11.5 to 12.9 µg/m³; the share of people above 15 µg/m³ is 0% on two fields and 17% on the third. The tail of the exposure distribution — the part health impact assessments care about — depends entirely on how much spatial detail the pollution field has.
Example 3 — save the exposure grid
out = np.zeros_like(pop, dtype="float32")
out[rows, cols] = cells.no2.to_numpy()
with rasterio.open("berlin_no2_on_pop_grid.tif", "w", driver="GTiff", height=out.shape[0],
width=out.shape[1], count=2, dtype="float32", crs="EPSG:4326",
transform=pop_tf, nodata=0) as dst:
dst.write(pop.astype("float32"), 1)
dst.write(out, 2)
dst.set_band_description(1, "GHSL 2020 population")
dst.set_band_description(2, "CAMS 2023 annual NO2, µg/m³")
print(f"wrote {out.shape[1]} × {out.shape[0]} cells, population and NO2 on one grid")wrote 82 × 42 cells, population and NO2 on one grid
Keeping population and concentration on one grid, in one file, makes every later weighting — by age group, by district, by deprivation — a one-line calculation.
Why weighting raises the figure
Why population weighting raises the figure
Emissions and population are correlated: dense districts have more traffic, more heating and less open space. Weighting by where people live gives those districts their proper share, and the result is almost always higher than the area mean.
Why population must be summed, not interpolated
A population grid is a count per cell. Resampling it bilinearly to a finer grid would create people; to a coarser grid by averaging would lose them. The safe direction is to bring the concentration to the population grid and leave the population untouched.
Why residence is only a proxy
People spend most of the day indoors, and many spend it somewhere other than home. Residential exposure is the standard in epidemiology because it is measurable, and it systematically misses commuting and workplace exposure — which is often higher for people who travel along busy roads.
GHSL values, boundary cells and year mismatches
- GHSL values are floats (fractional people per cell) and include tiny values in unpopulated cells; threshold or keep, but do not round cell by cell.
- Boundary cells are fully counted by
rasterio.mask; for small areas useall_touched=Falseand area fractions. - CAMS domains:
cams_europeis about 0.1°;cams_globalis much coarser and should not be used for a city. - Year mismatch: population for 2020 and pollution for 2023 is normal; state both.
- CAMS assimilates observations, so it is not independent of the station network; that helps accuracy and rules it out as an independent check.
- Health impact: exposure × concentration–response function × baseline rates gives attributable cases; that step has its own assumptions.
Internal links
- Exposure assessment explained — what "exposure" means and its proxies
- How to build a land-use regression model for NO2 — a pollution field with road detail
- How to interpolate air quality station data — the IDW field above
- How to compare CAMS with monitoring stations — how good the model field is
- Population-weighted centroids and dasymetric mapping — population grids in general
- How to calculate zonal statistics in Python — the district breakdown with rasters
- How to estimate population exposure to heat — the same weighting for surface heat
FAQ
What is population-weighted exposure?
The average concentration experienced by residents: each cell's concentration weighted by the number of people living in it.
Why is it higher than the average concentration?
People live where emissions are highest; weighting by population gives those places their share.
Which population data should I use?
GHSL or WorldPop globally, or a national census grid; note the year and resolution.
Should I resample population to the pollution grid?
No — bring the pollution values to the population grid. Population is a count and must be summed, not interpolated.
Does the resolution of the pollution data matter?
Yes. The mean changes little; the share of people above a threshold can change a lot.
Is residential exposure the same as personal exposure?
No. It ignores time spent commuting, at work and indoors, but it is the standard measurable proxy.