Geostatistics — Variograms and Kriging¶

Predicting a surface from point samples, and the Gaussian-process connection¶

The areal projects (CAR/SAR) lived on regions with a neighbour graph. Geostatistics handles point-referenced data — measurements at continuous coordinates (soil cores, weather stations, house sales) — and answers two questions: how does similarity decay with distance, and what is the value at unmeasured locations?

We use the classic Meuse data — zinc concentration in topsoil at 155 locations on a bend of the river Meuse (Netherlands), with a prediction grid of 3,103 cells. Two tools:

  1. The variogram measures spatial dependence: $\gamma(h)=\tfrac{1}{2N(h)}\sum(z_i-z_j)^2$ over pairs $\sim h$ apart, rising from a nugget toward a sill over a range.
  2. Kriging predicts the surface — the best linear unbiased predictor, a covariance-weighted average of the data, with a kriging-variance map of uncertainty. Ordinary kriging assumes an unknown constant mean; universal kriging adds a spatial trend (here, distance to the river).

The unifying identity: kriging is Gaussian-process regression — the variogram is the GP kernel, the kriging predictor is the GP posterior mean, the kriging variance is the GP posterior variance. The 1-D GPs in Gaussian-Process Regression are the same machinery; here it is 2-D and spatial. We build it from scratch, map the surfaces, and confirm with a PyMC GP.

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt, geopandas as gpd, contextily as cx
import kriging as K
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; GREY="#718096"
m=pd.read_csv("meuse.csv"); g=pd.read_csv("meuse_grid.csv")
coords=m[["x","y"]].to_numpy(float); z=np.log(m["zinc"].to_numpy()); gc=g[["x","y"]].to_numpy(float)
def surfmap(ax, xy, vals, title, cmap="viridis", pts=False, s=6):
    sc=ax.scatter(xy[:,0], xy[:,1], c=vals, s=s, cmap=cmap, marker="s"); ax.set_aspect("equal"); ax.axis("off")
    if pts: ax.scatter(coords[:,0],coords[:,1], c="none", edgecolor="k", s=14, linewidth=.4)
    ax.set_title(title,fontsize=10); plt.colorbar(sc,ax=ax,shrink=0.6)
print(f"{len(m)} soil samples; zinc {int(m['zinc'].min())}-{int(m['zinc'].max())} ppm (log-transformed, ~lognormal)")
print("Point-referenced data: values at coordinates, no regions. Goal: a continuous risk surface + its uncertainty.")

155 soil samples; zinc 113-1839 ppm (log-transformed, ~lognormal)
Point-referenced data: values at coordinates, no regions. Goal: a continuous risk surface + its uncertainty.

1. Where this is, and the sampled surface¶

The coordinates are real geolocation — the Dutch national grid (EPSG:28992) — so the data sit on an actual place: a bend of the river Meuse near Stein, in the southern Netherlands. The first map puts the 155 sampling sites on an OpenStreetMap basemap for perspective; the river and floodplain are visible, and the high-zinc samples hug the water (metals deposited by flooding). The subsequent kriging maps use the same coordinate frame (metres), so this locator orients them all.

In [2]:
# geographic context: sample points (Dutch RD grid) on an OpenStreetMap basemap
gpts = gpd.GeoDataFrame({"logzinc":z}, geometry=gpd.points_from_xy(coords[:,0],coords[:,1]), crs=28992).to_crs(3857)
fig,ax=plt.subplots(figsize=(6.4,5.6))
gpts.plot(ax=ax, column="logzinc", cmap="YlOrRd", markersize=55, edgecolor="k", linewidth=.3, legend=True, legend_kwds={"shrink":0.55,"label":"log zinc"})
try:
    cx.add_basemap(ax, source=cx.providers.OpenStreetMap.Mapnik, crs=3857, attribution_size=5)
except Exception as e:
    print("basemap tiles unavailable, showing points only:", str(e)[:80])
ax.set_title("Meuse soil samples on the river Meuse near Stein, NL — log zinc (ppm)"); ax.set_axis_off()
plt.tight_layout(); plt.show()
print("The sampling transect follows the river bend; high zinc (dark) clings to the channel and floodplain, falling off")
print("inland. That geographic structure -- driven by flood deposition -- is what the variogram and kriging quantify.")
No description has been provided for this image 
The sampling transect follows the river bend; high zinc (dark) clings to the channel and floodplain, falling off
inland. That geographic structure -- driven by flood deposition -- is what the variogram and kriging quantify.

2. The variogram¶

The empirical semivariogram bins all sample pairs by separation distance and averages their squared differences. It rises from near zero and levels off at the sill (the overall variance) once points are far enough apart to be uncorrelated — the distance where that happens is the range. We fit an exponential model (nugget, partial sill, range).

