7.4 Geostatistical Interpolation

Part: 7 — Spatial Interpolation (Chapter 7.4 of 8)

Dataset: meuse.csv (155 topsoil samples, Meuse floodplain, EPSG:28992) predicted onto meuse_grid.csv (3,103 cells, 40 m). Most kriging uses log(zinc). Cokriging adds log(lead); regression kriging uses distance to the river; indicator kriging stays on zinc in ppm.

Focus: the semivariogram; ordinary, simple, universal, co-, regression and indicator kriging; prediction variance; calibration; agreement with pykrige

Library: spatialstats.interpolate (PySpatialStats)

---

Table of Contents

1. Overview

2. Python Setup

3. The Semivariogram

4. Empirical and Fitted Variograms

5. Ordinary Kriging

6. The Kriging Prediction Variance

7. Calibration: Standardized LOO Errors

8. Agreement with pykrige

9. Simple Kriging

10. Universal Kriging

11. Cokriging

12. Regression Kriging

13. Indicator Kriging and the E-type Estimate

14. Summary and Conclusions

15. Resources

---

1. Overview

This chapter delivers on Chapter 7.1's promise: a method whose predictions come with an honest, checkable measure of confidence, derived from an explicit model of the field's spatial correlation rather than from distance alone.

Chapter 7.3 interpolated zinc in ppm. Zinc on this floodplain is strongly right-skewed (113 to 1,839 ppm), so the kriging in Sections 5–12 is done on $\log(\mathrm{zinc})$, the scale on which a constant-mean model is plausible. Coordinates are already in metres (Amersfoort / RD New, EPSG:28992). Predictions are made at the cell centres in meuse_grid.csv, the same 40 m floodplain lattice as in Chapter 7.3.

Ordinary kriging is the workhorse: an unknown constant mean, weights from the variogram, and a prediction variance. The later sections change one assumption at a time. Simple kriging treats the mean as known. Universal kriging lets the mean be a polynomial in the coordinates. Cokriging borrows a second metal, lead, measured at the same sites. Regression kriging borrows distance to the river, which is known at every grid cell, not only at the samples. Indicator kriging changes the target itself, from the mean of log zinc to the probability that zinc exceeds a cutoff, and integrates those probabilities into an E-type mean on the ppm scale.

Step	What we do
Theory	The semivariogram $\gamma(h)$; nugget, partial sill, range; three common models
Fit	empirical_variogram, fit_variogram, fit_variogram_auto on the 155 Meuse samples
Krige	OrdinaryKriging — the BLUP under an unknown constant mean
Trust it	The kriging variance, and a calibration check confirming it is not just a number but a correct one
Cross-check	Agreement with the independent, widely used pykrige package
Change the mean	SimpleKriging (known mean) and UniversalKriging (a plane in $x$ and $y$)
Borrow more data	Cokriging with lead; RegressionKriging with distance to the river
Change the target	IndicatorKriging for $P(Z \ge t)$, then the E-type estimate of the mean in ppm

---

2. Python Setup

import sys, warnings
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=DeprecationWarning)


def find_repo_root(start: Path = Path.cwd()) -> Path:
    for p in [start, *start.parents]:
        if (p / "pyproject.toml").exists() and (p / "data").is_dir():
            return p
    raise FileNotFoundError("Run this notebook from inside the PySpatialStats repository.")


ROOT = find_repo_root()
DATA = ROOT / "data"
try:
    import spatialstats
except ImportError:
    sys.path.insert(0, str(ROOT))
    import spatialstats

from spatialstats.interpolate import (
    empirical_variogram, fit_variogram, fit_variogram_auto, VariogramModel,
    OrdinaryKriging, SimpleKriging, UniversalKriging, TrendSurface,
    empirical_cross_variogram, fit_lmc, Cokriging,
    RegressionKriging, IndicatorKriging,
)

plt.rcParams.update({
    "figure.dpi": 100, "savefig.dpi": 100, "font.size": 10, "axes.titlesize": 11,
    "axes.titleweight": "bold", "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25, "legend.frameon": False,
})
COLORS = {"blue": "#2b6cb0", "red": "#c0392b", "green": "#2f855a", "orange": "#dd6b20",
          "grey": "#718096", "purple": "#6b46c1", "teal": "#2c7a7b"}

meuse = pd.read_csv(DATA / "meuse.csv")
grid = pd.read_csv(DATA / "meuse_grid.csv")
coords = meuse[["x", "y"]].to_numpy(dtype=float)
zinc = meuse["zinc"].to_numpy(dtype=float)
values = np.log(zinc)
grid_xy = grid[["x", "y"]].to_numpy(dtype=float)

_xs = np.sort(grid["x"].unique())
_ys = np.sort(grid["y"].unique())
_dx, _dy = float(_xs[1] - _xs[0]), float(_ys[1] - _ys[0])
EXTENT = (_xs[0] - _dx / 2, _xs[-1] + _dx / 2, _ys[0] - _dy / 2, _ys[-1] + _dy / 2)


