Geostatistical Modeling Spatial Statistics with Python

Mapping Kriging Variance Surfaces in Python

Every kriging call returns uncertainty for free: z_pred, ss = ok.execute('grid', gridx, gridy). The second array ss is the kriging variance in squared target units. Take sigma = np.sqrt(ss) for the standard error, build 95% prediction intervals as z_pred ± 1.96 * sigma, and map sigma with a sequential colormap. High values mark data-sparse, low-confidence zones.

Why This Matters

A kriging prediction without its variance is only half an answer. The kriging variance surface tells you where to trust the map — it pinpoints the data-sparse and extrapolation zones where predictions are least reliable, and it turns a single best-guess surface into calibrated prediction intervals for risk-aware decisions. This guide is the practical core of the Uncertainty & Variance Mapping workflow inside the broader Kriging, Interpolation & Surface Generation Techniques section, and it closes the loop on any interpolation you have already run.

Environment

bash
pip install pykrige==1.7.2 numpy==1.26.4 scipy==1.13.0 matplotlib==3.9.0
python
import numpy as np                          # 1.26.4
from pykrige.ok import OrdinaryKriging      # pykrige 1.7.2
from scipy.stats import norm                # scipy 1.13.0
import matplotlib.pyplot as plt             # 3.9.0

Python 3.10 or 3.11 recommended. Pin numpy<2 if your pykrige build predates NumPy 2 support.

Step-by-Step Implementation

Step 1 — Run Kriging and Capture the Variance Array

ok.execute returns a tuple. The first element is the predicted surface; the second, conventionally named ss, is the kriging variance at every grid cell. Keep it — this is the whole point.

python
rng = np.random.default_rng(3)
n = 80
x = rng.uniform(0, 10000, n)
y = rng.uniform(0, 10000, n)
z = 12 + 3 * np.sin(x / 2500) + rng.normal(0, 0.4, n)

gridx = np.linspace(0, 10000, 160)
gridy = np.linspace(0, 10000, 160)

ok = OrdinaryKriging(x, y, z, variogram_model="spherical",
                     nlags=15, weight=True, verbose=False)

z_pred, ss = ok.execute("grid", gridx, gridy)   # ss = kriging variance
print(f"Variance range: {ss.min():.4f} to {ss.max():.4f}  (units: z^2)")
What each half of the returned tuple depends on Four horizontal bands list the inputs to ordinary kriging: sample coordinates, observed values z, the variogram model, and the prediction grid. Two columns to the right show whether each input feeds z_pred and whether it feeds ss. Every input feeds both arrays except the observed values z, which feed z_pred but are marked with a red cross under ss. z_pred, ss = ok.execute('grid', gridx, gridy) — what feeds which array z_pred — prediction (160×160) ss — variance, units z² (160×160) Sample coordinates x, y  (n = 80) sets the kriging weights geometry is most of the answer Observed values z at those points sets every predicted value never enters the ss equations Variogram model (spherical, nlags = 15) supplies the covariances supplies the sill and the range Prediction grid, 160 × 160 cells where the surface is evaluated where the variance is evaluated Shuffle every z value and ss comes back identical: the variance array is built from the geometry and the variogram alone.

Step 2 — Convert Variance to Standard Error

The variance is in squared units of the target. The standard error, its square root, is in the target’s own units and is what prediction intervals need. Clip tiny negatives that can arise from floating-point rounding at sample points.

python
ss = np.asarray(ss)
ss = np.clip(ss, 0, None)          # guard against -1e-15 rounding artefacts
sigma = np.sqrt(ss)                # kriging standard error, units of z
print(f"Std error range: {sigma.min():.4f} to {sigma.max():.4f}  (units: z)")

Step 3 — Build 95% Prediction Intervals

Under an approximate Gaussian assumption, the 95% prediction interval is the prediction plus or minus 1.96 standard errors. The interval width is a direct uncertainty map.

python
z_score = norm.ppf(0.975)          # 1.9599...
lower = z_pred - z_score * sigma
upper = z_pred + z_score * sigma
pi_width = upper - lower           # = 2 * 1.96 * sigma

print(f"Mean 95% PI width: {pi_width.mean():.3f}")
print(f"Widest 95% PI    : {pi_width.max():.3f} (most uncertain cell)")

Step 4 — Map the Variance Surface

