Interpolating Environmental Exposure Surfaces

This guide is part of the Environmental Exposure & Risk Surface Modeling section. It covers the case where you have direct measurements of the exposure — a network of air quality monitors, well samples, noise loggers — and need a value everywhere, including the places nobody measured.

Concept & Epidemiological Alignment

Interpolation answers a narrow question: given measurements at $n$ locations, what is the most defensible estimate at an unmeasured location? Geostatistical methods answer it by first estimating how quickly the measured quantity decorrelates with distance, then using that model to weight nearby observations. The weighting is not a heuristic — for kriging it is the set of weights that minimises expected squared prediction error under the fitted covariance model, which is why kriging alone returns a prediction variance along with the prediction.

That variance is the reason this family matters to epidemiology. Every other way of producing a surface hands you a number with no stated precision, and a number with no stated precision cannot be carried honestly into a health model. Kriging hands you $\hat{Z}(s_0)$ and $\sigma^2(s_0)$ at every location, and the second quantity is what tells you whether the eastern half of the state is being described or guessed.

The method’s assumptions are strong and must be checked against surveillance reality rather than assumed. Second-order stationarity requires that the covariance between two points depends only on the vector between them, not on where they are — a state where one region is dominated by industrial sources and another by wildfire smoke violates it. Isotropy, the default in most fitting code, requires that decorrelation is the same in all directions, which prevailing winds routinely break. Adequate sampling means enough site pairs at short separations to estimate the variogram near the origin, where it matters most; below about thirty sites, the fitted model is largely an assumption in numerical clothing.

Interpolation is the right tool for regionally smooth exposures — fine particulate matter, ozone, temperature, groundwater chemistry. It is the wrong tool for anything whose variation is driven by a source at street scale, because the method’s premise is that nearby places are similar, and a roadside and the garden behind it are not. For those, Land-Use Regression for Air Pollution Exposure is the correct family.

Method Selection

The decision is mostly made by two facts you can measure before fitting anything: how many sites you have, and whether a plot of the measured values against easting and northing shows a systematic trend.

Spatial Data Prerequisites

• Geometry: monitor sites as points with a stable site identifier; the prediction target as either a regular grid or the centroids or geometries you will summarise to.
• Coordinate reference system: one projected metric system for sites and prediction grid, chosen for the study extent. Variogram distances are in the units of this system, so a geographic system silently produces a variogram in degrees.
• Sample size: at least 30 sites for a defensible variogram, preferably 50 or more; below that, use a published variogram from the literature and say so, or use a non-geostatistical method.
• Temporal alignment: one averaging window applied identically at every site, with a completeness rule and a record of which sites it removed.
• Value distribution: strongly right-skewed measurements — many environmental concentrations are — should be log-transformed before fitting, with back-transformation handled explicitly rather than by exponentiating the mean.
• Duplicate coordinates: two sites at the same location produce a singular kriging system. Deduplicate or jitter deliberately, and log the decision.

Production Implementation

The routine below fits an exponential variogram by weighted least squares, solves the ordinary kriging system at each prediction location using the nearest $k$ sites, and returns prediction and variance together. It keeps the neighbourhood local, which is what makes a state-wide grid tractable, and it logs the fitted parameters so the surface is reproducible.

# Ordinary kriging of a monitor network onto a prediction grid.
# Pinned: geopandas==0.14.4, numpy==1.26.4, scipy==1.13.1, pyproj==3.6.1
import hashlib
import json
import logging

import numpy as np
import geopandas as gpd
from scipy.optimize import curve_fit
from scipy.spatial import cKDTree

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


def exponential_model(h, nugget, sill, rng):
    """Exponential semivariogram. rng is the practical range / 3."""
    return nugget + (sill - nugget) * (1.0 - np.exp(-h / rng))


def fit_variogram(xy, z, n_lags=15, max_frac=0.5):
    """Empirical semivariogram + weighted least-squares exponential fit."""
    d = np.sqrt(((xy[:, None, :] - xy[None, :, :]) ** 2).sum(-1))
    iu = np.triu_indices_from(d, k=1)
    h, g = d[iu], 0.5 * (z[iu[0]] - z[iu[1]]) ** 2
    cutoff = max_frac * h.max()
    keep = h <= cutoff
    edges = np.linspace(0, cutoff, n_lags + 1)
    idx = np.digitize(h[keep], edges) - 1
    lag_h, lag_g, lag_n = [], [], []
    for b in range(n_lags):
        m = idx == b
        if m.sum() < 30:          # too few pairs to trust this lag
            continue
        lag_h.append(h[keep][m].mean())
        lag_g.append(g[keep][m].mean())
        lag_n.append(int(m.sum()))
    lag_h, lag_g = np.asarray(lag_h), np.asarray(lag_g)
    p0 = [lag_g.min(), lag_g.max(), cutoff / 3.0]
    # Weight each lag by its pair count: short lags decide the near-origin shape.
    popt, _ = curve_fit(exponential_model, lag_h, lag_g, p0=p0,
                        sigma=1.0 / np.sqrt(lag_n), maxfev=20000,
                        bounds=([0, 0, 1.0], [np.inf, np.inf, np.inf]))
    log.info("variogram nugget=%.4f sill=%.4f range=%.0f m lags=%d",
             popt[0], popt[1], popt[2] * 3.0, len(lag_h))
    return tuple(float(v) for v in popt)


