Fitting and Validating a Variogram for Sensor Fields #

The variogram is the only thing separating kriging from a weighted average with extra steps. It is where the data itself tells the interpolator how quickly values decorrelate with distance, and a badly fitted one produces weights no better than IDW’s fixed rule while looking authoritative enough that nobody checks. This guide fits one properly for a sensor network — including the part most tutorials skip, which is deciding whether the fit is stable enough to use at all. It is the prerequisite for spatial interpolation with kriging and IDW for sensor fields.

Reading a Variogram #

Plot semivariance — half the mean squared difference between value pairs — against the distance separating them. Three parameters describe the resulting curve.

Nugget is the intercept: semivariance at (extrapolated) zero distance, which in theory should be zero. It is not, because two sensors at the same place would still disagree, by their measurement uncertainty plus any variation at scales finer than your shortest sensor spacing. For a low-cost PM2.5 network with 6 µg/m³ uncertainty, a nugget near 36 is exactly what physics predicts.

Sill is the plateau: the total variance of the field. Beyond this distance, knowing one sensor’s value tells you nothing about another’s.

Range is where the curve reaches the sill: the distance over which spatial correlation persists. It is the single most useful number the exercise produces, because it tells you how far apart sensors can be before they stop informing each other — which is a network design answer as much as an interpolation parameter.

The ratio of nugget to sill is the diagnostic. Below about 0.25 the field has strong spatial structure and kriging will do well. Above 0.75, most of the variance is noise or sub-grid variation, and no interpolator can reconstruct a field that your network cannot see.

Production-Ready Implementation #

# python 3.11 · numpy==1.26.4 · scipy==1.12.0 · pyproj==3.6.1
from __future__ import annotations

from dataclasses import dataclass

import numpy as np
from scipy.optimize import curve_fit


@dataclass(frozen=True)
class VariogramModel:
    nugget: float
    sill: float
    range_m: float
    model: str = "exponential"

    @property
    def nugget_ratio(self) -> float:
        return self.nugget / self.sill if self.sill > 0 else 1.0

    def gamma(self, h: np.ndarray) -> np.ndarray:
        """Semivariance predicted at separation distances h."""
        h = np.asarray(h, dtype=float)
        partial = self.sill - self.nugget
        if self.model == "exponential":
            return self.nugget + partial * (1 - np.exp(-3.0 * h / self.range_m))
        if self.model == "spherical":
            out = self.nugget + partial * (1.5 * h / self.range_m - 0.5 * (h / self.range_m) ** 3)
            return np.where(h < self.range_m, out, self.sill)
        if self.model == "gaussian":
            return self.nugget + partial * (1 - np.exp(-3.0 * (h / self.range_m) ** 2))
        raise ValueError(f"unknown model: {self.model}")


