Spatial Epi Public Health GIS

Building Buffer Predictors for a Land-Use Regression Model

This guide, part of Land-Use Regression for Air Pollution Exposure, covers the step that consumes most of the wall-clock time and produces most of the silent errors: turning a road network and a land-cover layer into a rectangular table of predictors, first at sixty monitoring sites and then at two million prediction points.

Problem Context & Constraints

The naive implementation buffers each site, clips the road network to the buffer, and sums the segment lengths. It is correct, it is obvious, and it is unusable at prediction scale: with six buffer radii, two road classes and two land-cover classes, a metropolitan grid of two million points requires forty-eight million clip operations. A loop that takes 40 ms per site finishes the campaign in three seconds and the grid in three weeks.

The correctness traps are worse than the performance one because they are invisible. Clipping a line at a buffer edge is not the same as selecting lines that intersect the buffer. The latter counts a 5 km motorway segment in full because one end clips the 100 m buffer, inflating the predictor by an order of magnitude at exactly the sites where near-road exposure matters. A buffer built in a geographic coordinate system is an ellipse of varying size, so a 300 m buffer is 300 m of latitude and something else entirely of longitude. Land-cover polygons that overlap — common where a source layer has been edited over time — double-count area, pushing a proportion above one.

And one design constraint runs through the whole build: the predictor computed at a monitoring site and the predictor computed at a prediction point must come from the same code path. If sites are processed with an exact overlay and the grid with a raster approximation, the model is fitted on one definition and applied to another, and the mismatch appears as unexplained bias in exactly the places the campaign did not sample.

Clip the Segment, Do Not Select It Two panels showing the same monitoring site with a one hundred metre buffer and a long motorway segment passing nearby. In the left panel the road is clipped at the buffer boundary and only the portion inside is measured, giving a road-length predictor of about one hundred and eighty metres. In the right panel the whole segment is selected because it intersects the buffer, giving a predictor of five thousand two hundred metres for the same location. The caption notes the error is largest at the small buffers, which are precisely the predictors that carry the near-road signal. Same site, same road, two predictor values differing by 30 times Correct: clip at the boundary road_major_100 = 182 m Wrong: select whole segment road_major_100 = 5 210 m The error grows as the buffer shrinks, so it is worst in the predictors that carry the near-road gradient and it survives every downstream check, because the table is complete and the units look right

Prerequisites

  • python 3.11
  • geopandas 0.14.4, shapely 2.0.4, numpy 1.26.4, pandas 2.2.2
  • rasterio 1.3.10 and rasterstats 0.19.0 for the raster path
  • A road network with a functional class attribute and no duplicate geometries
  • A land-cover layer with non-overlapping polygons, or a categorical raster
  • Every layer in one projected coordinate reference system in metres, selected per Coordinate Reference Systems for Public Health

Step-by-Step Solution

The build below uses a spatial index to restrict the candidate geometries per point, clips rather than selects, and processes the largest buffer first so smaller buffers reuse the already-restricted candidate set.

# Buffer predictors for land-use regression, vectorised and deterministic.
# Pinned: geopandas==0.14.4, shapely==2.0.4, numpy==1.26.4, pandas==2.2.2
import logging

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.strtree import STRtree

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("lur.buffers")

BUFFERS_M = (25, 50, 100, 300, 1000, 3000)


def road_length_predictors(points: gpd.GeoDataFrame, roads: gpd.GeoDataFrame,
                           class_col: str = "fclass",
                           classes=("motorway", "trunk", "primary"),
                           buffers=BUFFERS_M) -> pd.DataFrame:
    """Clipped road length within each buffer, for every point.

    Uses one STRtree over the road geometries and queries the LARGEST buffer once
    per point, then clips the resulting candidate set against every smaller buffer.
    """
    _require_metric(points, roads)
    points = points.sort_values("point_id").reset_index(drop=True)
    sel = roads[roads[class_col].isin(classes)].reset_index(drop=True)
    geoms = sel.geometry.to_numpy()
    tree = STRtree(geoms)
    rmax = max(buffers)

    out = {f"road_major_{b}": np.zeros(len(points)) for b in buffers}
    for i, pt in enumerate(points.geometry):
        # One index query at the widest radius; everything else is arithmetic.
        cand_idx = tree.query(pt.buffer(rmax))
        if len(cand_idx) == 0:
            continue
        cand = geoms[cand_idx]
        for b in buffers:
            ring = pt.buffer(b)
            total = 0.0
            for g in cand:
                if not g.intersects(ring):
                    continue
                total += g.intersection(ring).length   # CLIP, never select whole
            out[f"road_major_{b}"][i] = total
        if i and i % 50_000 == 0:
            log.info("road predictors: %d/%d points", i, len(points))

    df = pd.DataFrame(out)
    df.insert(0, "point_id", points["point_id"].to_numpy())
    log.info("road predictors complete: %d points x %d buffers", len(df), len(buffers))
    return df