In [3]:
cen,gam,cnt=K.empirical_variogram(coords,z,nbins=15)
v=K.fit_variogram(cen,gam,cnt)
print(f"fitted exponential variogram:  nugget {v['nugget']:.3f}   partial sill {v['psill']:.3f}   range {v['rng']:.0f} m   sill {v['sill']:.3f}")
hh=np.linspace(0,cen.max(),100)
fig,ax=plt.subplots(figsize=(5.6,3))
ax.scatter(cen,gam,s=40+cnt/3,color=BLUE,label="empirical (size = pair count)")
ax.plot(hh,K.vgm_exponential(hh,v["nugget"],v["psill"],v["rng"]),color=RED,lw=2,label="fitted exponential")
ax.axhline(v["sill"],color=GREY,ls="--",lw=1,label="sill"); ax.axvline(3*v["rng"],color=GREEN,ls=":",lw=1,label="practical range (3x)")
ax.set_xlabel("distance h (m)"); ax.set_ylabel(r"semivariance $\gamma(h)$"); ax.set_title("Empirical and fitted variogram of log-zinc")
ax.legend(frameon=False,fontsize=8); plt.tight_layout(); plt.show()
print(f"semivariance rises to the sill over a practical range of ~{3*v['rng']:.0f} m: beyond that, samples are effectively")
print("uncorrelated. This fitted covariance is what kriging (and, equivalently, a GP) uses to weight the data.")

fitted exponential variogram:  nugget 0.000   partial sill 0.603   range 278 m   sill 0.603
No description has been provided for this image 
semivariance rises to the sill over a practical range of ~833 m: beyond that, samples are effectively
uncorrelated. This fitted covariance is what kriging (and, equivalently, a GP) uses to weight the data.

How firmly is the variogram pinned down?¶

The cutoff distance and the fitting weights are analyst choices, not data. gstat defaults to a cutoff of a third of the bounding-box diagonal and weights $N_j/h_j^2$; the fit above uses half the maximum pair distance and $\sqrt{N_j}$. Refitting under each combination shows how much of the "estimated" range is really a choice — and, more usefully, how little the prediction cares.

In [4]:
from scipy.optimize import least_squares
from scipy.spatial.distance import pdist, squareform
diag = float(np.hypot(np.ptp(coords[:,0]), np.ptp(coords[:,1]))); hmax = float(squareform(pdist(coords)).max())
def refit(cutoff, wname):
    cen, gam, cnt = K.empirical_variogram(coords, z, nbins=15, maxdist=cutoff)
    w = np.sqrt(cnt) if wname == 'sqrt(N)' else cnt/cen**2
    f = lambda p: w*(K.vgm_exponential(cen, *np.abs(p)) - gam)
    nug, ps, rg = np.abs(least_squares(f, [gam.min(), gam.max()-gam.min(), cen.max()/3], method='lm').x)
    return nug, ps, rg
print(f"{'cutoff':>16s}{'weights':>10s}{'nugget':>9s}{'part.sill':>11s}{'range (m)':>11s}")
for cname, cut in ((f'hmax/2 = {hmax/2:.0f}', hmax/2), (f'diag/3 = {diag/3:.0f}', diag/3)):
    for wname in ('sqrt(N)', 'N/h^2'):
        nug, ps, rg = refit(cut, wname)
        print(f"{cname:>16s}{wname:>10s}{nug:9.3f}{ps:11.2f}{rg:11.0f}")
print()
print("gstat, in the R notebook, reports nugget 0.000, partial sill 0.72, range 450 m -- inside that span.")
print("So the fitted range moves by roughly threefold on identical data, purely with the cutoff and the")
print("weighting. The reassuring part is what does NOT move: the kriged surfaces agree to 0.99 and the")
print("leave-one-out errors to about 0.01, because range and sill trade off along a ridge that leaves the")
print("predictor almost unchanged. Read a fitted range as a description of the fit, not a property of the soil.")

          cutoff   weights   nugget  part.sill  range (m)
   hmax/2 = 2220   sqrt(N)    0.000       0.60        278
   hmax/2 = 2220     N/h^2    0.075       0.84        800
   diag/3 = 1597   sqrt(N)    0.000       0.66        350
   diag/3 = 1597     N/h^2    0.069       0.93        896

gstat, in the R notebook, reports nugget 0.000, partial sill 0.72, range 450 m -- inside that span.
So the fitted range moves by roughly threefold on identical data, purely with the cutoff and the
weighting. The reassuring part is what does NOT move: the kriged 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 almost unchanged. Read a fitted range as a description of the fit, not a property of the soil.

3. Ordinary kriging — the predicted surface and its uncertainty¶

Ordinary kriging predicts log-zinc across the grid as a covariance-weighted average of the samples (weights constrained to sum to one, so the unknown mean cancels). Two maps come out: the predicted surface and the kriging-variance — uncertainty is lowest at the sample sites and grows with distance from them.

