Building an IDW Air Quality Surface from Point Sensors in Python #

Inverse distance weighting turns a handful of scattered PM2.5 sensors into a continuous surface by estimating each grid cell as a distance-weighted average of the nearest sensors: closer sensors count more, farther ones less, with the falloff set by a power exponent. It is the deterministic, fast, no-fitting-required member of the interpolation toolkit — the right first move for low-cost air-quality networks where sensors are too few or too clustered to fit a reliable variogram. This page gives a self-contained IDW implementation built on a scipy.spatial.cKDTree, with the grid construction, power and neighbour tuning, and sparse-area masking that separate a usable surface from a field of bullseyes. It sits under the broader spatial interpolation with kriging and IDW for sensor fields workflow.

How Inverse Distance Weighting Behaves #

IDW makes one assumption: the value at an un-sampled point is more like nearby sensors than distant ones. Each estimate is a weighted mean where the weight of sensor i is 1 / d_i**power and d_i is the distance from the target cell to that sensor. Two knobs control everything.

The power exponent sets how fast influence decays with distance. At power 1, influence falls off linearly and the surface is smooth and diffuse — distant sensors still pull noticeably. At power 2 (the inverse-square default) influence drops quickly and local sensors dominate. Push past power 3 and each cell is effectively its nearest sensor, producing hard concentric rings — bullseyes — around isolated nodes. For urban PM2.5, where concentrations spike sharply near roads and fall off within hundreds of metres, a power of 2 to 3 preserves those gradients; for a slowly varying regional haze, a lower power avoids fabricating structure.

The neighbour count k bounds how many sensors vote for each cell. A small k keeps estimates local but makes the surface sensitive to any single bad sensor; a large k smooths across the whole network and washes out genuine hotspots. Because low-cost air-quality networks are uneven — dense in city centres, sparse at the edges — a fixed k is more robust than a fixed search radius, which would leave rural cells with no neighbours at all. The gap is that a fixed k will happily pull in a sensor 15 km away to fill a rural cell, which is why the distance mask in the implementation below matters as much as the weighting itself.

Unlike ordinary kriging, IDW is an exact interpolator that reproduces sensor values at their own locations and never predicts outside the observed range — but it gives no uncertainty estimate, which is the central trade-off explored in ordinary kriging vs IDW for sparse sensor networks.

Production-Ready Implementation #

The function below is self-contained: give it projected sensor coordinates and values, a grid, and tuning parameters, and it returns a masked 2-D surface. It uses a cKDTree for the neighbour search, handles the coincident-sensor divide-by-zero case, and blanks cells beyond the support radius.

# python 3.11 · numpy==1.26.4 · scipy==1.12.0
from __future__ import annotations
import numpy as np
from scipy.spatial import cKDTree


def idw_air_quality_surface(
    sensor_xy: np.ndarray,        # shape (n, 2), projected metres
    sensor_values: np.ndarray,    # shape (n,), e.g. PM2.5 in ug/m3
    grid_x: np.ndarray,           # 2-D meshgrid of eastings
    grid_y: np.ndarray,           # 2-D meshgrid of northings
    power: float = 2.5,
    k: int = 8,
    support_radius_m: float | None = None,
    eps: float = 1e-9,
) -> np.ndarray:
    """
    Interpolate a continuous PM2.5 surface from point sensors with IDW.

    Parameters
    ----------
    sensor_xy : np.ndarray
        (n, 2) array of projected sensor coordinates in metres. Reproject
        from WGS84 to the local UTM zone before calling — IDW weights are
        distance-based and meaningless in degrees.
    sensor_values : np.ndarray
        (n,) measurand values aligned with ``sensor_xy``.
    grid_x, grid_y : np.ndarray
        Matching 2-D coordinate arrays from ``np.meshgrid``.
    power : float
        Distance-decay exponent. 2.0 is the inverse-square default; raise
        toward 3 for sharp roadside gradients, lower toward 1.5 to smooth.
    k : int
        Number of nearest sensors that contribute to each cell.
    support_radius_m : float, optional
        Cells whose nearest sensor is farther than this are set to NaN.
        Defaults to 3x the median nearest-neighbour spacing.
    eps : float
        Small constant preventing division by zero when a cell coincides
        with a sensor location.

    Returns
    -------
    np.ndarray
        2-D surface matching ``grid_x.shape``; unsupported cells are NaN.
    """
    if sensor_xy.shape[0] != sensor_values.shape[0]:
        raise ValueError("sensor_xy and sensor_values length mismatch")

    tree = cKDTree(sensor_xy)
    k_eff = min(k, sensor_xy.shape[0])

    # Auto support radius from the network's own density
    if support_radius_m is None:
        nn_dist, _ = tree.query(sensor_xy, k=2)   # column 1 = nearest other sensor
        support_radius_m = 3.0 * float(np.median(nn_dist[:, 1]))

    targets = np.column_stack([grid_x.ravel(), grid_y.ravel()])
    dist, idx = tree.query(targets, k=k_eff)

    # Ensure 2-D even when k_eff == 1
    if k_eff == 1:
        dist = dist[:, None]
        idx = idx[:, None]

    weights = 1.0 / np.power(dist + eps, power)
    surface = np.sum(weights * sensor_values[idx], axis=1) / np.sum(weights, axis=1)

    # Mask cells with no genuine support
    nearest = dist[:, 0]
    surface[nearest > support_radius_m] = np.nan

    return surface.reshape(grid_x.shape)


