Geostatistics — Variograms and Kriging

Points, Not Regions

The areal examples lived on regions with a neighbour graph. Geostatistics handles the other kind of spatial data — measurements at continuous coordinates, where distance is meaningful and the object of interest is the surface between the observations. The data are the Meuse floodplain: 155 soil samples with zinc from 113 to 1839 ppm, high along the river channel where flooding deposits it, falling off inland.

The variogram describes how similarity decays with distance, rising from a nugget (micro-scale noise) toward a sill (the overall variance) across a range. Fitting an exponential model gives a range of 278 m — a practical range near 833 m, since the exponential reaches 95% of its sill at roughly three times the parameter. That fitted covariance is what kriging uses to weight the data.

Predicting the surface, and its uncertainty

Kriging is the best linear unbiased predictor at an unmeasured location: a weighted average of the data, with weights from the covariance and a variance attached to every prediction — the honest "how sure are we here?" map, darkest at the samples and brightest in the gaps. Ordinary kriging assumes an unknown constant mean; universal kriging lets the mean be a regression on covariates. Adding distance-to-river as a trend lowers leave-one-out error from 0.402 to 0.383, and R's gstat agrees closely at 0.393 and 0.378 — the deterministic trend explaining part of what ordinary kriging left to pure smoothing.

How firmly is the variogram pinned down?

A caution that the cross-engine comparison forces into the open. The from-scratch fit reports a range of 278 m; gstat reports 450 m on the same 155 samples. Neither is wrong — the binning cutoff and the fitting weights are analyst choices, and moving them from (half the maximum pair distance, $\sqrt{N_j}$) to gstat's defaults (a third of the bounding-box diagonal, $N_j/h_j^2$) walks the fitted range from 278 m to 896 m. What does not move is the prediction: the surfaces agree to 0.99 and the leave-one-out errors to about 0.01, because range and sill trade off along a ridge that leaves the predictor nearly unchanged. A fitted range describes the fit, not the soil.

Kriging is a Gaussian process

The identity that ties this example to the rest of the collection: kriging is Gaussian-process regression. The covariance is the GP kernel, the kriging predictor is the posterior mean, and the kriging variance is the posterior variance — the 1-D GP example is the same machinery in one dimension. PyMC's GP reproduces the kriged surface at a correlation of 0.985, and once the parameterisations are matched the kernels nearly agree too: 335 m against 278 m, nugget 0.032 against 0.000. Most of the apparent disagreement was units — PyMC's Exponential kernel is $e^{-r/2\ell}$, so its lengthscale is twice a variogram range, and its eta and sigma are standard deviations where nugget and sill are variances. What genuinely differs is the process variance (1.21 against 0.60), largely the $\text{Gamma}(2,2)$ prior pulling the lengthscale up.

Notebooks

Downloads

Geostatistics Module — Source Code

"""
kriging.py -- GEOSTATISTICS: variograms and kriging for point-referenced data (from scratch).

Backs the notebooks in  "Geostatistics -- Variograms and Kriging".

Areal models (CAR/SAR) live on regions with a neighbour graph. GEOSTATISTICS instead handles
POINT-REFERENCED data -- measurements at continuous coordinates (soil samples, weather stations,
house locations) -- and asks two things: how does similarity decay with distance, and what is the
value at an UNMEASURED location?

The VARIOGRAM answers the first. The empirical semivariogram

    gamma(h) = (1 / 2 N(h)) * sum over pairs at distance ~ h of (z_i - z_j)^2

rises from a NUGGET (micro-scale noise) toward a SILL (the overall variance) over a RANGE (the
distance beyond which points are effectively uncorrelated). Fitting a model (exponential, Matern,
spherical) to it gives the covariance C(h) = sill - gamma(h) used for prediction.

KRIGING answers the second: it is the best linear unbiased predictor (BLUP) at a new location, a
weighted average of the data with weights from the covariance. ORDINARY kriging assumes an unknown
constant mean (a Lagrange constraint forces the weights to sum to one); UNIVERSAL kriging lets the
mean be a regression on covariates (a spatial trend). Each prediction comes with a KRIGING VARIANCE
-- the map of uncertainty, largest far from data.

The key identity for this arc: KRIGING IS GAUSSIAN-PROCESS REGRESSION. C(h) is the GP covariance
kernel, the kriging predictor is the GP posterior mean, and the kriging variance is the GP posterior
variance. The 1-D GP notebooks in "Gaussian-Process Regression" are the same machinery; here it
is 2-D and spatial, and the variogram is just how geostatisticians estimate the kernel.
"""

import numpy as np
from scipy.spatial.distance import cdist, pdist, squareform
from scipy.optimize import least_squares


