Spatial Epi Public Health GIS

Kriging PM2.5 Monitor Data for Tract Exposure

This guide, part of Interpolating Environmental Exposure Surfaces, walks the specific operational job of turning a year of state fine-particulate monitor readings into one exposure value and one uncertainty value per census tract — the form a small-area health analysis actually consumes.

Problem Context & Constraints

The naive version of this task is a three-line script: read the monitors, krige onto the tract centroids, join. It produces a table that looks correct and is wrong in three specific ways that a reviewer will find.

A tract is an area, not a point. Assigning the value at the centroid treats a 40 km² rural tract as though everyone in it lived at one spot, and in tracts with a strong internal gradient — anything containing a highway corridor or a city edge — the centroid value can miss the population-weighted mean by more than the exposure contrast the study is trying to detect. The tract value must be an average over the tract, and ideally an average weighted by where people actually live.

Annual means are not comparable unless completeness is enforced. A monitor that ran from March to October reports a warm-season mean, and fine particulate matter in most of North America peaks in winter. Blending that site with full-year sites shifts the whole surface downward in its neighbourhood, and the shift looks like real spatial structure.

The variance is not optional here. Tract-level exposure feeds directly into rate models, and in the parts of a state with no monitor within the fitted range, the kriged value is essentially the regional mean. If the exposure table does not carry a per-tract variance, the health model cannot distinguish a measured tract from a guessed one, and the guessed ones are systematically rural.

Centroid Value Versus Population-Weighted Tract Mean A single census tract crossed by a concentration gradient that runs from high near a highway corridor on one side to low on the far side. The geometric centroid sits in the low-density middle of the tract and takes a mid-range value. The residential population, drawn as a band of dwellings, is concentrated near the corridor where the concentration is high. The population-weighted mean is therefore substantially higher than the centroid value, and the difference is labelled as the error introduced by centroid assignment. The centroid is rarely where the people are one census tract high low concentration gradient across the tract dwellings centroid Centroid value mid-gradient, low population there Population-weighted mean pulled toward the dwelling band the gap is the assignment error, and it is not random

Prerequisites

  • python 3.11
  • geopandas 0.14.4 (pulls shapely 2.0.x, pyproj 3.6.x)
  • numpy 1.26.4, scipy 1.13.1
  • rasterio 1.3.10 for writing the surface

Input state assumed by the code below:

  • Daily monitor readings for one calendar year as a table with a site identifier, date, concentration and site coordinates.
  • Tract polygons for the same vintage as the denominators the study uses, in a projected coordinate reference system in metres. Tract and monitor frames must match; the selection logic is in Coordinate Reference Systems for Public Health.
  • A population raster or block-level population for weighting the within-tract average. Without it the code falls back to an unweighted areal mean and says so.
  • A fitted kriging routine — the ordinary kriging implementation from the parent guide is assumed available as ordinary_kriging.

Step-by-Step Solution

The run has four stages: reduce daily readings to site-year means under a completeness rule, krige onto a regular grid covering the state, aggregate grid cells to tracts with population weights, and write the table with its variance column.

# Annual PM2.5 monitor readings -> population-weighted tract exposure with variance.
# Pinned: geopandas==0.14.4, numpy==1.26.4, scipy==1.13.1, rasterio==1.3.10
import logging

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

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

MIN_DAYS = 274          # 75% of a 365-day year, the usual completeness convention
GRID_M = 1000           # 1 km prediction grid