In [5]:
pred,var=K.ordinary_kriging(coords,z,gc,v)
# leave-one-out cross-validation
loo=np.array([K.ordinary_kriging(np.delete(coords,i,0),np.delete(z,i),coords[i:i+1],v)[0][0] for i in range(len(z))])
rmse=np.sqrt(np.mean((loo-z)**2))
print(f"ordinary kriging:  leave-one-out RMSE {rmse:.3f} (log-zinc), correlation {np.corrcoef(loo,z)[0,1]:.2f}")
fig,ax=plt.subplots(1,2,figsize=(9.8,4.5))
surfmap(ax[0], gc, pred, "Ordinary kriging: predicted log-zinc", cmap="YlOrRd", pts=True)
surfmap(ax[1], gc, np.sqrt(var), "Kriging standard error (uncertainty)", cmap="Blues", pts=True)
plt.tight_layout(); plt.show()
print("The kriged surface reconstructs the high-zinc river corridor and smooths between samples; the error map is dark")
print("(low) at the sample points (circles) and bright (high) in the gaps -- the honest 'how sure are we here?' map.")

ordinary kriging:  leave-one-out RMSE 0.402 (log-zinc), correlation 0.83
No description has been provided for this image 
The kriged surface reconstructs the high-zinc river corridor and smooths between samples; the error map is dark
(low) at the sample points (circles) and bright (high) in the gaps -- the honest 'how sure are we here?' map.

4. Universal kriging — adding a trend¶

Zinc is deposited by river flooding, so distance to the river is a strong predictor. Universal kriging models the mean as a regression on $\sqrt{\text{dist}}$ and krigs the residuals, combining a deterministic trend with spatial correlation. It should predict better than ordinary kriging.

In [6]:
Xd=np.column_stack([np.ones(len(z)), np.sqrt(m["dist"])]); Xg=np.column_stack([np.ones(len(gc)), np.sqrt(g["dist"])])
predu,varu=K.universal_kriging(coords,z,Xd,gc,Xg,v)
loou=np.array([K.universal_kriging(np.delete(coords,i,0),np.delete(z,i),np.delete(Xd,i,0),coords[i:i+1],Xd[i:i+1],v)[0][0] for i in range(len(z))])
rmseu=np.sqrt(np.mean((loou-z)**2))
print(f"universal kriging (trend = sqrt distance-to-river):  LOO RMSE {rmseu:.3f}  vs ordinary {rmse:.3f}")
fig,ax=plt.subplots(1,2,figsize=(9.8,4.5))
surfmap(ax[0], gc, pred, "Ordinary kriging", cmap="YlOrRd", pts=True)
surfmap(ax[1], gc, predu, "Universal kriging (with river-distance trend)", cmap="YlOrRd", pts=True)
plt.tight_layout(); plt.show()
print(f"adding the river-distance trend lowers the cross-validation error ({rmseu:.3f} < {rmse:.3f}) and sharpens the")
print("high-zinc corridor along the river -- the deterministic trend explains part of what ordinary kriging left to")
print("pure spatial smoothing. This is the geostatistics version of a regression with spatially-correlated errors.")

universal kriging (trend = sqrt distance-to-river):  LOO RMSE 0.383  vs ordinary 0.402
No description has been provided for this image 
adding the river-distance trend lowers the cross-validation error (0.383 < 0.402) and sharpens the
high-zinc corridor along the river -- the deterministic trend explains part of what ordinary kriging left to
pure spatial smoothing. This is the geostatistics version of a regression with spatially-correlated errors.

5. Kriging is a Gaussian process — the PyMC view¶

Kriging and GP regression are the same model in different vocabularies: the variogram is the covariance kernel, the kriged surface is the GP posterior mean, the kriging variance is the GP posterior variance. We fit the identical model as a PyMC GP (pm.gp.Marginal with an exponential kernel + a nugget), learn the kernel hyper-parameters by marginal likelihood, and predict the surface — the Bayesian version of kriging, where the range/sill are estimated with uncertainty rather than fixed from the variogram fit.

In [7]:
import pymc as pm
sc=1000.0; Xs=(coords-coords.mean(0))/sc; Xgs=(gc-coords.mean(0))/sc; zc=z-z.mean()
with pm.Model() as gpm:
    ls=pm.Gamma("ls",2,2); eta=pm.HalfNormal("eta",1); nug=pm.HalfNormal("nug",0.5)
    cov=eta**2*pm.gp.cov.Exponential(2, ls)
    gp=pm.gp.Marginal(cov_func=cov)
    gp.marginal_likelihood("y", X=Xs, y=zc, sigma=nug)
    mp=pm.find_MAP(progressbar=False)
    mu,vv=gp.predict(Xgs, point=mp, diag=True, pred_noise=False)