Plot the standard error (or variance) with a sequential colormap and overlay the samples. The visual should show troughs of low uncertainty hugging the data and ridges of high uncertainty in the gaps.

python
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

im0 = axes[0].pcolormesh(gridx, gridy, z_pred, cmap="viridis", shading="auto")
axes[0].scatter(x, y, c="white", edgecolors="black", s=25, linewidth=0.4, zorder=5)
fig.colorbar(im0, ax=axes[0], label="Prediction (z)")
axes[0].set_title("Kriging prediction")

im1 = axes[1].pcolormesh(gridx, gridy, sigma, cmap="magma", shading="auto")
axes[1].scatter(x, y, c="cyan", edgecolors="black", s=25, linewidth=0.4, zorder=5)
fig.colorbar(im1, ax=axes[1], label="Std error (z)")
axes[1].set_title("Kriging standard error — uncertainty map")

for ax in axes:
    ax.set_xlabel("Easting (m)"); ax.set_ylabel("Northing (m)")
plt.tight_layout(); plt.savefig("kriging_variance.png", dpi=150); plt.show()

Step 5 — Flag High-Variance Zones

Threshold the standard error to delineate where the surface should not be trusted — for example, cells above the 90th percentile of standard error.

python
thr = np.percentile(sigma, 90)
low_conf = sigma > thr
frac = low_conf.mean()
print(f"Std-error 90th percentile threshold: {thr:.3f}")
print(f"Fraction of grid flagged low-confidence: {frac:.1%}")

Interpreting the Output

  • Kriging variance (ss): near zero at and around sampled points, rising with distance from data toward the sill in empty regions. It depends only on sample geometry and the variogram — never on the observed values — so it is a structural map of information density, not of prediction error per se.
  • Standard error (sigma): the variance in interpretable units. A standard error of 0.5 on a target measured in the same units means roughly a ±1.0\pm 1.0 swing at 95% confidence.
  • Prediction-interval width: the actionable uncertainty. Wide intervals mark places where a decision based on the point prediction carries real risk; narrow intervals mark well-constrained areas.
  • High-variance zones: these are recommendations for where to sample next. Densifying observations in flagged areas is the most efficient way to shrink overall uncertainty.
Kriging standard error along a transect A line chart of kriging standard error sigma against easting from 0 to 10000 metres. The curve drops to about 0.26 at each of eleven sample locations drawn as dots on the axis, stays low where samples are dense between 300 and 4600 metres, and rises to about 1.05 in the middle of the empty stretch from 4600 to 8300 metres. A dashed line at 0.86 marks the 90th-percentile threshold and a shaded strip marks the cells it flags. The curve also lifts at both map edges where kriging extrapolates. Where the uncertainty lives: standard error read along one transect of the grid kriging standard error σ 90th-percentile threshold (Step 5) sample location cells flagged low-confidence flagged: σ above the 90th percentile thr = 0.86 0 0.4 0.8 1.2 std error σ (units of z) 0 2 500 5 000 7 500 10 000 Easting along the transect (m) σ collapses toward the nugget on the data edge rise: kriging is extrapolating Two shapes carry the whole map: troughs pinned to the samples, and the ridge across the 3 700 m gap where kriging has nothing nearby to lean on.

Critical Best Practices

Kriging Variance Is Not Validated Error

The kriging variance is a model-internal quantity. It reflects the variogram and sampling layout, not whether your variogram is correct. Always cross-check it against cross-validated errors: compute standardized errors (zi−z^i)/σi(z_i - \hat{z}_i)/\sigma_i by leave-one-out; their variance should be close to 1 if the variance surface is calibrated. If it is far from 1, the variogram — and therefore every interval — is miscalibrated.

python
from sklearn.model_selection import LeaveOneOut
loo, se = LeaveOneOut(), []
idx = np.arange(n)
for i in idx:
    tr = idx != i
    ok_cv = OrdinaryKriging(x[tr], y[tr], z[tr], variogram_model="spherical",
                            nlags=12, weight=True, verbose=False)
    zp, sp = ok_cv.execute("points", x[i:i+1], y[i:i+1])
    se.append((float(zp[0]) - z[i]) / np.sqrt(max(float(sp[0]), 1e-12)))