def ordinary_kriging(sites: gpd.GeoDataFrame, value_col: str,
                     grid: gpd.GeoDataFrame, k: int = 24):
    """Predict value_col at every grid point with a local ordinary kriging system.

    Returns (prediction, variance) arrays aligned to `grid`. Both frames must
    already share one projected CRS in metres.
    """
    if sites.crs is None or grid.crs is None or sites.crs != grid.crs:
        raise ValueError("sites and grid must share one projected CRS")
    if sites.crs.is_geographic:
        raise ValueError("project to metres before kriging; degrees are not distances")

    sites = sites.sort_values("site_id").reset_index(drop=True)  # deterministic
    xy = np.c_[sites.geometry.x.values, sites.geometry.y.values]
    z = sites[value_col].to_numpy(float)
    nugget, sill, rng = fit_variogram(xy, z)

    tree = cKDTree(xy)
    gxy = np.c_[grid.geometry.x.values, grid.geometry.y.values]
    k = min(k, len(sites))
    dists, nbr = tree.query(gxy, k=k)
    pred = np.empty(len(gxy))
    var = np.empty(len(gxy))

    for i in range(len(gxy)):
        idx = nbr[i]
        h = np.sqrt(((xy[idx][:, None, :] - xy[idx][None, :, :]) ** 2).sum(-1))
        A = np.ones((k + 1, k + 1))
        A[:k, :k] = exponential_model(h, nugget, sill, rng)
        np.fill_diagonal(A[:k, :k], 0.0)
        A[k, k] = 0.0
        b = np.ones(k + 1)
        b[:k] = exponential_model(dists[i], nugget, sill, rng)
        try:
            w = np.linalg.solve(A, b)
        except np.linalg.LinAlgError:          # duplicate coordinates
            w = np.linalg.lstsq(A, b, rcond=None)[0]
        pred[i] = w[:k] @ z[idx]
        var[i] = max(w[:k] @ b[:k] + w[k], 0.0)

    cfg = {"model": "exponential", "nugget": nugget, "sill": sill,
           "range_m": rng * 3.0, "k": k, "crs": sites.crs.to_string()}
    log.info("kriged %d locations cfg_sha256=%s", len(gxy),
             hashlib.sha256(json.dumps(cfg, sort_keys=True).encode()).hexdigest()[:16])
    return pred, var

Three details in that code are the difference between a working surface and a plausible fiction. The variogram fit is weighted by pair count, so the short lags — which control the near-origin behaviour and therefore almost all of the prediction weight — dominate the fit instead of being outvoted by the noisy long lags. The kriging system is solved on the nearest $k$ sites rather than globally, which keeps the cost linear in grid size and also limits the damage from a distant outlier. And the variance is floored at zero, because a numerically negative kriging variance signals a non-positive-definite system rather than a physically meaningful result.

Parameter Selection & Tuning

The variogram model family matters less than the near-origin behaviour. Exponential and spherical models fit most environmental data adequately; a Gaussian model implies an unrealistically smooth field and frequently produces an unstable kriging system. Fit two families and compare their cross-validated error rather than arguing about which is theoretically preferable.

The lag binning is the parameter most people leave at a default and should not. Bins must be narrow enough to resolve the rise near the origin and wide enough that each holds enough pairs — thirty is a common minimum — for a stable mean. The maximum lag should stop at roughly half the study extent, because beyond that the pair count collapses and the empirical points are noise.

The neighbourhood size $k$ trades bias against cost. Values between 16 and 32 are typical; below about 12 the prediction becomes jumpy at cell boundaries, and above 40 the gain is negligible while the solve cost grows cubically. Check the choice by re-running cross-validation at two settings — if the error is insensitive, take the cheaper one.

Anisotropy should be tested, not assumed away. Compute directional variograms in four sectors; if the range differs by more than roughly a factor of two between the longest and shortest directions, an isotropic model is misallocating weight, and the surface will be smeared perpendicular to the true gradient.

The nugget should not be forced to zero to make the surface pass exactly through the observations. An exact interpolator is only appropriate when the measurements have no error, which is never true of environmental monitoring. Forcing a zero nugget produces the characteristic bullseye artefact around each site.

Edge Cases & Failure Modes

Fewer than thirty sites. The fitted variogram is dominated by sampling noise and the resulting surface has an unearned appearance of precision. Use a published variogram for the pollutant and region, declare that you did, and inflate the reported uncertainty accordingly.