def to_raster(pred) -> np.ndarray:
    """Place predictions for the meuse_grid cells on the 40 m lattice (NaN off the floodplain)."""
    raster = np.full((len(_ys), len(_xs)), np.nan)
    raster[np.searchsorted(_ys, grid_xy[:, 1]), np.searchsorted(_xs, grid_xy[:, 0])] = pred
    return raster


print(f"spatialstats : version {spatialstats.__version__}")
print(f"n = {len(coords)} samples, {len(grid_xy)} prediction cells at {_dx:.0f} m")
print(f"zinc: {zinc.min():.0f} to {zinc.max():.0f} ppm; log(zinc): {values.min():.2f} to {values.max():.2f}")

spatialstats : version 0.3.0
n = 155 samples, 3103 prediction cells at 40 m
zinc: 113 to 1839 ppm; log(zinc): 4.73 to 7.52

---

3. The Semivariogram

$$ \gamma(h) = \tfrac{1}{2}\,\mathrm{Var}\!\left[Z(\mathbf{s}+h) - Z(\mathbf{s})\right] = C(0) - C(h) $$

$\gamma(h)$ measures how dissimilar two observations typically are as a function of the distance $h$ between them, under second-order stationarity (Chapter 7.1, Section 4). Three parameters describe a fitted variogram model:

• nugget ($C_0$): dissimilarity that persists even at $h \to 0$ — measurement error plus any spatial structure finer than the sampling network can resolve.
• partial sill ($C$): the additional variance explained by spatial correlation, on top of the nugget.
• range: the distance beyond which $\gamma(h)$ flattens out — pairs farther apart than the range are (approximately) spatially uncorrelated. This library reports the practical range for the exponential and Gaussian models (where $\gamma(h)$ only asymptotes, never truly flattens): the distance at which 95% of the sill is reached.

Three common models, each with a different rate of approach to the sill:

$$ \gamma_{\text{spherical}}(h) = C_0 + C\left[\tfrac{3h}{2a} - \tfrac{1}{2}\left(\tfrac{h}{a}\right)^3\right] \;(h \le a), \quad \gamma_{\text{exponential}}(h) = C_0 + C\left[1 - e^{-h/a}\right], \quad \gamma_{\text{gaussian}}(h) = C_0 + C\left[1 - e^{-(h/a)^2}\right] $$

---

4. Empirical and Fitted Variograms

