How to Build a Land-Use Regression Model for NO2 in Python
On this page
Problem statement
Interpolating monitoring stations maps the urban background and misses the roads. A dispersion model maps the roads and needs emission inventories, traffic flows and meteorology that most analysts do not have. Land-use regression sits between the two: it explains each station's annual mean by what surrounds it — how much road is nearby, how far it is from the centre — and applies the fitted relationship everywhere those predictors can be computed.
The method is simple and easy to overfit. A city's network has tens of stations, the candidate predictors number in the dozens (every road class at every buffer radius), and a model chosen by searching them all will fit the stations beautifully and predict the rest of the city badly. The discipline is in the validation.
Quick answer
Compute predictors at each station, fit a small linear model, and validate by leaving each station out:
import io, requests, numpy as np, pandas as pd, geopandas as gpd, osmnx as ox
from shapely.geometry import Point
API = "https://eeadmz1-downloads-api-appservice.azurewebsites.net/ParquetFile/urls"
META = "https://discomap.eea.europa.eu/map/fme/metadata/PanEuropean_metadata.csv"
r = requests.post(API, json={"countries": ["DE"], "cities": ["Berlin"], "pollutants": ["NO2"],
"dataset": 2, "dateTimeStart": "2023-01-01T00:00:00Z",
"dateTimeEnd": "2023-12-31T23:00:00Z", "aggregationType": "hour"},
timeout=120)
rows = []
for url in [u.strip() for u in r.text.splitlines()[1:] if u.strip()]:
d = pd.read_parquet(io.BytesIO(requests.get(url, timeout=120).content))
d["Value"] = pd.to_numeric(d.Value, errors="coerce")
y = d[(d.Start >= "2023-01-01") & (d.Start < "2024-01-01") & (d.Validity > 0)]
if len(y) >= 0.9 * 8760:
rows.append((d.Samplingpoint.iloc[0].split("/")[1], y.Value.mean()))
meta = pd.read_csv(META, sep="\t", low_memory=False,
storage_options={"User-Agent": "spatialworkflow-examples/1.0"}).drop_duplicates("SamplingPoint") # EEA refuses Python's default agent
st = pd.DataFrame(rows, columns=["SamplingPoint", "no2"]).merge(meta, on="SamplingPoint")
st = gpd.GeoDataFrame(st, geometry=gpd.points_from_xy(st.Longitude, st.Latitude),
crs=4326).to_crs(25833)
ox.settings.http_user_agent = "spatialworkflow-examples/1.0"
roads = ox.features_from_bbox((13.08, 52.33, 13.77, 52.68),
{"highway": ["motorway", "trunk", "primary", "secondary"]})
roads = roads[roads.geom_type.isin(["LineString", "MultiLineString"])].to_crs(25833)
centre = gpd.GeoSeries([Point(13.4133, 52.5219)], crs=4326).to_crs(25833).iloc[0] # Alexanderplatz
st["road_50m"] = [roads.clip(p.buffer(50)).length.sum() for p in st.geometry]
st["km_centre"] = st.distance(centre) / 1000
X = np.c_[np.ones(len(st)), st.road_50m, st.km_centre]
yv = st.no2.to_numpy()
beta, *_ = np.linalg.lstsq(X, yv, rcond=None)
pred = X @ beta
loo = np.array([np.r_[1, st.road_50m.iloc[i], st.km_centre.iloc[i]] @
np.linalg.lstsq(np.delete(X, i, 0), np.delete(yv, i), rcond=None)[0]
for i in range(len(st))])
r2 = lambda p: 1 - ((yv - p) ** 2).sum() / ((yv - yv.mean()) ** 2).sum()
print(f"NO2 = {beta[0]:.1f} + {beta[1]:.4f} × metres of road within 50 m "
f"{beta[2]:+.2f} × km from the centre")
print(f"fit on all {len(st)} stations: R² {r2(pred):.2f}, RMSE {np.sqrt(((yv - pred) ** 2).mean()):.1f}")
print(f"leave-one-out: R² {r2(loo):.2f}, RMSE {np.sqrt(((yv - loo) ** 2).mean()):.1f} µg/m³")NO2 = 20.0 + 0.0603 × metres of road within 50 m -0.74 × km from the centre
fit on all 17 stations: R² 0.77, RMSE 3.6
leave-one-out: R² 0.68, RMSE 4.3 µg/m³
Two predictors, one for the road increment and one for the city-scale gradient. The leave-one-out numbers are the ones to quote; the fitted R² is always better.
Every station, predictors at several scales, a small model
1. Use every station, of every type
Unlike interpolation, LUR wants the kerbside stations: they are the observations that tie high concentrations to road predictors. A model fitted only to background stations has never seen a road.
2. Compute candidate predictors at several scales
Road length within 25, 50, 100, 300, 500 and 1,000 m, by road class; distance to the nearest major road; population or built-up area within buffers; distance to the centre. Each is a hypothesis about how pollution disperses.
3. Keep the model small
A rule of thumb is one predictor per ten sites. With seventeen stations, two predictors is already generous. Choose them by physical sense and by leave-one-out performance, not by maximising the fitted R².
4. Check the signs
Road length should increase NO2; distance from the centre should decrease it. A predictor with the "wrong" sign is fitting noise and will produce nonsense when applied across the city.
5. Validate by leaving stations out
Leave-one-out for small networks; spatially blocked cross-validation for larger ones. If predictor selection was part of fitting, repeat the selection inside each fold, or the validation is optimistic.
6. Predict on a grid, within the range of the stations
Compute the same predictors on a grid and apply the model. Where a grid cell's predictors fall outside the range seen at stations — a motorway interchange, say — the prediction is an extrapolation and should be capped or flagged.
Buffer radii scanned honestly, three methods compared
The data used on this page
Berlin's 2023 hourly NO₂ from the European Environment Agency's air-quality download service and station metadata, as in the quick answer; the roads come from OpenStreetMap through osmnx, via the maps.mail.ru Overpass mirror while overpass-api.de was refusing connections. The setup only points osmnx at the mirror. The OpenStreetMap answers are read as the Overpass mirror gave them on 10 October 2026.
import io
import os
import osmnx as ox
import requests
import zipfile
# OpenStreetMap as Overpass returned it on 10 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-land-use-regression-no2-python-2026-10-10.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" # overpass-api.de refused us
print(ox.__version__)2.1.1
Example 1 — scanning buffer radii honestly
for radius in (25, 50, 100, 300, 1000):
x = np.array([roads.clip(p.buffer(radius)).length.sum() for p in st.geometry])
Xr = np.c_[np.ones(len(st)), x, st.km_centre]
loo_r = np.array([Xr[i] @ np.linalg.lstsq(np.delete(Xr, i, 0), np.delete(yv, i),
rcond=None)[0] for i in range(len(st))])
print(f"road length within {radius:>4} m: r = {np.corrcoef(x, yv)[0, 1]:+.2f}, "
f"LOO RMSE with distance-to-centre {np.sqrt(((yv - loo_r) ** 2).mean()):.1f}")road length within 25 m: r = +0.61, LOO RMSE with distance-to-centre 5.3
road length within 50 m: r = +0.75, LOO RMSE with distance-to-centre 4.3
road length within 100 m: r = +0.68, LOO RMSE with distance-to-centre 5.6
road length within 300 m: r = +0.49, LOO RMSE with distance-to-centre 6.5
road length within 1000 m: r = +0.64, LOO RMSE with distance-to-centre 6.3
50 m is the sweet spot: 25 m misses roads whose mapped centreline lies a little further from the monitor than the kerb does, and wider buffers start measuring the neighbourhood instead of the street. Picking the radius from this table and then quoting its LOO error is still slightly optimistic — the choice was made on the same stations — which is why five candidates is better than fifty.
Example 2 — the same stations, three methods
bg = (st.AirQualityStationType == "background").to_numpy()
xy = np.c_[st.geometry.x, st.geometry.y]
def idw(xk, zk, xq):
d = np.hypot(*(xq[:, None, :] - xk[None, :, :]).transpose(2, 0, 1))
w = 1 / np.maximum(d, 1.0) ** 2
return (w * zk).sum(1) / w.sum(1)
idw_loo = np.array([idw(np.delete(xy, i, 0)[np.delete(bg, i)], np.delete(yv, i)[np.delete(bg, i)],
xy[[i]])[0] for i in range(len(st))])
for label, p in (("IDW of background stations", idw_loo), ("land-use regression", loo)):
e = p - yv
print(f"{label:<27} RMSE at traffic sites {np.sqrt((e[~bg] ** 2).mean()):5.1f}, "
f"at background sites {np.sqrt((e[bg] ** 2).mean()):4.1f}")IDW of background stations RMSE at traffic sites 11.2, at background sites 4.0
land-use regression RMSE at traffic sites 6.2, at background sites 2.0
The background interpolation misses the kerbside stations by 11.2 µg/m³ on average — it has no idea a road is there. The regression halves that, because the road is one of its inputs, and it also beats interpolation at the background stations (2.0 against 4.0): the distance-to-centre term captures the city gradient better than ten scattered points do.
Example 3 — a prediction map
import rasterio.features
from rasterio.transform import from_origin
from scipy.signal import fftconvolve
x0, y0 = 385000, 5810000
res, n = 10, 2400 # a 24 km square, 10 m cells
tf = from_origin(x0, y0 + n * res, res, res)
road_px = rasterio.features.rasterize(((g, 1) for g in roads.geometry), out_shape=(n, n),
transform=tf, all_touched=False).astype("float32")
k = int(50 / res)
yy, xx = np.mgrid[-k:k + 1, -k:k + 1]
disc = (np.hypot(xx, yy) * res <= 50).astype("float32")
road50 = fftconvolve(road_px, disc, mode="same") * res # ≈ metres of road within 50 m
cx, cy = np.meshgrid(x0 + (np.arange(n) + 0.5) * res, y0 + n * res - (np.arange(n) + 0.5) * res)
kmc = np.hypot(cx - centre.x, cy - centre.y) / 1000
road50 = np.clip(road50, 0, st.road_50m.max()) # stay inside the fitted range
no2_map = beta[0] + beta[1] * road50 + beta[2] * kmc
print(f"{n} × {n} cells: median {np.median(no2_map):.1f}, 99th percentile "
f"{np.percentile(no2_map, 99):.1f}, max {no2_map.max():.1f} µg/m³")
print(f"cells within 50 m of a major or secondary road: {(road50 > 0).mean():.0%}")2400 × 2400 cells: median 12.9, 99th percentile 27.1, max 31.5 µg/m³
cells within 50 m of a major or secondary road: 58%
The map has ridges along roads rather than discs around monitors, and its highest values sit wherever road density is highest, monitor or not — which is the point, and the reason the predictor range is clipped: beyond what the stations cover, the model has no evidence.
Why small-buffer road length works
Why road length in a small buffer works
NO2 from traffic falls steeply with distance from the kerb, halving within tens of metres. Road length within 50 m is a crude but effective measure of how much traffic emission is close enough to matter. Larger buffers pick up the neighbourhood's general density, which the distance-to-centre term already carries.
Why so few predictors
Every added predictor fits one more direction of the station noise. With seventeen points, a five-predictor model can reach a fitted R² near 1 and a cross-validated one near zero. Keeping two predictors with the right signs and a physical story is what lets the map be trusted away from the stations.
Why OSM road classes are a proxy
The model really wants traffic volume, especially of diesel vehicles. OSM road class is a stand-in: most primary roads are busy and most residential streets are not. Where traffic counts exist, road length weighted by volume is a better predictor.
Temporal scope, transferability and canyons
- Temporal scope: LUR predicts annual means. It says nothing about hours or episodes.
- Transferability: a model fitted in Berlin does not carry over to another city without refitting.
- Street canyons raise concentrations beyond what road length predicts; building height near the road is a common extra predictor.
- Sampling campaigns: many LUR studies deploy 40–80 passive samplers for a few weeks to get enough sites; the regulatory network alone is thin.
- Spatial autocorrelation in residuals suggests a missing predictor; kriging the residuals is the standard hybrid.
Internal links
- Interpolated pollution map shows bullseyes — why interpolation fails on roads
- How to interpolate air quality station data — the background part
- Air quality monitoring networks explained — traffic and background stations
- How to calculate population-weighted exposure — using the prediction map
- Spatial cross-validation in Python — validation beyond leave-one-out
- How to download OpenStreetMap data in Python — the road predictors
FAQ
What is land-use regression?
A regression of measured pollution on the characteristics of each site's surroundings — roads, land use, population — applied to predict concentrations everywhere those characteristics are known.
How many monitoring sites do I need?
At least ten per predictor as a rule of thumb. Regulatory networks are often too thin; dedicated sampling campaigns add sites.
Should traffic stations be included?
Yes. They are what links high concentrations to road predictors.
How do I validate an LUR model?
Leave-one-out or spatially blocked cross-validation, repeating any predictor selection inside each fold.
Which predictors work best for NO2?
Road length or traffic within small buffers (25–100 m) and a measure of the urban gradient, such as distance from the centre or built-up density.
Can I apply an LUR model to another city?
Not reliably without refitting; the coefficients encode local fleet, climate and street geometry.