How to Cross-Validate an Interpolated Surface
On this page
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:
random, 500 samples clustered, 377 samples
leave-one-out CV RMSE 55.8 m 37.0 m
true RMSE, every cell 56.0 m 185.3 m
5 spatial blocks CV 188.3 m 196.4 m
On a well-spread sample LOO was right to within half a percent. On a clustered sample of the same terrain it was five times too optimistic. Block cross-validation came close on the clustered sample and was three times too pessimistic on the random one (Example 1). No single scheme 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".
What you get — from The data used on this page
DEM 52 to 1072 m; 500 samples, 12,000 held out
What LOO measures, and block CV for clustered samples
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 — much the problem the average cell of the map faces. That is why LOO works on a well-spread design: 55.8 m against a true 56.0 m here.
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 the clustered design here, the true RMSE was 185.3 m against 56.0 m for the random one — and the clustered design's LOO reported 37.0 m, a better score than the random design'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". On the clustered design that is the right question — five blocks gave 196.4 m against a true 185.3 m. On the random design it is the wrong one: holding out a fifth of the map makes every prediction bridge a hole a kilometre or more wide, and five blocks gave 188.3 m against a true 56.0 m. Match the scheme to how the map will be used.
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, on the clustered design (Example 2):
0 - 100 m 9,680 cells 21.2 m
100 - 250 m 14,020 cells 50.3 m
250 - 500 m 15,753 cells 101.6 m
500 - 1,000 m 25,193 cells 142.2 m
over 1,000 m 46,243 cells 258.7 m
One number would say "about 185 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.
LOO, k-fold and block CV side by side
The data used on this page
The surface is real: the Copernicus DEM GLO-30 over a 10 × 10 km window centred on the summit of Snowdon (Yr Wyddfa), warped to a 30 m grid in EPSG:32630 (x 422,880–432,870, y 5,875,440–5,885,430; 333 × 333 cells, 99.8 km²). The 500 samples are cell centres drawn at random with a fixed seed; 12,000 further cells, none of them samples, are held out as test_points/test_values. grid_points/truth are every cell centre and its elevation, and truth_fn reads the DEM bilinearly at any location, which the error tables score against. idw is the function from Spatial interpolation explained, and surface is its power-2, 12-neighbour surface on that 30 m grid. This page validates both the random sample and a clustered design: 8 centres placed at random, 50 samples round each with a 400 m standard deviation (seed 1), snapped to cells — 377 samples once those outside the window are dropped — as cpoints/cvalues.
import os
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.warp import reproject
from scipy.interpolate import RegularGridInterpolator
bounds = (422_880.0, 5_875_440.0, 432_870.0, 5_885_430.0)
transform = from_origin(bounds[0], bounds[3], 30, 30)
shape = (333, 333)
if not os.path.exists("snowdon_dem_utm30.tif"):
COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
"Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
dem = np.full(shape, np.nan, "float32")
for tile in ("W005", "W004"):
part = np.full(shape, np.nan, "float32")
with rasterio.open(COP.format(tile)) as src:
reproject(rasterio.band(src, 1), part, dst_transform=transform, dst_crs="EPSG:32630",
dst_nodata=np.nan, resampling=Resampling.bilinear)
dem = np.where(np.isnan(dem), part, dem)
with rasterio.open("snowdon_dem_utm30.tif", "w", driver="GTiff", width=shape[1], height=shape[0], count=1,
dtype="float32", crs="EPSG:32630", transform=transform, nodata=np.nan) as dst:
dst.write(dem, 1)
dem = rasterio.open("snowdon_dem_utm30.tif").read(1)
xs = bounds[0] + (np.arange(shape[1]) + 0.5) * 30
ys = bounds[3] - (np.arange(shape[0]) + 0.5) * 30
gx, gy = np.meshgrid(xs, ys)
grid_points = np.column_stack([gx.ravel(), gy.ravel()])
truth = dem.ravel().astype("float64")
rng = np.random.default_rng(0)
order = rng.permutation(truth.size)
sample_idx, test_idx = order[:500], order[500:12_500]
points, values = grid_points[sample_idx], truth[sample_idx]
test_points, test_values = grid_points[test_idx], truth[test_idx]
_bilinear = RegularGridInterpolator((ys[::-1], xs), dem[::-1], bounds_error=False, fill_value=None)
truth_fn = lambda xy: _bilinear(np.column_stack([xy[:, 1], xy[:, 0]]))
from scipy.spatial import cKDTree
def idw(points, values, targets, power=2, k=12):
tree = cKDTree(points)
distance, index = tree.query(targets, k=min(k, len(points)))
distance, index = distance.reshape(len(targets), -1), index.reshape(len(targets), -1)
weights = 1.0 / np.maximum(distance, 1e-9) ** power
predicted = (weights * values[index]).sum(1) / weights.sum(1)
exact = distance[:, 0] < 1e-9
predicted[exact] = values[index[exact, 0]]
return predicted, distance[:, 0]
surface = idw(points, values, grid_points)[0].reshape(shape)
crng = np.random.default_rng(1)
centres = np.column_stack([crng.uniform(bounds[0], bounds[2], 8), crng.uniform(bounds[1], bounds[3], 8)])
raw = (centres[:, None, :] + crng.normal(0, 400, (8, 50, 2))).reshape(-1, 2)
raw = raw[(raw[:, 0] > bounds[0]) & (raw[:, 0] < bounds[2]) & (raw[:, 1] > bounds[1]) & (raw[:, 1] < bounds[3])]
cidx = np.unique(((bounds[3] - raw[:, 1]) // 30).astype(int) * shape[1] + ((raw[:, 0] - bounds[0]) // 30).astype(int))
cpoints, cvalues = grid_points[cidx], truth[cidx]
print(f"DEM {np.nanmin(dem):.0f} to {np.nanmax(dem):.0f} m; {len(points)} samples, {len(test_points):,} held out")DEM 52 to 1072 m; 500 samples, 12,000 held out
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}
predict = lambda p, v, t: idw(p, v, t)[0]
true_error = lambda p, v: float(np.sqrt(np.mean((predict(p, v, grid_points) - truth) ** 2)))
for name, p, v in (("random", points, values), ("clustered", cpoints, cvalues)):
print(f"{name} ({len(p)} samples), true error over every cell {true_error(p, v):.2f}")
_ = compare_validation(p, v, predict)
print("loo_cv on the random sample:", {k: round(x, 3) for k, x in loo_cv(points, values, predict).items()})random (500 samples), true error over every cell 55.99
leave-one-out 55.82
random 5-fold 62.90
5 spatial blocks 188.25
block / LOO ratio 3.37
clustered (377 samples), true error over every cell 185.27
leave-one-out 37.04
random 5-fold 40.82
5 spatial blocks 196.40
block / LOO ratio 5.30
loo_cv on the random sample: {'rmse': 55.824, 'mae': 38.26, 'bias': 2.644, 'r2': 0.92}
The ratio on its own does not diagnose the design: it was 3.4 on the random sample and 5.3 on the clustered one. What differs is which score matches the truth — LOO on the random design, the blocks on the clustered one. Read the coverage statistics in Example 3 to know which situation you are in.
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 rows
_ = error_by_distance(cpoints, cvalues, grid_points, truth, predict) 0- 100 m 9,680 cells RMSE 21.20
100- 250 m 14,020 cells RMSE 50.34
250- 500 m 15,753 cells RMSE 101.63
500- 1000 m 25,193 cells RMSE 142.21
1000- inf m 46,243 cells RMSE 258.68
This 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 report
_ = validation_report(points, values, predict, bounds) n_samples 500
area_km2 99.8
ideal_spacing_m 447
mean_gap_m 228
p95_gap_m 449
max_gap_m 757
block_cv_rmse 188.25
block_cv_bias 34.59
n_blocks 5
caveat block CV RMSE is measured at sample locations; error grows with distance from the nearest sample
Shipping 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.
When LOO is right and when it misleads
Why LOO was right on one design and wrong on the other
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 close to the problem the average cell faces — hence 55.8 m against a true 56.0 m.
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 reported 37.0 m for a surface whose RMSE was 185.3 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 for clustered data
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. On a well-spread sample that is true and large — 188.3 m against a true 56.0 m here — so use LOO there. On clustered data it is the only score that came close, and a surface whose block CV is acceptable is a surface you can defend.
Why R² is a poor choice here
RMSE in the variable's own units is more informative — "±56 m of elevation" is a statement anyone can act on, and "R² = 0.92" (the LOO value here) 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.
Self-prediction, fold variograms and area weighting
- 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
- A species distribution model scores well and predicts badly — why random folds flatter a spatial model
- Interpolated pollution map shows bullseyes — when validation cannot choose the parameter
FAQ
How do I validate an interpolated surface?
Hold out spatial blocks rather than individual points when the samples are clustered (on a well-spread sample, leave-one-out is closer), 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 matched the true error within half a percent on a random sample here, and was five times too optimistic on a clustered one, 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 21.2 m within 100 m of a sample to 258.7 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.