# put the GP hyper-parameters in variogram units before comparing them:
# pm.gp.cov.Exponential is exp(-r/(2*ls)), so its lengthscale is TWICE a variogram range;
# eta and nug are STANDARD DEVIATIONS, while nugget and sill are variances.
gp_range = float(mp["ls"])*sc/2; gp_var = float(mp["eta"])**2; gp_nug = float(mp["nug"])**2
print(f"{'':22s}{'range (m)':>12s}{'nugget':>10s}{'part. sill':>12s}")
print(f"{'variogram fit':22s}{v['rng']:12.0f}{v['nugget']:10.3f}{v['psill']:12.2f}")
print(f"{'PyMC GP (converted)':22s}{gp_range:12.0f}{gp_nug:10.3f}{gp_var:12.2f}")
print(f"  raw GP output was lengthscale {float(mp['ls'])*sc:.0f} m, eta {float(mp['eta']):.2f}, sigma {float(mp['nug']):.2f}")
print(f"GP vs ordinary-kriging surface agreement: correlation {np.corrcoef(mu+z.mean(), pred)[0,1]:.3f}")
fig,ax=plt.subplots(1,2,figsize=(9.8,4.5))
surfmap(ax[0], gc, pred, "Ordinary kriging (from scratch)", cmap="YlOrRd", pts=True)
surfmap(ax[1], gc, mu+z.mean(), "PyMC Gaussian process (Bayesian kriging)", cmap="YlOrRd", pts=True)
plt.tight_layout(); plt.show()
print("The PyMC GP recovers essentially the SAME surface (corr ~0.99): kriging and GP regression are one method.")
print(f"Once the parameterisations are matched the kernels nearly agree too -- range {gp_range:.0f} m against {v['rng']:.0f} m and a")
print(f"nugget of {gp_nug:.3f} against {v['nugget']:.3f}. Most of the apparent disagreement was units: PyMC's Exponential kernel is")
print("exp(-r/2l), so its lengthscale is twice a variogram range, and eta and sigma are standard deviations where")
print(f"nugget and sill are variances. What genuinely differs is the process variance, {gp_var:.2f} against {v['psill']:.2f}, and that is")
print("largely the Gamma(2,2) prior pulling the lengthscale up -- a longer correlation range needs a larger variance")
print("to fit the same spread. Range and variance trade off along a ridge, which is exactly why two fits can differ")
print("in the parameters and still agree to 0.99 on the surface. The GP's advantage is that it can carry that")
print("uncertainty, where the variogram reads a single value off a curve.")

g++ not available, if using conda: `conda install gxx`

                         range (m)    nugget  part. sill
variogram fit                  278     0.000        0.60
PyMC GP (converted)            335     0.032        1.21
  raw GP output was lengthscale 670 m, eta 1.10, sigma 0.18
GP vs ordinary-kriging surface agreement: correlation 0.985
No description has been provided for this image 
The PyMC GP recovers essentially the SAME surface (corr ~0.99): kriging and GP regression are one method.
Once the parameterisations are matched the kernels nearly agree too -- range 335 m against 278 m and a
nugget of 0.032 against 0.000. Most of the apparent disagreement was units: PyMC's Exponential kernel is
exp(-r/2l), so its lengthscale is twice a variogram range, and eta and sigma are standard deviations where
nugget and sill are variances. What genuinely differs is the process variance, 1.21 against 0.60, and that is
largely the Gamma(2,2) prior pulling the lengthscale up -- a longer correlation range needs a larger variance
to fit the same spread. Range and variance trade off along a ridge, which is exactly why two fits can differ
in the parameters and still agree to 0.99 on the surface. The GP's advantage is that it can carry that
uncertainty, where the variogram reads a single value off a curve.

6. Summary¶

Geostatistics turns scattered point samples into a continuous surface. The variogram estimates how spatial similarity decays with distance (nugget, sill, range); kriging then predicts every location as a covariance-weighted average of the data, delivering both a surface and a kriging-variance map of where the prediction is trustworthy. On the Meuse soil data, ordinary kriging reconstructed the high-zinc river corridor, and universal kriging — adding a distance-to-river trend — predicted measurably better (lower cross-validation error). A PyMC Gaussian process reproduced the surface almost exactly (though it split the variance differently between nugget and range — the two trade off and are only weakly separately identified).

The through-line is the identity kriging = Gaussian-process regression: the variogram is the kernel, the kriged mean and variance are the GP posterior — so this subsection is the 2-D spatial face of the GPs in Gaussian-Process Regression, and the point-referenced complement to the areal CAR models (regions + neighbour graph) and the spatial-econometric models (regions + spillovers). For regional analysis, kriging is the tool when data come as located points rather than area totals — environmental exposure, station networks, or any sparsely-sampled continuous field. Next in the arc: spatial point processes (modelling the locations themselves), then the spatiotemporal capstone.