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 pillar, and it closes the loop on any interpolation you have already run.

Environment #

pip install pykrige==1.7.2 numpy==1.26.4 scipy==1.13.0 matplotlib==3.9.0

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.

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)")

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.

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.

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.

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.

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 $\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.

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 $(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.

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)")

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 $\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 #

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

• Uncertainty & Variance Mapping — parent guide on predictive uncertainty
• Step-by-Step Ordinary Kriging with PyKrige — where the variance array is produced
• Regression Kriging — combining residual variance with trend uncertainty

← Back to Uncertainty & Variance Mapping