print(f"Std. error variance: {np.var(se):.3f}  (target ~1.0)")
Reading the leave-one-out standardized-error variance Three side-by-side panels each plot the distribution of leave-one-out standardized errors as a solid curve over a dashed standard normal. In the left panel the variance is 0.42 and the curve is tall and narrow, so sigma is overstated and the nominal 95 percent interval really covers 99.7 percent. In the middle panel the variance is 1.03 and the two curves coincide, giving about 95 percent coverage. In the right panel the variance is 2.70 and the curve is low and wide, so sigma is understated and coverage falls to about 77 percent. Is the variance surface calibrated? Read the spread of the leave-one-out standardized errors Solid: distribution of (z − prediction)/σ per held-out sample. Dashed: N(0, 1), the shape a calibrated variance surface must produce. np.var(se) = 0.42 np.var(se) = 1.03 np.var(se) = 2.70 −2 0 2 −2 0 2 −2 0 2 σ overstated 95% intervals really cover ~99.7% safe ground gets flagged as risky cause: sill fitted too high calibrated ±1.96σ covers ~95%, as advertised the variance map can be published cause: the variogram fits σ understated 95% intervals really cover ~77% sparse zones look safe — they are not cause: nugget or sill fitted too low The kriging variance is computed from the variogram, so a wrong variogram still returns a smooth, confident, wrong uncertainty surface. Accept roughly 0.8–1.25 from the loop above; outside that band, re-fit the variogram before you publish an interval.

Square-Root Before You Interval

Prediction intervals use the standard error, not the variance. Forgetting the square root inflates interval widths by orders of magnitude. Convert with np.sqrt(ss) and clip negatives from floating-point noise first.

The Gaussian Interval Assumes Normal Errors

The ±1.96σ\pm 1.96\sigma interval is exact only if kriging errors are Gaussian. For skewed data (concentrations, rainfall), transform (log or normal-score) before kriging and back-transform the interval, or use quantile-based approaches. Check the standardized-error histogram for symmetry.

Variance Ignores Trend Uncertainty

For a trend-based model, the kriging variance of the residual step does not include the trend model’s own uncertainty. When combining a trend with kriged residuals, as in regression kriging, the residual variance understates total uncertainty; add the trend prediction variance before forming intervals.

Troubleshooting

Symptom Likely cause Fix
Tiny negative variance values Floating-point rounding near samples np.clip(ss, 0, None) before sqrt
Variance uniformly near the sill Range too short relative to spacing Refit the variogram; range should span sample gaps
Intervals absurdly wide Used variance instead of standard error Take sqrt(ss) before multiplying by 1.96
Standardized-error variance far from 1 Miscalibrated variogram (sill/nugget wrong) Re-fit variogram; validate with LOOCV
Variance surface looks flat Nugget dominates Inspect nugget/sill ratio; reconsider model
Skewed target, poor interval coverage Non-Gaussian errors Normal-score or log transform, then back-transform intervals

Next Steps

Return to the Uncertainty & Variance Mapping guide for the broader treatment of predictive uncertainty, and revisit where the variance array comes from in Step-by-Step Ordinary Kriging with PyKrige.

Frequently Asked Questions

What does the ss array from PyKrige execute contain?

ok.execute returns a tuple (z, ss). The second array, ss, is the kriging variance at each prediction location, in squared units of the target variable. It is near zero at sampled points and rises toward the sill in data-sparse areas. Take its square root to get the kriging standard error in the target’s own units, which is what you use for prediction intervals.

How do I build a 95% prediction interval from kriging variance?

Convert variance to standard error with sigma = sqrt(ss), then form the interval as prediction plus or minus 1.96 * sigma under a Gaussian assumption. The interval width 2 * 1.96 * sigma is a direct map of predictive uncertainty. This assumes the kriging errors are approximately normal; validate that assumption with cross-validated standardized errors before trusting the coverage.

Why is kriging variance high in some areas but not others?

Kriging variance depends only on the sampling geometry and the variogram, not on the observed values. It is low where prediction locations are close to many samples and high where samples are sparse or the location sits beyond the sampled area. High-variance zones flag where you should collect more data or treat predictions cautiously, since the interpolator is extrapolating.

Related

← Back to Uncertainty & Variance Mapping

Related

Building Prediction Intervals from Kriging Variance Sequential Gaussian Simulation for Uncertainty in Python