def experimental_variogram(
    xy: np.ndarray, values: np.ndarray, *, n_lags: int = 15,
    max_distance: float | None = None, min_pairs: int = 30,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Empirical semivariance by lag bin.

    Bins with fewer than `min_pairs` are dropped rather than plotted: a lag
    estimated from four pairs is noise, and fitting through it is how an
    unstable variogram gets a confident-looking curve.
    """
    n = len(values)
    iu = np.triu_indices(n, k=1)
    d = np.hypot(xy[iu[0], 0] - xy[iu[1], 0], xy[iu[0], 1] - xy[iu[1], 1])
    sq = 0.5 * (values[iu[0]] - values[iu[1]]) ** 2

    cutoff = max_distance or float(np.percentile(d, 50))    # half the max distance is convention
    keep = d <= cutoff
    d, sq = d[keep], sq[keep]

    edges = np.linspace(0, cutoff, n_lags + 1)
    idx = np.digitize(d, edges) - 1
    lags, gammas, counts = [], [], []
    for b in range(n_lags):
        sel = idx == b
        if sel.sum() < min_pairs:
            continue
        lags.append(float(d[sel].mean()))
        gammas.append(float(sq[sel].mean()))
        counts.append(int(sel.sum()))
    return np.array(lags), np.array(gammas), np.array(counts)


def fit_variogram(lags: np.ndarray, gammas: np.ndarray, counts: np.ndarray,
                  model: str = "exponential") -> VariogramModel:
    """Weighted least-squares fit, weighting each lag by its pair count.

    Unweighted fitting lets a far lag estimated from 40 pairs pull the curve as
    hard as a near lag estimated from 4 000 — and the near lags are the ones that
    actually determine the kriging weights.
    """
    sill0, range0 = float(np.max(gammas)), float(np.max(lags) / 2)
    nugget0 = float(min(gammas[0], 0.5 * sill0))

    def f(h, nugget, sill, rng):
        return VariogramModel(nugget, sill, rng, model).gamma(h)

    bounds = ([0, 1e-9, 1e-6], [sill0 * 1.5, sill0 * 3, float(np.max(lags)) * 3])
    popt, _ = curve_fit(f, lags, gammas, p0=[nugget0, sill0, range0],
                        sigma=1.0 / np.sqrt(counts), bounds=bounds, maxfev=20_000)
    return VariogramModel(nugget=popt[0], sill=popt[1], range_m=popt[2], model=model)

Choosing between models is an empirical question, answered by comparing cross-validated error rather than by preference:

def choose_model(xy, values, candidates=("exponential", "spherical", "gaussian")) -> VariogramModel:
    """Fit each model and keep the one with the lowest leave-one-out RMSE."""
    lags, gammas, counts = experimental_variogram(xy, values)
    best, best_rmse = None, float("inf")
    for name in candidates:
        try:
            fitted = fit_variogram(lags, gammas, counts, model=name)
        except RuntimeError:
            continue                       # this model would not converge on this data
        rmse = loo_rmse(xy, values, fitted)
        if rmse < best_rmse:
            best, best_rmse = fitted, rmse
    if best is None:
        raise RuntimeError("no variogram model converged — use IDW instead")
    return best

Parameter Tuning Guide #

Verification and Testing #

Cross-validation is the only real test: refit without each sensor and predict it.

# python 3.11 · numpy==1.26.4 · pytest==8.2.0
def loo_rmse(xy: np.ndarray, values: np.ndarray, model: VariogramModel) -> float:
    """Leave-one-out kriging RMSE under a fitted variogram."""
    errors = []
    for i in range(len(values)):
        mask = np.arange(len(values)) != i
        pred = ordinary_kriging(xy[mask], values[mask], xy[i], model)
        errors.append(pred - values[i])
    return float(np.sqrt(np.mean(np.square(errors))))


def test_fit_recovers_a_known_synthetic_variogram():
    """Simulate a field with a known range; the fit must find it within 30%."""
    xy, values = simulate_field(n=120, nugget=4.0, sill=40.0, range_m=8000.0, seed=5)
    lags, gammas, counts = experimental_variogram(xy, values)
    fitted = fit_variogram(lags, gammas, counts)
    assert 0.7 * 8000 < fitted.range_m < 1.3 * 8000
    assert 0.5 * 4.0 < fitted.nugget < 2.0 * 4.0


def test_unstable_fit_is_detected_rather_than_used():
    """With 8 sensors the fit is arbitrary — the pipeline must fall back to IDW."""
    xy, values = simulate_field(n=8, nugget=4.0, sill=40.0, range_m=8000.0, seed=9)
    lags, gammas, counts = experimental_variogram(xy, values, min_pairs=30)
    assert len(lags) < 4                              # too few usable bins to fit


def test_kriging_beats_idw_only_when_structure_exists(strong_field, weak_field):
    for field, expect_kriging_better in ((strong_field, True), (weak_field, False)):
        xy, values = field
        model = choose_model(xy, values)
        assert (loo_rmse(xy, values, model) < idw_loo_rmse(xy, values)) is expect_kriging_better

That last test encodes the honest position: kriging is not universally better, and the variogram is what tells you which case you are in.

Gotchas #

Fitting on geographic coordinates. Distances in degrees are not distances — a degree of longitude is 111 km at the equator and 71 km at 50° north, so a variogram fitted in degrees is anisotropic by construction. Project to a metric CRS first, as the CRS handling stage describes.

Ignoring anisotropy. A pollution field along a valley or a prevailing wind has a longer range in one direction. An isotropic variogram averages the two and under-fits both. Compute directional variograms if your residuals show orientation.

Fitting through under-populated far lags. They have the widest error bars and the greatest leverage on the fitted range. Cut at half the maximum separation and weight by pair count.

Treating the nugget as a nuisance to be set to zero. For low-cost sensor networks the nugget is mostly real measurement uncertainty. Forcing it to zero makes kriging interpolate exactly through every noisy observation, producing a surface with a spike at every sensor.

FAQ #

Related #

• Spatial Interpolation with Kriging and IDW for Sensor Fields — the interpolation stage the variogram feeds
• Ordinary Kriging vs IDW for Sparse Sensor Networks — the decision this fit either supports or rules out
• Building an IDW Air Quality Surface from Point Sensors in Python — the fallback when the variogram will not fit stably