A gap wider than the range. Inside such a gap, kriging returns the global mean with a variance near the sill — mathematically honest, cartographically misleading, because the colour ramp shows a smooth field. Mask or flag every prediction whose variance exceeds a stated fraction of the sill rather than shipping it as if it were informed.

Skewed data and naive back-transformation. Kriging the log of a concentration and exponentiating the prediction returns a median, not a mean, and understates the mean by a factor that depends on the variance. If your health model needs a mean, apply the lognormal correction explicitly and note it.

Preferential monitor siting. Regulatory networks are sited where problems are expected, which means the sample is not representative of the region. Kriging assumes nothing about representativeness for the prediction itself, but it does mean that the study-area average computed from the surface is biased. Report area-weighted summaries from the grid, not the mean of the monitors.

Temporal gaps masquerading as spatial structure. A site that reported for eight months of a year, kept by a lenient completeness rule, contributes a value from a different season than its neighbours. This appears as a spurious short-range nugget. Enforce the completeness rule before fitting, not after.

Memory on a fine grid. A 30 m grid over a state is hundreds of millions of cells. Predict in tiles with a halo wider than the range, write each tile to a cloud-optimized GeoTIFF, and never assemble the whole array in memory — the patterns for this are covered in Cloud-Native Geospatial Storage for Surveillance.

The variance band is where most of those failure modes become visible at once. A map of prediction variance beside the prediction itself is the single most useful diagnostic an interpolated surface can ship with:

Compliance & Audit Controls

A kriged surface is a derived product with a long downstream life, so its provenance must travel with it. Record the source monitor dataset and vintage, the averaging window and completeness rule, the site count after the rule, the variogram family and fitted nugget, sill and range, the neighbourhood size, the projected coordinate reference system, and the cross-validated error and bias. Hash that configuration and write the digest into the raster metadata so a given file maps to exactly one run.

Determinism requires a stable site sort before fitting and a fixed optimiser seed where the fit is stochastic. The surface itself is not protected data — it contains no personal information — but the moment it is joined to residences it enters the protected environment, and the join is a separately logged step. Publish the surface and its variance band openly where policy allows: an exposure surface that external analysts can re-run is the cheapest form of peer review available.

Production Implementation Checklist

• Sites and prediction grid share one projected metric coordinate reference system
• Completeness rule applied before fitting and the number of dropped sites recorded
• Site count confirmed adequate for variogram estimation, or a published variogram used and declared
• Values checked for skew and log-transformed with an explicit back-transformation where needed
• Duplicate site coordinates resolved before the kriging system is built
• Empirical variogram binned with a minimum pair count per lag and a cutoff near half the extent
• Directional variograms inspected and anisotropy modelled if the ranges differ materially
• Nugget fitted rather than forced to zero
• Neighbourhood size chosen by cross-validation rather than left at a default
• Leave-one-site-out validation reported as error, bias and fit in exposure units
• Predictions with variance above the agreed threshold masked or flagged, not published as informed
• Prediction variance written as a band alongside the prediction, not discarded
• Variogram parameters and configuration hash recorded in the output metadata

Frequently Asked Questions

Is inverse distance weighting ever acceptable? Yes, for screening and for networks too small to fit a variogram, provided you do not attach an error bar to it. Its weakness is that the power parameter is chosen rather than estimated, so the smoothness of the surface is an author’s decision instead of a property of the data, and it gives you no way to say where the surface is uninformed.

How do I krige a whole time series without refitting every day? Fit the variogram on a pooled or seasonal basis rather than daily, because a single day’s network rarely has enough pairs. Check the pooled model against a handful of individual days, and refit when the cross-validated error drifts rather than on a fixed schedule.

Why does my kriged surface have rings around the monitors? Almost always a nugget forced to zero, which makes kriging an exact interpolator and forces the surface to pass through each measurement. Fit the nugget instead. A secondary cause is a neighbourhood so small that the set of contributing sites changes abruptly between adjacent cells.

Should I krige the pollutant or its logarithm? Krige whichever is closer to symmetric, which for most concentration data is the logarithm. Then decide deliberately what the health model needs: exponentiating gives a median, and recovering a mean requires the lognormal correction using the local kriging variance.

What do I do about the prediction variance in the health model? At minimum, report it and restrict sensitivity analyses to well-informed areas. Better, carry it forward — either by multiple imputation from the predictive distribution or by a measurement-error model. The methods are set out in Exposure Measurement Error & Uncertainty.

Related Topics

• Environmental Exposure & Risk Surface Modeling — the parent section and where interpolation sits among the four modeling families.
• Kriging PM2.5 Monitor Data for Tract Exposure — the end-to-end run from a state monitor pull to a tract-level exposure table.
• Cross-Validating an Exposure Surface With Held-Out Monitors — the validation protocol that reports the error which actually applies.
• Exposure Measurement Error & Uncertainty — carrying kriging variance into the health effect estimate instead of discarding it.