empirical_variogram bins observation pairs by distance (Matheron's method-of-moments estimator) into n_lags lags; fit_variogram fits a chosen model to those bins by Cressie-weighted least squares (down-weighting noisy, sparsely populated far lags); fit_variogram_auto tries several models and returns the best one by fit quality.

empirical = empirical_variogram(coords, values, n_lags=15)
print(empirical.head(8))

vm_auto, empirical_auto = fit_variogram_auto(coords, values, models=("spherical", "exponential", "gaussian"))
print(f"\nauto-selected model: {vm_auto.model}")
print(f"nugget={vm_auto.nugget:.3f}   partial sill={vm_auto.psill:.3f}   range={vm_auto.range:,.0f} m")

           lag     gamma  n_pairs
0   112.027584  0.149697      158
1   226.807778  0.272436      518
2   372.695422  0.382132      659
3   515.797630  0.518178      722
4   665.862112  0.580283      799
5   813.100442  0.622774      803
6   960.810149  0.678342      779
7  1109.813743  0.676410      714

auto-selected model: spherical
nugget=0.033   partial sill=0.570   range=807 m

fig, ax = plt.subplots(figsize=(8, 5.5))
vm_auto.plot(empirical, ax=ax)
ax.set_title(f"Auto-fitted variogram: {vm_auto.model}\n"
            f"nugget={vm_auto.nugget:.3f}, partial sill={vm_auto.psill:.3f}, range={vm_auto.range:.0f} m")
ax.set_xlabel("distance (m)")
plt.tight_layout()
plt.savefig("variogram_fit.png", dpi=100, bbox_inches="tight")
plt.show()

The auto-fitted model is spherical, with a nugget of about 0.033, a partial sill of about 0.57, and a practical range of about 810 m. The nugget is much smaller than the partial sill, so most of the variance in log zinc is spatially structured rather than measurement noise or variation finer than the sample spacing. Samples more than about 800 m apart are essentially uninformative about each other. The floodplain is a few kilometres long, so that range covers a neighbourhood, not the whole map.

---

5. Ordinary Kriging

Ordinary kriging predicts as a weighted average $\hat{Z}(\mathbf{s}_0) = \sum_i \lambda_i Z(\mathbf{s}_i)$, with weights $\boldsymbol{\lambda}$ solved from the variogram under the constraint $\sum_i \lambda_i = 1$ (so the predictor is unbiased without needing to know the field's true mean in advance). Because the weights come from the variogram, not raw distance, clustered samples are automatically down-weighted relative to isolated ones carrying unique information — something IDW's purely distance-based weights cannot do.

ok = OrdinaryKriging(variogram="auto")
ok.fit(coords, values)
print(f"fitted variogram: {ok.variogram_.model}, nugget={ok.variogram_.nugget:.3f}, "
      f"psill={ok.variogram_.psill:.3f}, range={ok.variogram_.range:,.0f} m")

pred_at_data, std_at_data = ok.predict(coords, return_std=True)
print(f"exact at data points: max |pred - observed| = {np.max(np.abs(pred_at_data - values)):.2e}")

fitted variogram: spherical, nugget=0.033, psill=0.570, range=807 m
exact at data points: max |pred - observed| = 1.69e-14

Because $\gamma(0) = 0$ in the kriging system, ordinary kriging reproduces every sample exactly — confirmed above to machine precision — even though the fitted nugget is positive. The nugget is the jump in $\gamma(h)$ as $h \to 0$, not a noise term placed on the diagonal.

---

6. The Kriging Prediction Variance

pred_grid, std_grid = ok.predict(grid_xy, return_std=True)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
im0 = axes[0].imshow(to_raster(pred_grid), extent=EXTENT, origin="lower", cmap="viridis")
axes[0].scatter(coords[:, 0], coords[:, 1], s=10, color="white", edgecolor="k", linewidth=0.3)
axes[0].set_title("Kriging prediction (mean)")
axes[0].set_aspect("equal")
plt.colorbar(im0, ax=axes[0], fraction=0.046, label="log zinc")

im1 = axes[1].imshow(to_raster(std_grid), extent=EXTENT, origin="lower", cmap="magma")
axes[1].scatter(coords[:, 0], coords[:, 1], s=10, color="white", edgecolor="k", linewidth=0.3)
axes[1].set_title("Kriging standard deviation")
axes[1].set_aspect("equal")
plt.colorbar(im1, ax=axes[1], fraction=0.046, label="std (log zinc)")
plt.tight_layout()
plt.savefig("kriging_mean_variance.png", dpi=100, bbox_inches="tight")
plt.show()

print(f"prediction range: {pred_grid.min():.2f} to {pred_grid.max():.2f}")
print(f"kriging std range on the grid: {std_grid.min():.2f} to {std_grid.max():.2f}")

prediction range: 4.75 to 7.47
kriging std range on the grid: 0.25 to 0.71

The standard-deviation map is the payoff Chapter 7.1 promised: low near the samples, higher in the gaps of the floodplain — visibly independent of the log-zinc surface, and entirely a function of the sampling geometry and the variogram. On the 40 m grid the standard deviation runs from about 0.25, beside a sample, to about 0.71 at the least-supported cells. It is zero only at the samples themselves. This is not something any of Chapter 7.3's deterministic methods can produce, because none of them has an explicit model of spatial correlation to draw it from.

---

7. Calibration: Standardized LOO Errors

A prediction variance is only useful if it is correctly calibrated — if the true error, divided by the predicted standard deviation, behaves like a standard normal variable. Leave-one-out is the natural check: hold out each sample, krige it from the rest, and standardize.

n = len(coords)
std_errors = np.empty(n)
for i in range(n):
    mask = np.ones(n, dtype=bool)
    mask[i] = False
    m = OrdinaryKriging(variogram="auto")
    m.fit(coords[mask], values[mask])
    pred_i, std_i = m.predict(coords[i:i + 1], return_std=True)
    std_errors[i] = (values[i] - pred_i[0]) / std_i[0]

print(f"standardized LOO error: mean={std_errors.mean():.3f}, sd={std_errors.std():.3f}")
print(f"fraction within +-1.96 (nominal 95%): {(np.abs(std_errors) < 1.96).mean()*100:.1f}%")

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].hist(std_errors, bins=20, color=COLORS["blue"], edgecolor="white", density=True)
xs = np.linspace(-3.5, 3.5, 200)
axes[0].plot(xs, stats.norm.pdf(xs), color=COLORS["red"], linewidth=2, label="N(0,1)")
axes[0].set_title("Standardized LOO errors vs. N(0,1)")
axes[0].legend()
stats.probplot(std_errors, dist="norm", plot=axes[1])
axes[1].set_title("Normal Q-Q plot")
plt.tight_layout()
plt.savefig("kriging_calibration.png", dpi=100, bbox_inches="tight")
plt.show()

standardized LOO error: mean=-0.000, sd=0.948
fraction within +-1.96 (nominal 95%): 95.5%

Mean essentially 0 and standard deviation about 0.95, with about 95% of standardized errors inside $\pm 1.96$ — this is a genuinely calibrated uncertainty estimate on the log-zinc scale, not just a number that happens to be called "variance." Chapter 7.5's cross-validation extends this same idea (held-out prediction accuracy) to every method in this Part, not just kriging's own calibration.

---

8. Agreement with pykrige

As an independent correctness check, pykrige — a widely used, separately implemented kriging package — should agree with this library's own kriging system exactly, once fed the identical variogram parameters. One important gotcha: pykrige takes the total sill (nugget + partial sill), not the partial sill this library reports separately. The check below predicts at a sample of meuse_grid.csv cell centres.

from pykrige.ok import OrdinaryKriging as PyKrigeOK

vg = ok.variogram_
pk = PyKrigeOK(coords[:, 0], coords[:, 1], values, variogram_model="spherical",
              variogram_parameters={"nugget": vg.nugget, "sill": vg.nugget + vg.psill, "range": vg.range})

rng = np.random.RandomState(0)
test_pts = grid_xy[rng.choice(len(grid_xy), 20, replace=False)]

pred_ours = ok.predict(test_pts)
pred_pk, _ = pk.execute("points", test_pts[:, 0], test_pts[:, 1])
pred_pk = np.asarray(pred_pk).ravel()

print(f"max |ours - pykrige|: {np.max(np.abs(pred_ours - pred_pk)):.2e}")
print(f"ours   : {pred_ours[:5]}")
print(f"pykrige: {pred_pk[:5]}")

max |ours - pykrige|: 1.20e-13
ours   : [5.32580672 4.89096774 6.05453787 5.35674574 5.18858733]
pykrige: [5.32580672 4.89096774 6.05453787 5.35674574 5.18858733]

Agreement to machine precision confirms this library's kriging implementation is solving the same linear system pykrige does, not merely producing plausible-looking numbers — a genuine, independent correctness check on the core method this whole chapter relies on.

---

9. Simple Kriging

Ordinary kriging never asks for the mean. The constraint $\sum_i \lambda_i = 1$ estimates it locally, which is why a prediction far from every sample falls back to a neighbourhood mean rather than to a number you supplied. Simple kriging is the version in which the mean $m$ is known — from a process model, or from a survey much larger than the map being drawn. The weights then come from the covariance alone and need not sum to 1:

$$ \hat Z(\mathbf{s}_0) = m + \sum_i \lambda_i \bigl(Z(\mathbf{s}_i) - m\bigr), \qquad \mathbf{C}\,\boldsymbol{\lambda} = \mathbf{c}. $$

Past the variogram range every covariance is zero, so the prediction is exactly $m$ and the variance is the sill. Here $m$ is the mean of the 155 log-zinc values, and the variogram is the one already fitted for ordinary kriging. The only thing that changes is the treatment of the mean.

sk = SimpleKriging(variogram=ok.variogram_, mean=float(values.mean()))
sk.fit(coords, values)
sk_at_data = sk.predict(coords)
print(f"known mean m = {sk.mean_:.3f}")
print(f"exact at data points: max |pred - observed| = {np.max(np.abs(sk_at_data - values)):.2e}")

sk_grid = sk.predict(grid_xy)
print(f"max |simple − ordinary| on the grid: {np.max(np.abs(sk_grid - pred_grid)):.3f}")

lo = min(pred_grid.min(), sk_grid.min())
hi = max(pred_grid.max(), sk_grid.max())
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
for ax, surface, title in (
    (axes[0], pred_grid, "Ordinary kriging"),
    (axes[1], sk_grid, "Simple kriging"),
):
    im = ax.imshow(to_raster(surface), extent=EXTENT, origin="lower", cmap="viridis", vmin=lo, vmax=hi)
    ax.scatter(coords[:, 0], coords[:, 1], s=10, color="white", edgecolor="k", linewidth=0.3)
    ax.set_title(title)
    ax.set_aspect("equal")
    plt.colorbar(im, ax=ax, fraction=0.046, label="log zinc")
plt.tight_layout()
plt.savefig("simple_kriging.png", dpi=100, bbox_inches="tight")
plt.show()

known mean m = 5.886
exact at data points: max |pred - observed| = 1.78e-15
max |simple − ordinary| on the grid: 0.089

On this floodplain the two maps are almost the same: the largest gap is under a tenth of a log-zinc unit. That is the expected result, not a failed comparison. Every grid cell sits within a couple of variogram ranges of a sample, so ordinary kriging's locally estimated mean and the single global mean tell nearly the same story. Simple kriging separates from ordinary kriging where the data stop — off the edge of the sampled region — because only simple kriging is required to return exactly $m$ once the covariances vanish. Both remain exact at the samples.

---

10. Universal Kriging

A constant mean is the wrong model when the field has a trend in the coordinates. Universal kriging replaces the constant with a polynomial drift. Degree 1 is a plane in $x$ and $y$; degree 0 is ordinary kriging again, which is checked on the same 20 grid cells used in Section 8.

The variogram has to describe what is left after the drift, not the raw field. A degree-1 TrendSurface is fitted first, and fit_variogram_auto is applied to its residuals. That model is what UniversalKriging treats as the covariance of the stochastic part. The coordinates are centered and scaled inside the solver; the column space is still the plane $a + bx + cy$.

trend = TrendSurface(degree=1).fit(coords, values)
resid = values - trend.predict(coords)
vm_uk, _ = fit_variogram_auto(coords, resid)
print(f"degree-1 trend R² = {trend.r2_:.2f}")
print(f"residual variogram: {vm_uk.model}, nugget={vm_uk.nugget:.3f}, "
      f"psill={vm_uk.psill:.3f}, range={vm_uk.range:,.0f} m")
print(f"raw-field sill was {ok.variogram_.sill:.3f}; residual sill is {vm_uk.sill:.3f}")

uk = UniversalKriging(variogram=vm_uk, degree=1).fit(coords, values)
uk_grid = uk.predict(grid_xy)
uk0 = UniversalKriging(variogram=ok.variogram_, degree=0).fit(coords, values)
gap = np.max(np.abs(uk0.predict(test_pts) - ok.predict(test_pts)))
print(f"degree 0 vs ordinary kriging, max |diff| on 20 cells: {gap:.2e}")
print(f"max |universal − ordinary| on the grid: {np.max(np.abs(uk_grid - pred_grid)):.3f}")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
im0 = axes[0].imshow(to_raster(uk_grid), extent=EXTENT, origin="lower", cmap="viridis")
axes[0].scatter(coords[:, 0], coords[:, 1], s=10, color="white", edgecolor="k", linewidth=0.3)
axes[0].set_title("Universal kriging (degree 1)")
axes[0].set_aspect("equal")
plt.colorbar(im0, ax=axes[0], fraction=0.046, label="log zinc")

im1 = axes[1].imshow(to_raster(uk_grid - pred_grid), extent=EXTENT, origin="lower", cmap="coolwarm")
axes[1].scatter(coords[:, 0], coords[:, 1], s=10, color="k", linewidth=0.3)
axes[1].set_title("Universal − ordinary")
axes[1].set_aspect("equal")
plt.colorbar(im1, ax=axes[1], fraction=0.046, label="log zinc")
plt.tight_layout()
plt.savefig("universal_kriging.png", dpi=100, bbox_inches="tight")
plt.show()

degree-1 trend R² = 0.27
residual variogram: spherical, nugget=0.086, psill=0.368, range=1,070 m
raw-field sill was 0.603; residual sill is 0.454
degree 0 vs ordinary kriging, max |diff| on 20 cells: 2.31e-14
max |universal − ordinary| on the grid: 0.519

The plane accounts for about a quarter of the variance in log zinc ($R^2 \approx 0.27$), and the residual sill drops from about 0.60 to about 0.45. Degree 0 reproduces ordinary kriging to machine precision. Degree 1 does not: the two surfaces differ by up to about half a log unit where the plane is pulling the prediction and the local sample mean is not. That difference is the drift. It is also a warning. A coordinate polynomial will keep sloping outside the convex hull of the samples, so universal kriging is a model for a trend you believe, not a way to make the edges of the map look smoother.

---

11. Cokriging

Lead is measured at the same 155 sites as zinc, and the two logs move together (correlation about 0.97). Ordinary cokriging predicts zinc from both. There are two unbiasedness constraints: the zinc weights sum to 1, and the lead weights sum to 0, so lead contributes only through its cross-covariance with zinc. Lead is not on meuse_grid.csv. It does not need to be. The secondary enters as a set of samples, not as a map.

fit_lmc fits the two direct variograms and the cross-variogram with one shared practical range, then clips the cross-sill by the Cauchy–Schwarz bound $|c_{12}| \le \sqrt{c_1 c_2}$ so the coregionalization matrices can be positive semidefinite. info["lmc_check"] reports that check.

lead = np.log(meuse["lead"].to_numpy(dtype=float))
print(f"corr(log zinc, log lead) = {np.corrcoef(values, lead)[0, 1]:.3f}")

emp_cross = empirical_cross_variogram(coords, values, lead, n_lags=12)
vgs, cross, info = fit_lmc(
    coords, [values, lead], names=["log zinc", "log lead"], model="spherical", n_lags=12,
)
cv = cross[(0, 1)]
print(f"shared range = {info['range']:,.0f} m; LMC valid: {info['lmc_check']['valid']}")
print(f"cross-variogram: nugget={cv.nugget:.3f}, partial sill={cv.psill:.3f}")

ck = Cokriging(vgs[0], vgs[1], cv).fit(coords, values, coords, lead)
ck_at = ck.predict(coords)
print(f"exact at zinc samples: max |pred - observed| = {np.max(np.abs(ck_at - values)):.2e}")
ck_grid = ck.predict(grid_xy)
print(f"max |cokriging − ordinary| on the grid: {np.max(np.abs(ck_grid - pred_grid)):.3f}")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
axes[0].plot(emp_cross["lag"], emp_cross["gamma"], "o", color=COLORS["blue"], label="empirical")
hh = np.linspace(0, emp_cross["lag"].max(), 200)
axes[0].plot(hh, cv.gamma(hh), color=COLORS["red"], linewidth=2, label="fitted")
axes[0].set_xlabel("distance (m)")
axes[0].set_ylabel("cross-semivariance")
axes[0].set_title("Zinc–lead cross-variogram")
axes[0].legend()

im = axes[1].imshow(to_raster(ck_grid), extent=EXTENT, origin="lower", cmap="viridis", vmin=lo, vmax=hi)
axes[1].scatter(coords[:, 0], coords[:, 1], s=10, color="white", edgecolor="k", linewidth=0.3)
axes[1].set_title("Cokriging prediction")
axes[1].set_aspect("equal")
plt.colorbar(im, ax=axes[1], fraction=0.046, label="log zinc")
plt.tight_layout()
plt.savefig("cokriging.png", dpi=100, bbox_inches="tight")
plt.show()

corr(log zinc, log lead) = 0.967
shared range = 827 m; LMC valid: True
cross-variogram: nugget=0.032, partial sill=0.496
exact at zinc samples: max |pred - observed| = 2.84e-14
max |cokriging − ordinary| on the grid: 0.077

The linear model of coregionalization is valid, the shared range is about 830 m, and cokriging reproduces every zinc sample. On the grid it barely moves ordinary kriging: the largest difference is under a tenth of a log unit. Lead is informative, but it is measured at the same places as zinc, and zinc is already sampled densely relative to its range. Cokriging earns its keep when the secondary is denser than the primary — a cheap variable everywhere, the expensive one at a subset of sites — not when the two are duplicates of the same network. The cross-variogram is still worth fitting: it is the check that the two variables really share a spatial structure, rather than a correlation that is only global.

---

12. Regression Kriging

Distance to the river, dist, is a different kind of extra information. It is known at the 155 samples and at all 3,103 grid cells. Regression kriging uses it in two steps. A linear trend fits $\log(\mathrm{zinc})$ on dist. Ordinary kriging then interpolates the residuals, and the two are added back together. The standard deviation from predict(..., return_std=True) is the kriging standard deviation of those residuals. It does not include uncertainty in the regression coefficients.

Because the trend is known everywhere, this can change the map in the gaps between samples. Cokriging could not: lead stopped where the samples stopped.

x_dist = meuse[["dist"]].to_numpy(dtype=float)
x_grid = grid[["dist"]].to_numpy(dtype=float)
rk = RegressionKriging(regressor="linear", variogram="spherical", n_lags=12)
rk.fit(x_dist, values, coords)
slope = float(rk.regressor.coef_[0])
print(f"trend: log(zinc) = {rk.regressor.intercept_:.2f} + ({slope:.2f}) · dist")
print(f"residual variogram: nugget={rk.variogram_.nugget:.3f}, psill={rk.variogram_.psill:.3f}, "
      f"range={rk.variogram_.range:,.0f} m")
print(f"residual sill {rk.variogram_.sill:.3f} vs raw-field sill {ok.variogram_.sill:.3f}")

rk_grid, rk_std = rk.predict(x_grid, grid_xy, return_std=True)
print(f"max |regression − ordinary| on the grid: {np.max(np.abs(rk_grid - pred_grid)):.3f}")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
im0 = axes[0].imshow(to_raster(rk_grid), extent=EXTENT, origin="lower", cmap="viridis")
axes[0].scatter(coords[:, 0], coords[:, 1], s=10, color="white", edgecolor="k", linewidth=0.3)
axes[0].set_title("Regression kriging")
axes[0].set_aspect("equal")
plt.colorbar(im0, ax=axes[0], fraction=0.046, label="log zinc")

im1 = axes[1].imshow(to_raster(rk_std), extent=EXTENT, origin="lower", cmap="magma")
axes[1].scatter(coords[:, 0], coords[:, 1], s=10, color="white", edgecolor="k", linewidth=0.3)
axes[1].set_title("Residual kriging standard deviation")
axes[1].set_aspect("equal")
plt.colorbar(im1, ax=axes[1], fraction=0.046, label="std (log zinc)")
plt.tight_layout()
plt.savefig("regression_kriging.png", dpi=100, bbox_inches="tight")
plt.show()

trend: log(zinc) = 6.53 + (-2.70) · dist
residual variogram: nugget=0.072, psill=0.194, range=689 m
residual sill 0.267 vs raw-field sill 0.603
max |regression − ordinary| on the grid: 0.894

The fitted slope is negative, about $-2.7$ per unit of dist: log zinc falls as one moves away from the river, which is the floodplain pattern the coordinate plane in Section 10 was only approximating. The residual sill is about 0.27, less than half the sill of the raw field, so distance to the river is carrying a large share of the spatial structure. The residual standard-deviation map is still low near the samples and higher in the gaps — that part is ordinary kriging, unchanged in character, now applied to a smaller residual field. On the grid the regression-kriging surface departs from ordinary kriging by up to about 0.9 log units. That is the covariate doing work in places where there is no nearby zinc sample but there is a known distance to the river.

In-sample, regression kriging matches the data exactly, because ordinary kriging of the residuals is exact and the residuals were computed at those same points. Exactness here is not a measure of predictive skill. Chapter 7.5's cross-validation is the right check.

---

13. Indicator Kriging and the E-type Estimate

Sections 5–12 estimate a mean. Indicator kriging estimates a probability. At a cutoff $t$, each sample is coded $I(Z \ge t)$ and that indicator is ordinary-kriged. The result is read as $P(Z \ge t)$ and clipped to $[0, 1]$. The variogram is the variogram of the indicator, whose sill is at most $p(1-p) \le 1/4$, not the variogram of zinc. Each cutoff is kriged on its own, so a set of probabilities is not forced to decrease as $t$ increases.

Three cutoffs, the quartiles of zinc in ppm, are enough to draw a coarse conditional distribution. For a non-negative variable the mean is the integral of the survival function,

$$ \mathrm{E}[Z] = \int_0^\infty P(Z \ge z)\, dz. $$

IndicatorKriging.etype integrates from z_min with the trapezoid rule, assuming probability 1 at z_min and probability 0 at z_max. Zinc cannot fall below the smallest sample, so that call returns $\mathrm{E}[Z] - z_{\min}$, and $z_{\min}$ is added back. The search uses the 12 nearest samples, the usual local indicator-kriging neighbourhood.

cuts = [float(c) for c in np.quantile(zinc, [0.25, 0.5, 0.75])]
ik = IndicatorKriging(variogram="spherical", n_neighbors=12, n_lags=12).fit(coords, zinc)
probs = ik.predict(grid_xy, cuts)
z_min, z_max = float(zinc.min()), float(zinc.max())
etype = ik.etype(grid_xy, cuts, z_min=z_min, z_max=z_max) + z_min

hi = cuts[-1]
vg_hi = ik.variograms_[hi]
print(f"cutoffs (ppm): {', '.join(f'{c:.0f}' for c in cuts)}")
print(f"indicator variogram at {hi:.0f} ppm: nugget={vg_hi.nugget:.3f}, "
      f"psill={vg_hi.psill:.3f}, range={vg_hi.range:,.0f} m")
print(f"P(zinc ≥ {hi:.0f}) on the grid: {probs[hi].min():.2f} to {probs[hi].max():.2f}, "
      f"mean {probs[hi].mean():.2f}")
print(f"E-type median on the grid: {np.median(etype):.0f} ppm")
print(f"median of exp(ordinary kriging): {np.median(np.exp(pred_grid)):.0f} ppm")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
im0 = axes[0].imshow(to_raster(probs[hi]), extent=EXTENT, origin="lower", cmap="cividis", vmin=0, vmax=1)
axes[0].scatter(coords[:, 0], coords[:, 1], s=10, color="white", edgecolor="k", linewidth=0.3)
axes[0].set_title(f"P(zinc ≥ {hi:.0f} ppm)")
axes[0].set_aspect("equal")
plt.colorbar(im0, ax=axes[0], fraction=0.046, label="probability")

im1 = axes[1].imshow(to_raster(etype), extent=EXTENT, origin="lower", cmap="viridis")
axes[1].scatter(coords[:, 0], coords[:, 1], s=10, color="white", edgecolor="k", linewidth=0.3)
axes[1].set_title("E-type estimate")
axes[1].set_aspect("equal")
plt.colorbar(im1, ax=axes[1], fraction=0.046, label="zinc (ppm)")
plt.tight_layout()
plt.savefig("indicator_kriging.png", dpi=100, bbox_inches="tight")
plt.show()

cutoffs (ppm): 198, 326, 674
indicator variogram at 674 ppm: nugget=0.086, psill=0.115, range=745 m
P(zinc ≥ 674) on the grid: 0.00 to 0.96, mean 0.19
E-type median on the grid: 374 ppm
median of exp(ordinary kriging): 263 ppm

The upper-quartile probability map is a different object from every surface above it. It is high where the local samples themselves sit above about 674 ppm, and it is not a rescaling of the ordinary-kriging mean. The E-type map is back on the ppm scale. Its grid median is several hundred ppm, while the median of $\exp(\hat Z)$, the ordinary-kriging surface exponentiated, is lower. Those two numbers answer different questions. Exponentiating a kriging of log zinc estimates the median of a lognormal field. The E-type integral estimates the mean of the conditional distribution built from the indicator probabilities, tails included. With only three cutoffs the distribution is coarse — a linear interpolation between quartiles, not a full conditional law — but it is the mean on the scale the measurements were made in, which $\exp(\hat Z)$ is not.

---

14. Summary and Conclusions

This chapter introduced the semivariogram, fit it to log zinc on the Meuse floodplain, kriged onto meuse_grid.csv, checked ordinary kriging two different ways, and then changed one assumption at a time: the mean, the extra data, and the target of the prediction.

What we did

1. Explained the semivariogram, its nugget/partial-sill/range parameters, and three common models.
2. Auto-fit a variogram on the 155 samples — spherical, nugget $\approx$ 0.033, partial sill $\approx$ 0.57, range $\approx$ 810 m.
3. Fit ordinary kriging, confirmed it is exact at the samples, and mapped both its prediction and its prediction standard deviation on the 40 m floodplain grid. The standard deviation is low only near samples.
4. Confirmed the kriging variance is genuinely calibrated: standardized leave-one-out errors have mean near 0, standard deviation near 1, and about 95% coverage at the nominal 95% level.
5. Confirmed agreement with pykrige to machine precision at floodplain grid cells, after converting this library's partial sill to pykrige's total-sill convention.
6. Simple-kriged with the sample mean as the known mean. On this grid the surface matches ordinary kriging; the two methods part company where covariances go to zero and only simple kriging must return $m$.
7. Universal-kriged with a degree-1 drift, after fitting the variogram to the trend residuals. Degree 0 reproduces ordinary kriging. The plane changes the map by up to about half a log unit.
8. Cokriged log zinc with log lead under a shared-range linear model of coregionalization. The predictor is exact at the zinc samples. Because lead is colocated with zinc, the grid barely moves.
9. Regression-kriged on distance to the river. The residual sill falls by more than half, and the surface departs from ordinary kriging in the gaps, where dist is known and zinc is not.
10. Indicator-kriged zinc at its three quartiles and integrated the probabilities into an E-type mean in ppm. That mean is not $\exp$ of the ordinary-kriging prediction.

Practical takeaways

• Always check nugget vs. partial sill: a nugget-dominated variogram means there is little spatially structured signal left to krige, whatever the fitted range says.
• Never trust a prediction variance without a calibration check like Section 7's; a wrong variogram produces a confidently wrong uncertainty estimate, not merely an inaccurate mean prediction.
• Remember pykrige's total-sill convention (sill = nugget + partial sill) when cross-checking against it or any other kriging package.
• Choose the kriging by the assumption you are willing to defend. A known mean is simple kriging. A trend in the coordinates is universal kriging, with the variogram fitted to the residuals, and with no faith placed in the polynomial outside the data. A second variable measured on a denser network is cokriging. A covariate known at the prediction locations is regression kriging. A probability, or a mean on the original scale of a skewed variable, is indicator kriging and the E-type integral.
• Do not read exactness at the samples as predictive accuracy. Ordinary, simple, universal and regression kriging are all exact there. Held-out error is a different question, taken up in Chapter 7.5.

Next steps

Chapter 7.5 (Model Validation) puts every method from Chapters 7.3–7.4 through the same cross-validation protocol this chapter previewed with leave-one-out — including the spatial-block variant that, unlike LOO, does not flatter clustered samples.

---

15. Resources

Textbooks

Reference	Notes
Cressie, N. (1993). Statistics for Spatial Data (rev. ed.). Wiley.	The Cressie-weighted least squares fitting used by fit_variogram
Goovaerts, P. (1997). Geostatistics for Natural Resources Evaluation. Oxford University Press.	Ordinary, simple, universal, co- and indicator kriging, and the E-type estimate
Matheron, G. (1963). Principles of geostatistics. Economic Geology, 58, 1246–1266.	The original semivariogram formulation
Journel, A. G. (1983). Nonparametric estimation of spatial distributions. Mathematical Geology, 15, 445–468.	Indicator kriging as a route to the conditional distribution, not only the mean

Key papers

Paper	Contribution
Cressie, N. (1985). Fitting variogram models by weighted least squares. Journal of the International Association for Mathematical Geology, 17, 563–586.	The specific weighting scheme this library's fit_variogram uses

In this library

Object	Role
empirical_variogram, fit_variogram, fit_variogram_auto, VariogramModel	Sections 3–4
OrdinaryKriging	Sections 5–8
SimpleKriging	Section 9
TrendSurface, UniversalKriging	Section 10
empirical_cross_variogram, fit_lmc, Cokriging	Section 11
RegressionKriging	Section 12
IndicatorKriging, IndicatorKriging.etype	Section 13
Chapter 7.5	Extends this chapter's LOO calibration idea to a full cross-validation comparison across every method