def site_year_means(daily: pd.DataFrame, sites: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Collapse daily readings to one mean per site, dropping incomplete records."""
    counts = daily.groupby("site_id")["pm25"].agg(["count", "mean"])
    keep = counts[counts["count"] >= MIN_DAYS]
    dropped = sorted(set(counts.index) - set(keep.index))
    log.info("sites kept=%d dropped_for_completeness=%d ids=%s",
             len(keep), len(dropped), dropped[:10])
    out = sites.merge(keep["mean"].rename("pm25"), left_on="site_id", right_index=True)
    return out.sort_values("site_id").reset_index(drop=True)


def build_grid(extent: gpd.GeoDataFrame, step: int = GRID_M) -> gpd.GeoDataFrame:
    """Regular point grid over the study extent, in the extent's CRS."""
    minx, miny, maxx, maxy = extent.total_bounds
    xs = np.arange(minx + step / 2, maxx, step)
    ys = np.arange(miny + step / 2, maxy, step)
    xx, yy = np.meshgrid(xs, ys)
    pts = [Point(x, y) for x, y in zip(xx.ravel(), yy.ravel())]
    grid = gpd.GeoDataFrame({"cell_id": np.arange(len(pts))},
                            geometry=pts, crs=extent.crs)
    log.info("grid cells=%d step=%dm", len(grid), step)
    return grid


def tract_exposure(grid: gpd.GeoDataFrame, pred: np.ndarray, var: np.ndarray,
                   tracts: gpd.GeoDataFrame,
                   pop_col: str | None = "cell_pop") -> gpd.GeoDataFrame:
    """Aggregate kriged grid cells to tracts, weighting by population where available."""
    g = grid.copy()
    g["pm25"], g["pm25_var"] = pred, var
    joined = gpd.sjoin(g, tracts[["geoid", "geometry"]], how="inner", predicate="within")

    if pop_col and pop_col in joined.columns and joined[pop_col].sum() > 0:
        joined["w"] = joined[pop_col].clip(lower=0)
        weighting = "population"
    else:
        joined["w"] = 1.0
        weighting = "area"
        log.warning("no population weights available — falling back to areal mean")

    def agg(df):
        w = df["w"].to_numpy(float)
        if w.sum() == 0:
            w = np.ones_like(w)
        w = w / w.sum()
        return pd.Series({
            "pm25": float(w @ df["pm25"].to_numpy(float)),
            # Variance of a weighted mean of correlated cells is bounded below by
            # the weighted mean of cell variances; report that bound, not zero.
            "pm25_var": float(w @ df["pm25_var"].to_numpy(float)),
            "n_cells": int(len(df)),
        })

    out = joined.groupby("geoid").apply(agg, include_groups=False).reset_index()
    out["weighting"] = weighting
    out["pm25_se"] = np.sqrt(out["pm25_var"])
    missing = set(tracts["geoid"]) - set(out["geoid"])
    if missing:
        log.warning("tracts with no grid cell inside: %d (smaller than one cell)",
                    len(missing))
    return tracts.merge(out, on="geoid", how="left")

Two choices in that code deserve defending. The tract value is a weighted mean of grid cells rather than a centroid lookup, which is what makes it an areal quantity. And the reported tract variance is the weighted mean of the cell kriging variances — a conservative figure, because neighbouring cells are strongly correlated and the true variance of their mean is closer to the cell variance than to the variance of an independent average. Reporting the smaller, independence-based number would understate uncertainty exactly where the network is thin.

Tracts smaller than one grid cell fall through the spatial join entirely. Dense urban tracts are frequently under 1 km², so this is not an edge case but a routine outcome, and it must be handled explicitly rather than left as a silent null.

The Four Stages and What Each One Drops A left-to-right pipeline in four stages. Stage one holds daily readings from all monitor sites. Stage two applies the seventy-five per cent completeness rule and reports how many sites were removed. Stage three fits the variogram and krige onto a one-kilometre grid, producing a prediction and a variance for every cell. Stage four aggregates cells to census tracts with population weights and flags tracts smaller than one cell. Beneath each stage a note records the quantity that must be logged at that point: the raw record count, the dropped site identifiers, the variogram parameters, and the count of tracts that received no cell. Every stage discards something — log what, or the table is unexplainable 1 · Daily readings all sites, all days one calendar year 2 · Completeness at least 274 days per site or dropped 3 · Krige to grid 1 km cells, prediction and variance per cell 4 · Tract rollup population-weighted mean per tract What the audit log must carry from each stage record count · dropped site ids · nugget, sill, range · tracts with no cell inside A tract with a null exposure is a finding, not a gap to fill with the state mean and a tract smaller than one grid cell needs a point lookup, not a join the same four stages re-run identically from the logged configuration

Validation & Edge Cases

Validate before the table leaves the script, not after it has been merged into a health analysis.

Leave-one-site-out error, in micrograms per cubic metre. Rebuild the surface without each site and predict at that site. For a typical state network and an annual mean, a root mean squared error in the region of one to two micrograms is ordinary; an error approaching the between-site standard deviation means the surface carries no spatial information at all and should not be used.

def loo_error(sites, value_col, krige_fn):
    """Leave-one-site-out: rebuild without each site and predict at it."""
    resid = []
    for i in range(len(sites)):
        train = sites.drop(sites.index[i])
        target = sites.iloc[[i]][["geometry"]].copy()
        pred, _ = krige_fn(train, value_col, target)
        resid.append(float(sites.iloc[i][value_col] - pred[0]))
    resid = np.asarray(resid)
    return {"rmse": float(np.sqrt((resid ** 2).mean())),
            "bias": float(resid.mean()),
            "n": len(resid)}

Tracts that received no grid cell. Predict directly at those tract centroids instead of leaving them null, and mark the row so a reviewer can see which values came from a point lookup.

Tracts whose variance is near the sill. These sit in a network gap; the exposure value is close to the regional mean regardless of what the map shows. Flag them and consider excluding them from any dose-response analysis rather than letting them dilute the estimate.

A completeness rule that removes a whole region. If the dropped sites cluster geographically — a common outcome when one agency’s sub-network had an outage — the surviving network has a new hole that the variogram was not fitted to. Plot the dropped sites before accepting the fit.

INFO pm25.tracts sites kept=41 dropped_for_completeness=7 ids=['060010011', '060130002', ...]
INFO exposure.kriging variogram nugget=0.7412 sill=4.1108 range=148000 m lags=11
INFO pm25.tracts grid cells=428913 step=1000m
WARNING pm25.tracts tracts with no grid cell inside: 236 (smaller than one cell)

Negative predictions. Kriging is not constrained to be positive and can return a small negative concentration in a low corner of the state. Do not clip silently — count them, and if there are more than a handful, refit on the log scale.

Tracts smaller than the grid cell deserve their own picture, because they are concentrated in exactly the dense urban areas where most cases live:

Where Small Tracts Fall Through the Grid A one-kilometre prediction grid drawn as a lattice of cell centre points over part of a city. Large suburban tracts each contain several cell centres and receive a proper weighted mean. Several small downtown tracts lie entirely between cell centres and contain none, so a within-join assigns them nothing. The fallback for those tracts is a point lookup at the tract's population-weighted centre, flagged in the output so a reviewer can see which values came from the fallback. Dense downtown tracts are routinely smaller than one grid cell suburban tract 4 tracts, 0 cell centres Join finds no cell silent null if not handled Point lookup fallback flagged in the output row Count these per run; a sudden rise means the grid resolution changed

Compliance Notes

The exposure table produced here contains no personal data and should be published alongside the analysis where policy allows; the disclosure question arises only when it is joined to case records. Three artifacts make the table defensible on its own: the monitor dataset identifier and vintage together with the completeness rule and the list of dropped sites; the fitted variogram parameters, grid resolution and neighbourhood size; and the leave-one-site-out error and bias.

Two columns must survive into any downstream use. The variance column is what lets a health model distinguish a well-measured urban tract from a guessed rural one. The weighting flag records whether the tract value is population-weighted or a plain areal mean, and the two are not interchangeable — an areal mean in a tract with a strong internal gradient is a different quantity, and a reviewer comparing your result with another study needs to know which one you produced.