# --------------------------------------------------------------------------- #
#  variogram                                                                    #
# --------------------------------------------------------------------------- #

def empirical_variogram(coords, z, nbins=15, maxdist=None):
    """binned empirical semivariogram. Returns (bin centres, semivariance, pair counts)."""
    z = np.asarray(z, float); D = squareform(pdist(coords)); iu = np.triu_indices(len(z), 1)
    h = D[iu]; g = 0.5 * (z[iu[0]] - z[iu[1]]) ** 2
    if maxdist is None:
        maxdist = h.max() / 2
    edges = np.linspace(0, maxdist, nbins + 1); cen = 0.5 * (edges[:-1] + edges[1:])
    gam = np.full(nbins, np.nan); cnt = np.zeros(nbins)
    for b in range(nbins):
        m = (h > edges[b]) & (h <= edges[b + 1])
        if m.sum() > 0:
            gam[b] = g[m].mean(); cnt[b] = m.sum()
    ok = ~np.isnan(gam)
    return cen[ok], gam[ok], cnt[ok]


def vgm_exponential(h, nugget, psill, rng):
    """exponential semivariogram: nugget at h>0, rising to nugget+psill over ~3*rng."""
    return nugget + psill * (1 - np.exp(-h / rng))


def fit_variogram(cen, gam, cnt):
    """weighted least-squares fit of an exponential variogram (nugget, partial sill, range)."""
    w = np.sqrt(cnt)
    def resid(p):
        nug, ps, rg = np.abs(p)
        return w * (vgm_exponential(cen, nug, ps, rg) - gam)
    p0 = [gam.min(), gam.max() - gam.min(), cen.max() / 3]
    r = least_squares(resid, p0, method="lm")
    nug, ps, rg = np.abs(r.x)
    return dict(nugget=nug, psill=ps, rng=rg, sill=nug + ps)


def _cov(D, v):
    """covariance matrix from an exponential variogram: C(h)=psill*exp(-h/rng); C(0)=psill+nugget."""
    C = v["psill"] * np.exp(-D / v["rng"])
    C[D == 0] = v["psill"] + v["nugget"]
    return C


# --------------------------------------------------------------------------- #
#  kriging                                                                       #
# --------------------------------------------------------------------------- #

def ordinary_kriging(coords, z, grid, v):
    """ordinary kriging (unknown constant mean). Returns prediction and kriging-variance arrays."""
    z = np.asarray(z, float); n = len(z)
    K = _cov(squareform(pdist(coords)), v)
    A = np.zeros((n + 1, n + 1)); A[:n, :n] = K; A[:n, n] = 1; A[n, :n] = 1
    Ai = np.linalg.inv(A)
    k0 = _cov(cdist(coords, grid), v)                              # (n, m)
    b = np.vstack([k0, np.ones((1, grid.shape[0]))])               # (n+1, m)
    W = Ai @ b                                                     # weights + Lagrange mult
    pred = W[:n].T @ z
    var = (v["psill"] + v["nugget"]) - np.sum(b * W, axis=0)
    return pred, np.maximum(var, 0)


def universal_kriging(coords, z, Xd, grid, Xg, v):
    """universal kriging with a linear trend. Xd:(n,p) trend basis at data (incl. intercept column),
    Xg:(m,p) at the grid. Returns prediction and kriging variance."""
    z = np.asarray(z, float); n, p = Xd.shape
    K = _cov(squareform(pdist(coords)), v)
    A = np.zeros((n + p, n + p)); A[:n, :n] = K; A[:n, n:] = Xd; A[n:, :n] = Xd.T
    Ai = np.linalg.inv(A)
    k0 = _cov(cdist(coords, grid), v)
    b = np.vstack([k0, Xg.T])                                      # (n+p, m)
    W = Ai @ b
    pred = W[:n].T @ z
    var = (v["psill"] + v["nugget"]) - np.sum(b * W, axis=0)
    return pred, np.maximum(var, 0)

References

• Matheron, G. (1963). Principles of geostatistics. Economic Geology 58(8), 1246–1266. — the variogram and kriging
• Cressie, N. (1993). Statistics for Spatial Data, revised ed. Wiley. — the standard treatment, including variogram-fitting weights
• Diggle, P. J. & Ribeiro, P. J. (2007). Model-based Geostatistics. Springer. — the likelihood view that makes the GP equivalence explicit
• Rasmussen, C. E. & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. — the same equations in the other vocabulary
• Pebesma, E. J. (2004). Multivariable geostatistics in S: the gstat package. Computers & Geosciences 30(7), 683–691. — the package used for the R cross-check
• Rikken, M. G. J. & Van Rijn, R. P. G. (1993). Soil Pollution with Heavy Metals. University of Utrecht. — the Meuse survey