How to Fit a Variogram and Krige a Surface in Python
Problem statement
Kriging is inverse distance weighting with the weights derived from measured spatial structure instead of from a typed-in exponent. That difference is worth about 25% accuracy on a well-behaved variable — measured on 500 elevation samples against 12,000 held-out cells:
kriging 40.2 m RMSE
TIN, linear 46.7 m
IDW, power 2 54.0 m
nearest neighbour 65.1 m
It also returns a variance for every prediction, which is the reason most people reach for it. That variance needs careful reading. On the same run:
mean kriging variance 4,536 m² -> predicted sd 67.3 m
actual RMSE 40.2 m
The kriging variance overstated the real error by 68%. It is a map of where the sample geometry is weak, not a confidence interval.
Quick answer
Ordinary kriging in about forty lines, with no dependency beyond NumPy and SciPy:
import numpy as np
from scipy.spatial import cKDTree
def spherical(h, nugget, sill, rng):
h = np.asarray(h, dtype=float)
inside = nugget + sill * (1.5 * h / rng - 0.5 * (h / rng) ** 3)
return np.where(h == 0, 0.0, np.where(h <= rng, inside, nugget + sill))
def ordinary_kriging(points, values, targets, params, k=24):
"""Predictions and variances from a fitted variogram model."""
tree = cKDTree(points)
_, neighbours = tree.query(targets, k=min(k, len(points)))
neighbours = neighbours.reshape(len(targets), -1)
k = neighbours.shape[1]
predicted = np.empty(len(targets))
variance = np.empty(len(targets))
cache = {}
for i, idx in enumerate(neighbours):
key = idx.tobytes()
if key not in cache:
local = points[idx]
distances = np.sqrt(((local[:, None, :] - local[None, :, :]) ** 2).sum(-1))
A = np.ones((k + 1, k + 1))
A[:k, :k] = spherical(distances, *params)
A[k, k] = 0.0 # Lagrange multiplier row
cache[key] = np.linalg.pinv(A)
b = np.ones(k + 1)
b[:k] = spherical(np.sqrt(((points[idx] - targets[i]) ** 2).sum(-1)), *params)
w = cache[key] @ b
predicted[i] = w[:k] @ values[idx]
variance[i] = w @ b
return predicted, variance
Step-by-step solution
1. Fit the variogram, and check the fit before using it
Kriging is only as good as the variogram. Fit a model, then look at three things: does the empirical curve reach a plateau inside your largest lag, does the fitted range sit inside it, and is the nugget a plausible fraction of the sill?
nugget 0 m² (0.0% of the total)
sill 51,039 m² data variance 41,950 m²
range 4,537 m largest observed lag 4,875 m
The range is inside the data, which is what you want. The sill exceeds the data variance by 22% — expected here, because the range is nearly half the study extent. See The variogram explained.
2. Use a local neighbourhood, not every sample
Global kriging inverts an n × n matrix. At 500 points that is fine; at 50,000 it is impossible, and it is also unnecessary — samples beyond the range contribute weights near zero.
Take the nearest 16 to 32. The 24 used here gave 40.2 m and ran 12,000 predictions in 1.3 seconds.
3. Understand the Lagrange row
The (k+1) × (k+1) system with a row and column of ones and a zero corner is not padding. It is the constraint that the weights sum to one, which is what makes ordinary kriging unbiased when the mean is unknown.
Drop it and you get simple kriging, which needs the true mean supplied and is biased if you supply the wrong one.
4. Do not count on caching the matrix inverse
The left-hand matrix depends only on the neighbour set, not on the target, so it is tempting to cache it — adjacent cells "obviously" share neighbours.
Measured on the full 30 m grid, 116,137 targets against 500 scattered samples with k=24, that assumption is nearly worthless:
116,137 targets -> 113,543 distinct neighbour sets (2.2% cache hits), 12.8 s
With 24 neighbours drawn from a scattered sample, the set changes almost every cell. Keep the cache — it is three lines and it does help when k is small or the samples are gridded — but the thing that actually makes kriging tractable is the local neighbourhood, not the cache.
5. Read the variance as geometry, not as accuracy
mean kriging variance 4,536 m² -> sd 67.3 m
actual RMSE 40.2 m
Kriging variance depends only on the sample positions and the variogram model. The measured values never enter it. Two datasets with identical sample layouts and completely different values produce identical variance maps.
It is therefore excellent for comparing across a map — this corner is better supported than that one — and unreliable as an absolute interval. Here it was 68% too wide, inherited partly from a sill 22% above the data variance.
Code examples
Example 1 — the whole pipeline, variogram to raster
import numpy as np
from scipy.optimize import curve_fit
from scipy.spatial import cKDTree
def krige_surface(points, values, targets, bin_width=250, max_lag=None, k=24):
"""Empirical variogram, spherical fit, then ordinary kriging."""
# ---- empirical variogram -------------------------------------------
d = np.sqrt(((points[:, None, :] - points[None, :, :]) ** 2).sum(-1))
gamma = 0.5 * (values[:, None] - values[None, :]) ** 2
iu = np.triu_indices(len(points), 1)
d, gamma = d[iu], gamma[iu]
max_lag = max_lag or float(np.percentile(d, 60))
edges = np.arange(0, max_lag + bin_width, bin_width)
lags, semis = [], []
for lo, hi in zip(edges[:-1], edges[1:]):
m = (d >= lo) & (d < hi)
if m.sum() >= 30:
lags.append((lo + hi) / 2)
semis.append(gamma[m].mean())
lags, semis = np.array(lags), np.array(semis)
# ---- fit -----------------------------------------------------------
params, _ = curve_fit(
spherical, lags, semis,
p0=[semis[0] * 0.1, semis.max(), lags.max() / 2],
bounds=([0, 0, lags[0]], [semis.max(), semis.max() * 3, lags.max() * 3]),
maxfev=20000,
)
nugget, sill, rng = params
print(f" nugget {nugget:,.0f} sill {sill:,.0f} range {rng:,.0f} m")
if rng > lags.max():
print(" ! range beyond the largest lag — treat it as a lower bound")
# ---- krige ---------------------------------------------------------
predicted, variance = ordinary_kriging(points, values, targets, params, k=k)
print(f" predicted {predicted.min():.1f} to {predicted.max():.1f}, "
f"mean sd {np.sqrt(variance.mean()):.1f}")
return predicted, variance, params nugget 0 sill 51,039 range 4,537 m
predicted 51.6 to 987.8, mean sd 67.2
Note that the surface reaches below the lowest sample (61.2 m). That is not a bug — see the explanation of negative weights below — but it is worth knowing before you assume a kriged surface is bounded by its inputs.
Example 2 — kriging in blocks, so it scales
import numpy as np
from scipy.spatial import cKDTree
def krige_blocks(points, values, targets, params, k=24, block=4096):
"""Krige a large target set in chunks, sharing the neighbour cache."""
tree = cKDTree(points)
predicted = np.empty(len(targets))
variance = np.empty(len(targets))
cache = {}
for start in range(0, len(targets), block):
chunk = targets[start:start + block]
_, neighbours = tree.query(chunk, k=min(k, len(points)))
neighbours = neighbours.reshape(len(chunk), -1)
for j, idx in enumerate(neighbours):
key = idx.tobytes()
if key not in cache:
local = points[idx]
dd = np.sqrt(((local[:, None, :] - local[None, :, :]) ** 2).sum(-1))
A = np.ones((len(idx) + 1, len(idx) + 1))
A[:len(idx), :len(idx)] = spherical(dd, *params)
A[len(idx), len(idx)] = 0.0
cache[key] = np.linalg.pinv(A)
b = np.ones(len(idx) + 1)
b[:len(idx)] = spherical(
np.sqrt(((points[idx] - chunk[j]) ** 2).sum(-1)), *params)
w = cache[key] @ b
predicted[start + j] = w[:len(idx)] @ values[idx]
variance[start + j] = w @ b
print(f" {start + len(chunk):,}/{len(targets):,} "
f"({len(cache):,} distinct neighbour sets cached)")
return predicted, variancePrint the cache size and be honest about what it buys. On the grid measured here it was 2.2% — the neighbour set changed almost every cell. Chunking still matters for memory, and the cache pays off when k is small or the samples are themselves on a grid.
Example 3 — validating the variogram by its own cross-validation
import numpy as np
def krige_diagnostics(points, values, params, k=24):
"""Leave-one-out kriging, checking both the predictions and the variances."""
n = len(points)
predicted = np.empty(n)
variance = np.empty(n)
for i in range(n):
keep = np.ones(n, bool); keep[i] = False
p, v = ordinary_kriging(points[keep], values[keep],
points[i:i + 1], params, k=k)
predicted[i], variance[i] = p[0], v[0]
residual = values - predicted
rmse = float(np.sqrt(np.mean(residual ** 2)))
# standardised residuals should have a variance of about 1 if the
# variogram is right — this is the real test of the model
standardised = residual / np.sqrt(np.maximum(variance, 1e-12))
print(f" LOO RMSE {rmse:8.2f}")
print(f" mean standardised residual {standardised.mean():8.3f} (want ~0)")
print(f" variance of standardised {standardised.var():8.3f} (want ~1)")
if standardised.var() < 0.7:
print(" ! variances are too large — the sill is probably overestimated")
elif standardised.var() > 1.4:
print(" ! variances are too small — the model is overconfident")
return {"rmse": rmse, "std_var": float(standardised.var())}The variance of the standardised residuals is the diagnostic almost nobody runs, and it is the one that tests the variogram rather than the interpolation. Below 1 means the model's uncertainty is inflated — which is exactly the 68% overstatement measured at the top of this page.
Explanation
Why kriging beats IDW
Two mechanisms, both absent from IDW.
Declustering. The left-hand matrix holds the variogram between samples. Two samples close together have low semivariance between them, and the solved weights account for that redundancy — the pair jointly gets roughly the weight one of them would get alone. IDW gives both full weight, so a cluster of samples outvotes an isolated one that carries more information.
Structure-derived decay. The rate at which influence falls off comes from the fitted variogram, not from a chosen exponent. If the field is smooth over kilometres, kriging spreads weight widely; if it decorrelates in 200 m, it concentrates.
Why the variance is not a confidence interval
Look at the kriging system: A holds variogram values between samples, b holds variogram values from samples to the target. Neither contains a single measured value. The variance w·b therefore depends only on where the samples are and what the variogram model says — not on what was measured there.
That has a useful consequence: you can compute the variance map before collecting the data, and use it to design the survey. It has an awkward one: it cannot detect that the model is wrong, and it inherits every error in the variogram.
Here the sill was 22% above the data variance, and the variance scales directly with the sill. That plus the smoothing effect gave a predicted standard deviation of 67.3 m against a real RMSE of 40.2 m.
Why the surface is smoother than the data — but not bounded by it
Kriging minimises expected squared error, and the estimator that minimises squared error is a conditional mean. Means are smoother than the things they average, so a kriged surface systematically under-represents extremes. Measured over the full grid: the kriged surface had a standard deviation of 193.8 m against the DEM's 204.8 m.
What surprises people is that the surface is nevertheless not bounded by the samples. It reached 51.6 m where the lowest sample was 61.2 m.
The reason is that kriging weights can be negative. On the grid measured here, every single target — 116,137 of 116,137 — had at least one negative weight. This is the screen effect: when a sample sits directly behind a nearer one along the same bearing, the system gives the far one a small negative weight, because the near sample already carries its information and the far one now mostly contributes redundancy.
Negative weights are what let kriging extrapolate a local gradient slightly beyond the data. They are also why a kriged concentration or rainfall surface can come out negative, which is physically impossible. If that matters, either constrain the weights to be non-negative, or transform the variable (log, logit) before kriging and back-transform after.
If you need a surface with realistic variability rather than the smoothest one consistent with the data, use conditional simulation.
When kriging is not worth it
- No spatial structure. A flat variogram means kriging degenerates to the mean.
- A poor variogram fit. Garbage in; IDW with sensible parameters beats kriging with a bad model.
- Very large sample sets without a local neighbourhood, where the matrix algebra dominates.
- Categorical variables. Indicator kriging exists but is a different method.
- A strong trend. Ordinary kriging assumes a locally constant mean. With a strong regional trend, remove it first and krige the residuals.
Edge cases or notes
- Duplicate coordinates make the matrix singular. Aggregate them, or use
pinvas above. np.linalg.pinvis a safety net, not a fix. If it is doing real work, the neighbourhood has duplicates.- Project the coordinates. Kriging weights by distance.
- A fitted range beyond your largest lag is a lower bound, not a measurement.
- Standardised residual variance well below 1 means inflated variances — usually an overestimated sill.
- Kriged surfaces are smoother than reality (sd 193.8 against the DEM's 204.8 here). Do not use them where variance matters; simulate instead.
- Kriging weights can be negative — all 116,137 targets had at least one here — so the surface can leave the range of the data and a non-negative variable can come out negative.
- Use 16–32 neighbours. More costs time and changes little.
- The nugget is a discontinuity at zero, which is why
spherical()special-casesh == 0. Omit that and the surface no longer honours the data.
Internal links
- The variogram explained — fitting and checking the model kriging needs
- My kriging fit fails or returns a flat surface — when the fit goes wrong
- IDW, kriging, splines or TIN? Choosing an interpolator — whether the 25% is worth it
- Spatial interpolation explained — the error-versus-distance picture
- How to cross-validate an interpolated surface — validating properly
- How to interpolate points to a grid with IDW in Python — the simpler baseline
- Sample design explained: where to measure — designing from the variance map
- Spatial autocorrelation explained — the structure kriging exploits
FAQ
How do I do kriging in Python?
Fit a variogram model to the empirical variogram, then solve the ordinary kriging system for each target using its nearest 16–32 samples. It is about forty lines with NumPy and SciPy.
Is kriging more accurate than IDW?
It was 25% better here — 40.2 m against 54.0 m on the same 500 samples and 12,000 held-out cells. The gain depends on having a well-fitted variogram.
What does the kriging variance mean?
How poorly the target is surrounded by samples, given the variogram. It never sees the measured values, so it is a geometry map rather than an error estimate — here it overstated the true error by 68%.
Why is my kriged surface flatter than my data?
Because kriging is a conditional-mean estimator, and means are smooth. The kriged surface here had a standard deviation of 193.8 m against the DEM's 204.8 m. Use conditional simulation if you need realistic variability.
Why did my kriged surface produce a negative value?
Kriging weights can be negative — the screen effect — so the surface is not bounded by the data. Every target on the grid measured here had at least one negative weight. Constrain the weights or transform the variable if negatives are impossible.
How many neighbours should I use?
16 to 32. Global kriging inverts an n × n matrix for no benefit, because samples beyond the range get near-zero weight anyway.
Why is my kriging matrix singular?
Almost always duplicate sample coordinates. Aggregate them before kriging.
How do I know the variogram model is right?
Compute leave-one-out standardised residuals — the residual divided by the kriging standard deviation. Their variance should be about 1. Well below 1 means the variances are inflated.