Purpose: interpolation that reports what it doesn't know. The difference between a professional product and a pretty raster is the uncertainty surface and an honest cross-validation — both are non-optional here.
Method selection
| Situation | Method |
|---|---|
| Dense, smooth phenomenon, quick look | IDW (report power parameter; test 1-3) |
| Physical phenomenon with spatial structure, need uncertainty | Ordinary kriging (default professional choice) |
| Clear trend (elevation gradient, coastal effect) | Universal kriging or regression kriging on covariates |
| Strong covariates available (DEM, land cover, distances) | Regression kriging / random-forest residual kriging |
| Categorical target | Indicator kriging |
| Honeycomb-free tessellation, no extrapolation wanted | Natural neighbor |
IDW is a reasonable baseline but has no error model and bullseyes around extremes; say so when delivering IDW-only products.
Exploratory phase (before any interpolation)
- Map the points with values; look for duplicates at identical coordinates (average or offset them — kriging matrices go singular otherwise).
- Histogram + skew: strongly skewed variables (rainfall, concentrations) usually want a log/normal-score transform; back-transform predictions properly (bias correction for lognormal kriging).
- Trend check: regress value on x, y, and candidate covariates; visible trend → universal/regression kriging path.
- Declustering if sampling is preferential (dense where values are high).
Variogram discipline
The variogram is a MODELING decision, not an auto-fit output:
import gstools as gs
bin_center, gamma = gs.vario_estimate((x, y), values, max_dist=dmax) # dmax ≈ half extent
model = gs.Exponential(dim=2)
model.fit_variogram(bin_center, gamma, nugget=True)
print(model) # report: nugget, sill, range — plus the fitted plot
- Max lag ≈ half the domain diameter; ≥30 pairs per bin.
- Check anisotropy with directional variograms (0/45/90/135°); geological and meteorological fields are often anisotropic — fit an anisotropic model rather than ignoring it.
- Interpret and report the parameters in words: nugget (measurement error + micro-scale variance), range (correlation distance), sill. A nugget near the sill means the data barely support interpolation — say that honestly.
- Never interpolate meaningfully beyond the variogram range from the nearest sample; mask or flag those cells.
Kriging execution
from pykrige.ok import OrdinaryKriging
ok = OrdinaryKriging(x, y, values, variogram_model="exponential",
variogram_parameters={"sill": s, "range": r, "nugget": n},
coordinates_type="euclidean")
z, ss = ok.execute("grid", gridx, gridy) # ss = kriging VARIANCE — keep it!
- Work in a projected CRS (kriging distances in degrees are wrong except with explicitly geographic-aware models).
- Deliver TWO rasters: prediction AND kriging standard deviation
(
np.sqrt(ss)). The uncertainty map drives where to sample next and where not to trust the map. - Search neighborhood: 16-32 nearest points typical; document it.
Validation — leave-one-out and beyond
- LOOCV for small n: report ME (bias ≈ 0?), RMSE, and standardized RMSE (should be ≈ 1 if kriging variances are honest).
- For clustered samples, LOOCV flatters; also run spatial block CV
(canonical recipe:
ml-experiment-standards→references/spatial-cv-protocol.md) and report both. - Scatter plot observed vs predicted with 1:1 line; map the CV residuals — spatially clustered residuals mean missing trend/covariate.
- Compare against the dumb baseline (global mean, IDW): kriging must earn its complexity.
Reporting template
## Surface: <variable>
- n = <>, extent, CRS, transform applied: <log/none>
- Variogram: <model>, nugget/sill/range = <...>, anisotropy: <...>
- Method: <OK/UK/RK + covariates>
- LOOCV: ME <>, RMSE <>, RMSE_std <>; block-CV RMSE <>
- Uncertainty: kriging SD raster delivered; area beyond reliable range masked
Pitfalls checklist
- Auto-fitted variogram accepted without looking at the plot.
- Kriging in EPSG:4326 over large extents.
- Prediction raster delivered without its variance/SD companion.
- Log-kriging back-transformed by naive
exp()(bias!). - Extrapolation far beyond the data hull presented at equal confidence.
- Duplicated station coordinates crashing or silently distorting the fit.
Execution contract
- Workflow: inspect sampling and transform needs; model spatial dependence; fit candidate surfaces; validate against simple baselines; quantify uncertainty; mask unsupported extrapolation.
- Decision rules: choose IDW or trend methods for transparent baselines, kriging when a defensible variogram exists, and regression kriging when covariates improve blocked validation.
- Verification protocol: examine variogram diagnostics, LOOCV and spatial block-CV residuals, baseline uplift, uncertainty calibration, and behavior beyond the sample hull.
- Failure modes: withhold a confident surface when sampling is too sparse, dependence is absent, duplicates dominate, distance units are invalid, or uncertainty is uncalibrated.
- Deliverables: prediction surface, uncertainty surface, variogram and parameters, validation table, sampling-support mask, and method limitations.
- Source freshness: consult the authoritative source registry before relying on package APIs or defaults and record the checked date.