How to Cross-Validate an Interpolated Surface
Problem statement
You interpolate a surface and want a number for how good it is. The standard answer is leave-one-out cross-validation: remove each sample, predict it from the rest, and take the RMSE of the residuals.
That number is not the error of your surface. It is the error of your surface at the sample locations, and the sample locations are the best-supported places on the map.
Measured against a real elevation model, where the true error over every cell is known:
leave-one-out CV RMSE at the 500 sample points 60.8 m
true RMSE over all 116,137 cells 54.0 m
Here LOO was 13% pessimistic. Under clustered sampling the sign flips and it becomes wildly optimistic. Neither direction is safe to assume.
Quick answer
import numpy as np
from scipy.spatial import cKDTree
def loo_cv(points, values, predict, **kwargs):
"""Leave-one-out cross-validation for any interpolator."""
predicted = np.empty(len(points))
for i in range(len(points)):
keep = np.ones(len(points), bool)
keep[i] = False
predicted[i] = predict(points[keep], values[keep],
points[i:i + 1], **kwargs)[0]
residual = values - predicted
return {
"rmse": float(np.sqrt(np.mean(residual ** 2))),
"mae": float(np.mean(np.abs(residual))),
"bias": float(np.mean(residual)),
"r2": float(1 - residual.var() / values.var()),
}Then report it honestly: "leave-one-out RMSE at the sample locations", not "the accuracy of the map".
Step-by-step solution
1. Understand what LOO actually measures
Removing one point from a well-spread sample leaves a hole roughly one sample-spacing wide, and the prediction has to bridge it. That is a harder problem than the average cell faces, because the average cell has samples on all sides.
Hence the 13% pessimism measured here: 60.8 m from LOO against a true 54.0 m.
2. Understand why clustering flips the sign
In a clustered design, removing one point from a cluster leaves its neighbours a few metres away. The prediction is then almost exact, and LOO reports an error close to the measurement noise.
Meanwhile the surface between clusters is terrible. Measured on a 500-point clustered design, the true RMSE was 215.7 m against 58.7 m for the same number of randomly placed points — and the clustered design's LOO would report a far better number than the random one's.
Cross-validation on clustered data is worse than no validation, because it is confidently wrong in the reassuring direction.
3. Use spatial (block) cross-validation instead
Hold out whole regions, not individual points, so the held-out points do not have their own neighbours to lean on:
def block_folds(points, n_blocks=5, seed=0):
"""Assign points to spatial blocks by k-means on their coordinates."""
from sklearn.cluster import KMeans
return KMeans(n_clusters=n_blocks, random_state=seed,
n_init=10).fit_predict(points)Block CV answers "how well does this surface do somewhere I have not sampled", which is the question a map user is asking. It is systematically more pessimistic than LOO, and that pessimism is the point.
4. Report the error as a function of distance, not as one number
A single RMSE hides the thing that matters most. Measured error by distance to the nearest sample:
0 - 100 m 6,959 cells 18.9 m
100 - 250 m 7,286 cells 54.0 m
250 - 500 m 9,479 cells 106.1 m
500 - 1,000 m 19,802 cells 169.8 m
over 1,000 m 72,611 cells 334.3 m
One number would say "about 250 m". The table says the map is excellent in some places and useless in others, and exactly where.
5. Check the residuals for structure
residual = values - predictedResiduals should look like noise. If they are spatially clustered — all positive in one region, all negative in another — there is a trend the interpolator is not capturing, and the fix is to model the trend rather than to tune the interpolator.
Code examples
Example 1 — LOO, k-fold and block CV side by side
import numpy as np
from scipy.spatial import cKDTree
def compare_validation(points, values, predict, n_blocks=5, k_folds=5, seed=0):
"""The same interpolator scored three ways, so the gap is visible."""
rng = np.random.default_rng(seed)
n = len(points)
def score(fold_ids):
predicted = np.empty(n)
for fold in np.unique(fold_ids):
test = fold_ids == fold
predicted[test] = predict(points[~test], values[~test], points[test])
return float(np.sqrt(np.mean((values - predicted) ** 2)))
loo = score(np.arange(n))
kfold = score(rng.integers(0, k_folds, n))
from sklearn.cluster import KMeans
blocks = KMeans(n_clusters=n_blocks, random_state=seed,
n_init=10).fit_predict(points)
block = score(blocks)
print(f" leave-one-out {loo:7.2f}")
print(f" random {k_folds}-fold {kfold:7.2f}")
print(f" {n_blocks} spatial blocks {block:7.2f}")
print(f" block / LOO ratio {block / loo:7.2f}")
return {"loo": loo, "kfold": kfold, "block": block}A block-to-LOO ratio near 1 means the sampling is well spread and LOO is trustworthy. A ratio of 2 or more means the samples are clustered enough that LOO is measuring redundancy, not skill.
Example 2 — error as a function of support distance
import numpy as np
from scipy.spatial import cKDTree
def error_by_distance(points, values, targets, truth, predict,
bands=((0, 100), (100, 250), (250, 500),
(500, 1000), (1000, np.inf))):
"""The table that a single RMSE hides."""
predicted = predict(points, values, targets)
nearest, _ = cKDTree(points).query(targets, k=1)
rows = []
for lo, hi in bands:
band = (nearest >= lo) & (nearest < hi)
if band.sum() < 20:
continue
rmse = float(np.sqrt(np.mean((predicted[band] - truth[band]) ** 2)))
label = f"{lo:>5.0f}-{'inf' if np.isinf(hi) else int(hi):>5}"
print(f" {label} m {int(band.sum()):7,} cells RMSE {rmse:7.2f}")
rows.append({"lo": lo, "hi": None if np.isinf(hi) else hi,
"cells": int(band.sum()), "rmse": rmse})
return rowsThis needs ground truth over the whole area, which you usually do not have — but you can get it for a pilot: take a dense dataset, thin it to a realistic survey, interpolate, and compare. Doing that once for your variable and landscape calibrates your intuition permanently.
Example 3 — a validation report that will not mislead the reader
import numpy as np
from scipy.spatial import cKDTree
def validation_report(points, values, predict, bounds, n_blocks=5):
"""Everything a reader needs to judge the surface, in one dict."""
from sklearn.cluster import KMeans
n = len(points)
left, bottom, right, top = bounds
area_km2 = (right - left) * (top - bottom) / 1e6
# coverage first: it determines everything else
rng = np.random.default_rng(0)
probes = np.column_stack([rng.uniform(left, right, 20000),
rng.uniform(bottom, top, 20000)])
gap, _ = cKDTree(points).query(probes, k=1)
blocks = KMeans(n_clusters=n_blocks, random_state=0, n_init=10).fit_predict(points)
predicted = np.empty(n)
for fold in np.unique(blocks):
test = blocks == fold
predicted[test] = predict(points[~test], values[~test], points[test])
residual = values - predicted
report = {
"n_samples": n,
"area_km2": round(area_km2, 1),
"ideal_spacing_m": round(float(np.sqrt(area_km2 * 1e6 / n))),
"mean_gap_m": round(float(gap.mean())),
"p95_gap_m": round(float(np.percentile(gap, 95))),
"max_gap_m": round(float(gap.max())),
"block_cv_rmse": round(float(np.sqrt(np.mean(residual ** 2))), 2),
"block_cv_bias": round(float(residual.mean()), 2),
"n_blocks": n_blocks,
"caveat": "block CV RMSE is measured at sample locations; error grows "
"with distance from the nearest sample",
}
for key, value in report.items():
print(f" {key:18} {value}")
return reportShipping the gap percentiles beside the RMSE is what stops the RMSE being quoted alone. The caveat string is not padding — it is the sentence that will otherwise be omitted from every summary of your work.
Explanation
Why LOO was pessimistic here and optimistic elsewhere
The removed point's prediction quality depends on what is left nearby.
Under even sampling, removing a point leaves a gap about one spacing wide, which is a slightly harder problem than the average cell faces — hence 60.8 m against a true 54.0 m, 13% pessimistic.
Under clustered sampling, removing a point leaves its cluster-mates metres away, which is a far easier problem than the average cell faces. LOO then reports near-perfect skill for a surface with an RMSE of 215.7 m.
Both are the same mechanism. The bias direction is set by how the sample spacing at the held-out point compares with the typical distance from a map cell to its nearest sample.
Why block CV is the right default
Holding out a spatial block forces predictions across a real gap. The remaining points are not neighbours of the held-out ones, so the score reflects what a map user experiences: standing somewhere nobody sampled.
The usual objection is that block CV is pessimistic because each fold trains on less area. That is true, and it is the safer error to make. A surface whose block CV is acceptable is a surface you can defend.
Why R² is a poor choice here
R² compares your errors against predicting the mean everywhere. On a variable with strong spatial structure, predicting the mean is a very weak baseline, so R² comes out high even for a mediocre surface.
RMSE in the variable's own units is more informative — "±54 m of elevation" is a statement anyone can act on, and "R² = 0.93" is not.
Why validation cannot rescue a bad design
Cross-validation scores the interpolator. It cannot score the survey, because it only ever asks questions at locations the survey chose.
If the samples cover 60% of the area and cluster in the accessible parts, no validation scheme computed from those samples alone will reveal the state of the other 40%. The only instruments that will are the coverage statistics — mean, 95th-percentile and maximum gap — which need no ground truth at all.
Edge cases or notes
- Never let a point predict itself. With a k-d tree, the nearest neighbour of a sample is itself; drop the first column.
- LOO with
nsamples costsnfits. For kriging, refit the weights, not the variogram, or it becomes very slow. - Refit the variogram inside the fold if you want a strictly honest score; fitting it on all the data leaks information.
- Random k-fold behaves like LOO for this purpose — both leave the held-out points with nearby neighbours.
- Report bias as well as RMSE. A surface can have small RMSE and a systematic offset.
- Check residual autocorrelation. Spatially structured residuals mean an unmodelled trend.
- Weight the score by area, not by sample, if the design is clustered — otherwise the dense regions dominate.
- Report the coverage statistics alongside, always.
Internal links
- Spatial interpolation explained — the error-versus-distance table
- Sample design explained: where to measure — why clustering corrupts validation
- IDW, kriging, splines or TIN? Choosing an interpolator — running the bake-off honestly
- Spatial leakage explained: why random cross-validation lies — the same problem in machine learning
- How to do spatial cross-validation in Python — block folds in a modelling pipeline
- How to fit a variogram and krige a surface in Python — validating the variogram, not just the surface
- My interpolated surface extends far beyond the data — the region validation cannot see
- The variogram explained — structure that residuals should not still contain
FAQ
How do I validate an interpolated surface?
Hold out spatial blocks rather than individual points, report RMSE in the variable's units, and publish the coverage statistics — mean and 95th-percentile distance to the nearest sample — alongside.
Is leave-one-out cross-validation reliable?
Only when the samples are well spread. It was 13% pessimistic on an evenly sampled dataset here, and it becomes strongly optimistic on clustered data because held-out points keep their nearby neighbours.
What is block cross-validation?
Splitting the samples into spatial groups and holding out whole groups, so predictions must cross a real gap rather than borrow from a neighbour a few metres away.
Why is my cross-validation score better than my map looks?
Almost certainly clustered sampling. The score is measured at sample locations, which are the best-supported places, and clusters make them better supported still.
Should I use R² or RMSE?
RMSE, in the variable's own units. R² compares against predicting the mean, which is a weak baseline for a spatially structured variable.
How do I know my error where I have no data?
You cannot measure it, but you can bound it by the distance to the nearest sample. Error grew from 18.9 m within 100 m of a sample to 334.3 m beyond a kilometre in the measurement here.
Should I refit the variogram inside each fold?
For a strictly honest score, yes. Fitting it once on all the data leaks information from the held-out fold into the model.