def build_meshgrid(
    sensor_xy: np.ndarray,
    cell_size_m: float = 200.0,
    pad_m: float = 500.0,
) -> tuple[np.ndarray, np.ndarray]:
    """Regular grid over the sensor extent plus a padding margin (metres)."""
    xmin, ymin = sensor_xy[:, 0].min() - pad_m, sensor_xy[:, 1].min() - pad_m
    xmax, ymax = sensor_xy[:, 0].max() + pad_m, sensor_xy[:, 1].max() + pad_m
    xs = np.arange(xmin, xmax + cell_size_m, cell_size_m)
    ys = np.arange(ymin, ymax + cell_size_m, cell_size_m)
    return np.meshgrid(xs, ys)

The KD-tree query costs O(log n) per cell, so the whole surface is O(g · k · log n) for g grid cells and n sensors — a 200×200 grid over 50 sensors completes in well under a second. Memory is O(g + n): only the grid and the sensor arrays are held, never an n×n matrix, which is why IDW scales to large grids where dense kriging cannot.

Parameter Tuning Guide #

Power and k interact with sensor spacing and the sharpness of the pollutant field. These starting points come from low-cost urban and regional deployments; confirm with leave-one-out cross-validation before publishing a surface.

Grid resolution: keep cell_size_m at or above the median sensor spacing. A 50 m grid over sensors 2 km apart looks detailed but every sub-kilometre feature is interpolation fiction. Auto support radius: the 3× median nearest-neighbour default adapts to density — a dense downtown cluster gets a tight radius, a sparse rural fringe a wider one, but truly isolated cells still get masked.

Verification and Testing #

Two properties make IDW easy to test: it reproduces sensor values exactly at sensor locations, and it never leaves the observed range. The test below asserts both, plus that a sensor far outside the network gets masked.

import numpy as np
import pytest


def test_idw_reproduces_values_and_respects_range():
    rng = np.random.default_rng(42)
    sensor_xy = rng.uniform(0, 10_000, size=(25, 2))       # metres
    sensor_values = rng.uniform(5.0, 60.0, size=25)         # PM2.5 ug/m3

    grid_x, grid_y = build_meshgrid(sensor_xy, cell_size_m=200.0)
    surface = idw_air_quality_surface(
        sensor_xy, sensor_values, grid_x, grid_y, power=2.0, k=8
    )

    # 1. Exact interpolation: value at a sensor location equals its reading
    sx, sy = sensor_xy[0]
    single = idw_air_quality_surface(
        sensor_xy, sensor_values,
        np.array([[sx]]), np.array([[sy]]),
        power=2.0, k=8,
    )
    assert single[0, 0] == pytest.approx(sensor_values[0], abs=1e-3)

    # 2. Range preservation: no cell exceeds the observed min/max
    finite = surface[np.isfinite(surface)]
    assert finite.min() >= sensor_values.min() - 1e-6
    assert finite.max() <= sensor_values.max() + 1e-6


def test_idw_masks_unsupported_cells():
    sensor_xy = np.array([[0.0, 0.0], [500.0, 0.0], [0.0, 500.0]])
    sensor_values = np.array([10.0, 20.0, 30.0])
    # A grid cell 50 km away must be masked
    gx = np.array([[50_000.0]])
    gy = np.array([[50_000.0]])
    out = idw_air_quality_surface(sensor_xy, sensor_values, gx, gy)
    assert np.isnan(out[0, 0])

The exact-reproduction check is the single most useful regression test for an IDW pipeline: if a refactor breaks the eps handling or the weight normalisation, this assertion fails immediately. For a field comparison, cross-validate against a co-located reference-grade monitor and confirm the residual RMSE is within the tolerance your reporting requires.

Gotchas #

Coincident grid cell and sensor. When a grid cell falls exactly on a sensor, the distance is zero and the naive weight 1/0 is infinite. The eps term handles it, but confirm your eps is small relative to your cell size — too large and it dampens genuine near-sensor weighting.

A fixed k reaches too far in sparse areas. Nearest-k guarantees every cell an estimate, but in a rural corner it will pull in sensors tens of kilometres away and produce a confident-looking value with no real support. The distance mask is not optional — always blank cells beyond the support radius.

Outliers become bullseyes. IDW has no outlier defence. A single sensor stuck at a calibration-reset maximum paints a ring of inflated values across every cell it influences. Clean and range-check sensor values before interpolating, not after.

Interpolating in degrees. Passing raw longitude/latitude into the KD-tree makes distances a latitude-dependent mix of degrees, skewing weights east–west versus north–south. Reproject to the local UTM zone first — this is the most common silent error in IDW pipelines.

Related #

• Spatial Interpolation with Kriging and IDW for Sensor Fields — the parent workflow covering projection, gridding, and validation for both interpolators
• Ordinary Kriging vs IDW for Sparse Sensor Networks — when to reach past IDW for kriging’s variance estimate
• Writing Cloud-Optimized GeoTIFFs from Sensor Grids — the next step: write the finished IDW surface to a COG for tiling and archive