def landcover_area_predictors(points: gpd.GeoDataFrame, landcover: gpd.GeoDataFrame,
                              class_col: str = "class",
                              classes=("industrial", "forest"),
                              buffers=BUFFERS_M) -> pd.DataFrame:
    """Area proportion of each land-cover class within each buffer."""
    _require_metric(points, landcover)
    points = points.sort_values("point_id").reset_index(drop=True)
    out = {}
    for cls in classes:
        lc = landcover[landcover[class_col] == cls]
        # Dissolve first: overlapping source polygons would double-count area.
        merged = lc.geometry.union_all()
        tree = STRtree(np.asarray([merged]))
        for b in buffers:
            col = np.zeros(len(points))
            area_b = np.pi * b * b
            for i, pt in enumerate(points.geometry):
                ring = pt.buffer(b)
                if len(tree.query(ring)) == 0:
                    continue
                col[i] = merged.intersection(ring).area / area_b
            out[f"{'green' if cls == 'forest' else cls}_{b}"] = col
    df = pd.DataFrame(out)
    df.insert(0, "point_id", points["point_id"].to_numpy())
    return df


def _require_metric(*frames: gpd.GeoDataFrame) -> None:
    crs = {f.crs for f in frames}
    if len(crs) != 1 or next(iter(crs)) is None:
        raise ValueError("all layers must share one CRS")
    if next(iter(crs)).is_geographic:
        raise ValueError("buffers must be built in a projected CRS measured in metres")

Three implementation choices carry the correctness of this build. The intersection(...).length call is the clip, and replacing it with a length of the whole candidate geometry is the error shown in the diagram above. The union_all() dissolve before the land-cover intersection removes double counting from overlapping source polygons, which otherwise pushes a proportion above 1.0 and produces a predictor the regression cannot interpret. And the coordinate reference system guard fails loudly rather than silently producing degree-based buffers.

For a prediction grid of millions of points, the pure-Python inner loop is still too slow, and the right move is to switch the road predictors to a rasterised path: burn road length per cell once at a fine resolution, then compute a moving-window sum at each buffer radius with an integral image. The result is numerically close but not identical to the exact clip, so the sites must be recomputed the same way — fit on the raster definition if you predict on the raster definition.

Two Build Paths, One Rule Two computation paths for buffer predictors. The vector path clips geometries exactly at each buffer boundary and is accurate but scales poorly, suiting the tens of monitoring sites. The raster path burns road length into fine cells once and computes moving-window sums with an integral image, scaling to millions of prediction points but approximating the buffer edge. A rule beneath both states that whichever path produces the prediction grid must also produce the site predictors, because fitting on one definition and predicting on the other introduces bias that appears only where the campaign did not sample. Pick a path for accuracy or for scale — then use it for both tables Vector clip exact length inside the buffer one index query per point milliseconds per point fits tens to thousands of sites too slow for a full grid Raster moving window road length burned to fine cells integral image per radius microseconds per point scales to millions of points approximates the buffer edge The rule: the model is fitted on whichever definition the grid will use mixing them puts bias exactly where the campaign has no site to reveal it Record the path, the cell size and the buffer list in the model provenance

Validation & Edge Cases

Monotonicity across the ladder. Road length within 1000 m must be at least road length within 300 m at the same point. A violation means the clip is wrong or the candidate set was restricted per buffer instead of at the widest radius. Assert it.

Proportions inside the unit interval. Land-cover area proportion must lie between 0 and 1. Values above 1 mean overlapping source polygons survived the dissolve; values that sum above 1 across mutually exclusive classes mean the same.

Points near the study-area edge. A buffer at the boundary extends into territory with no road data, so the predictor is an undercount rather than a true low value. Pad the input layers by the widest buffer radius beyond the prediction extent, or flag the affected points.

Empty candidate sets. A point with no road within 3 km legitimately gets zero, but a systematic block of zeros usually means a coordinate reference system mismatch that the guard did not catch because both layers were wrong in the same way. Spot-check a handful of known locations against a map.

INFO lur.buffers road predictors: 500000/2118400 points
INFO lur.buffers road predictors complete: 2118400 points x 6 buffers
ERROR lur.buffers monotonicity violated at 14 points: road_major_300 > road_major_1000

Duplicate road geometries. Networks assembled from several sources routinely contain the same carriageway twice. Every duplicate doubles its contribution to every buffer containing it. Deduplicate on geometry before building anything.

Sites and grid built on different vintages. If the campaign predictors were built in March against one road extract and the grid in June against another, a new bypass appears in one table and not the other. Pin the extract, and record its identifier.

The two invariants in that list are cheap to assert and catch most build errors, because a correct predictor table has a shape you can predict in advance:

A Correct Ladder Never Goes Down Road length plotted against buffer radius for one monitoring site. In a correct build the value rises monotonically, because each larger buffer contains every smaller one. In a faulty build where the candidate set was restricted separately for each buffer, the value at three hundred metres exceeds the value at one kilometre, which is geometrically impossible and exposes the bug immediately. The monotonicity assertion catches this at build time without needing any expected values. Bigger buffers contain smaller ones, so the ladder can only rise correct build 1 km below 300 m: impossible 25 50 100 300 1000 3000 buffer radius in metres, one monitoring site road length

Compliance Notes

The predictor build is a deterministic transformation and must be recorded as one. Store the road and land-cover source identifiers with their extract dates, the functional classes included in each aggregate, the exact buffer list, the computation path (vector clip or raster window) with the cell size where applicable, the dissolve step for overlapping polygons, and the padding applied at the study-area edge.

Record the assertion results too, not just the table. A build that passed the monotonicity and unit-interval checks is a different artifact from one where the checks were never run, and the difference is invisible in the output file. Where the same predictor table is reused across studies, the checksum of the table plus the configuration digest is what lets a second team confirm they are working from the same inputs rather than a rebuild that quietly used a newer network.