How to Evaluate a Spatial Model Honestly
On this page
Problem statement
A single score is not an evaluation of a map. It is one number summarising a model's performance on the data it was validated against, and every map has parts where it is much better and much worse than that number.
Measured on a real problem — Snowdon NDVI predicted from terrain, ten samples of 3,000 cells — the same model reported four ways:
random 5-fold CV, inside the sampled area +0.339 sd 0.024
spatial block CV, 5 blocks +0.102 sd 0.106
true error over the whole unsampled map +0.344 sd 0.008
the model applied to the adjacent half -0.409
All four are correct answers to different questions. Reporting one without saying which question it answers is the most common failure in applied spatial modelling.
Quick answer
An honest evaluation has five parts:
report = {
"scheme": "spatial block CV, 10 blocks, 5 folds, 5 layouts",
"metric": {"r2": 0.245, "r2_sd": 0.008, "mae": 0.0866},
"deployment": "prediction within the sampled extent",
"applicability": "93.6% of the map is inside the area of applicability",
"residuals": {"spatial_autocorrelation": 0.554, "bias": 0.007},
}These are the values Examples 1 to 3 produce for a random forest on 3,000 Snowdon cells. The scheme and the deployment are what make the metric interpretable. Without them, R² = 0.245 is not information.
What you get — from The data used on this page
S2C_30UVD_20250712_0_L2A: 543 x 509 cells of 20 m, 270,536 usable (97.9%); training on 3,000
Scheme, spread, units, residual maps, applicability
1. State the validation scheme and the deployment task
These are the two facts that let a reader interpret everything else. R² = 0.374, random 5-fold, for gap-filling within the sampled extent is a result; R² = 0.374 is not.
2. Report the spread, not a point estimate
Five-block CV over ten different samples of the same landscape gave a standard deviation of 0.106, ranging from −0.057 to +0.258. Quoting one of those to three decimal places is false precision by a wide margin. Re-seeding the blocks on one sample is not enough to see it: five ten-block layouts of one sample varied only from +0.234 to +0.258 (Example 1).
Repeat over samples, over layouts and over seeds, and report the mean and spread.
3. Use metrics in the variable's units
MAE and RMSE in the variable's own units — "±0.087 NDVI" here — are directly interpretable and comparable across studies.
Report bias separately: a model can have a good RMSE and a systematic offset.
4. Map the residuals
A single error number hides where the errors are. Two checks:
- Spatial autocorrelation of residuals. Structured residuals mean an unmodelled spatial process. Moran's I on the residuals is the standard test.
- Error against covariate values. Systematically larger errors at high elevation or steep slope means the model fits one part of the range and not another.
5. Report the area of applicability
The share of the map where the covariates resemble the training data. Outside it, the model is extrapolating and the metric does not apply at all.
A map that is 60% inside its area of applicability with an RMSE of 0.1 is a very different product from one that is 98% inside with the same RMSE.
A full evaluation report with residual diagnostics
The data used on this page
A 10 km window on Snowdon, 4.16–4.00° W, 53.03–53.12° N. The target is NDVI from the least cloudy Sentinel-2 L2A scene of June–September 2025 on Earth Search, on the scene's 20 m grid; the features are terrain from the Copernicus DEM GLO-30 resampled to the same grid — elevation, slope, northness, eastness and a topographic position index over 300 m. Cells the scene classifies as water, cloud, cloud shadow, snow or no data are dropped. X_grid and y_grid hold every remaining cell — the map to predict — and X, y and coords a random 3,000 of them, the training sample.
import numpy as np
import rasterio
import rasterio.merge
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.warp import reproject, transform_bounds
from scipy import ndimage
AOI = (-4.16, 53.03, -4.00, 53.12)
items = Client.open("https://earth-search.aws.element84.com/v1").search(
collections=["sentinel-2-l2a"], bbox=list(AOI), datetime="2025-06-01/2025-10-01",
query={"eo:cloud_cover": {"lt": 60}}).items()
item = min(items, key=lambda i: i.properties["eo:cloud_cover"])
with rasterio.open(item.assets["scl"].href) as ds: # the 20 m grid
win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
transform=ds.transform).round_offsets().round_lengths()
scl, transform, crs = ds.read(1, window=win), ds.window_transform(win), ds.crs
bounds = rasterio.windows.bounds(win, ds.transform)
shape = scl.shape
band = {}
for name in ("red", "nir"): # 10 m bands averaged onto it
with rasterio.open(item.assets[name].href) as ds:
band[name] = ds.read(1, window=rasterio.windows.from_bounds(*bounds, transform=ds.transform),
out_shape=shape, resampling=Resampling.average).astype("float32") / 1e4
ndvi = (band["nir"] - band["red"]) / (band["nir"] + band["red"])
tiles = [rasterio.open(f"https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{t}_00_DEM/"
f"Copernicus_DSM_COG_10_N53_00_{t}_00_DEM.tif") for t in ("W005", "W004")]
mosaic, mosaic_transform = rasterio.merge.merge(tiles, bounds=(-4.18, 53.02, -3.98, 53.13))
elev = np.full(shape, np.nan, dtype="float32")
reproject(mosaic[0], elev, src_transform=mosaic_transform, src_crs=tiles[0].crs,
dst_transform=transform, dst_crs=crs, resampling=Resampling.bilinear)
dz_dy, dz_dx = np.gradient(elev, -20, 20) # rows run south
slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
aspect = np.arctan2(-dz_dx, -dz_dy) # direction the slope faces, from north
tpi = elev - ndimage.uniform_filter(elev, size=15) # against the 300 m neighbourhood
valid = np.isin(scl, (2, 4, 5, 7)) & np.isfinite(ndvi) & np.isfinite(elev)
rows, cols = np.nonzero(valid)
xs, ys = rasterio.transform.xy(transform, rows, cols)
coords = np.column_stack([xs, ys])
X = np.column_stack([elev[valid], slope[valid], np.cos(aspect[valid]), np.sin(aspect[valid]), tpi[valid]])
y = ndvi[valid]
FEATURES = ["elevation", "slope", "northness", "eastness", "tpi"]
X_grid, y_grid, coords_grid = X, y, coords
sample = np.random.default_rng(0).choice(len(y_grid), 3000, replace=False)
X, y, coords = X_grid[sample], y_grid[sample], coords_grid[sample]
print(f"{item.id}: {shape[1]} x {shape[0]} cells of 20 m, {valid.sum():,} usable "
f"({valid.mean():.1%}); training on {len(y):,}")S2C_30UVD_20250712_0_L2A: 543 x 509 cells of 20 m, 270,536 usable (97.9%); training on 3,000
Example 1 — a full evaluation report
import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import mean_absolute_error, r2_score
from sklearn.model_selection import GroupKFold, KFold
def evaluate(model_factory, X, y, coords, feature_names,
deployment="prediction within the sampled extent",
n_blocks=10, repeats=5):
"""Scheme, metrics with spread, residual diagnostics, applicability."""
def cv(splitter):
predicted = np.empty(len(y))
for train, test in splitter:
predicted[test] = model_factory().fit(X[train], y[train]) \
.predict(X[test])
return predicted
random_pred = cv(KFold(5, shuffle=True, random_state=0).split(X))
block_scores, block_pred = [], None
for seed in range(repeats):
blocks = KMeans(n_blocks, random_state=seed,
n_init=10).fit_predict(coords)
predicted = cv(GroupKFold(5).split(X, y, blocks))
block_scores.append(r2_score(y, predicted))
if block_pred is None:
block_pred = predicted
block_scores = np.array(block_scores)
residual = y - block_pred
print(f" deployment task : {deployment}")
print(f" random 5-fold : R2 {r2_score(y, random_pred):+.3f}, "
f"MAE {mean_absolute_error(y, random_pred):.4f}")
print(f" block CV : R2 {block_scores.mean():+.3f} "
f"± {block_scores.std():.3f} over {repeats} layouts "
f"(range {block_scores.min():+.3f} to {block_scores.max():+.3f})")
print(f" MAE (block) : {mean_absolute_error(y, block_pred):.4f} "
f"in units of the target")
print(f" bias : {residual.mean():+.4f}")
print(f" target sd : {y.std():.4f} (the baseline RMSE)")
if block_scores.std() > 0.05:
print(" ! block CV is unstable — report the range, not a point value")
return {"random_r2": float(r2_score(y, random_pred)),
"block_r2_mean": float(block_scores.mean()),
"block_r2_sd": float(block_scores.std()),
"mae": float(mean_absolute_error(y, block_pred)),
"bias": float(residual.mean()), "deployment": deployment,
"residual": residual}
from sklearn.ensemble import RandomForestRegressor
forest = lambda: RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1, random_state=0)
metrics = evaluate(forest, X, y, coords, FEATURES)
print(f" MAE of predicting the mean everywhere: {mean_absolute_error(y, np.full_like(y, y.mean())):.4f}") deployment task : prediction within the sampled extent
random 5-fold : R2 +0.374, MAE 0.0778
block CV : R2 +0.245 ± 0.008 over 5 layouts (range +0.234 to +0.258)
MAE (block) : 0.0866 in units of the target
bias : +0.0070
target sd : 0.1459 (the baseline RMSE)
MAE of predicting the mean everywhere: 0.1047
Printing the target's standard deviation next to the MAE gives the reader the baseline. Here the block-CV MAE is 0.087 NDVI, against 0.105 for predicting the mean everywhere (the target sd is 0.146) — a modest improvement, which is what the R² of +0.245 says less usefully.
Example 2 — residual diagnostics
import numpy as np
from scipy.spatial import cKDTree
def residual_diagnostics(residual, coords, X=None, feature_names=None,
k=8):
"""Is the error spatially structured, or related to a covariate?"""
tree = cKDTree(coords)
_, neighbours = tree.query(coords, k=k + 1)
lag = residual[neighbours[:, 1:]].mean(axis=1)
centred = residual - residual.mean()
lag_centred = lag - lag.mean()
morans = float((centred * lag_centred).sum() /
np.sqrt((centred ** 2).sum() * (lag_centred ** 2).sum()))
print(f" residual mean {residual.mean():+.4f}, sd {residual.std():.4f}")
print(f" spatial autocorrelation of residuals (k={k}): {morans:+.3f}")
if abs(morans) > 0.2:
print(" ! residuals are spatially structured — an unmodelled "
"spatial process remains")
if X is not None and feature_names is not None:
print(" error against each covariate (correlation with |residual|):")
for i, name in enumerate(feature_names):
correlation = float(np.corrcoef(X[:, i], np.abs(residual))[0, 1])
flag = " <- heteroscedastic" if abs(correlation) > 0.3 else ""
print(f" {name:16} {correlation:+.3f}{flag}")
return {"morans_i": morans, "bias": float(residual.mean())}
diagnostics = residual_diagnostics(metrics["residual"], coords, X, FEATURES) residual mean +0.0070, sd 0.1255
spatial autocorrelation of residuals (k=8): +0.554
! residuals are spatially structured — an unmodelled spatial process remains
error against each covariate (correlation with |residual|):
elevation +0.108
slope +0.101
northness -0.082
eastness +0.014
tpi +0.007
The correlation between a covariate and the absolute residual is the heteroscedasticity check. A model that is accurate on gentle slopes and poor on steep ones has a real limitation that RMSE averages away. Here no covariate passed ±0.11, but the residuals' spatial autocorrelation was +0.554: nearby errors are much alike, so something that varies smoothly across the window — soil, land use, exposure — is missing from the terrain covariates.
Example 3 — the report that ships with the map
import json
def evaluation_report(metrics, diagnostics, applicability_share,
scheme, deployment, features, n_train, path=None):
"""Everything a reader needs, in one serialisable object."""
report = {
"deployment_task": deployment,
"validation_scheme": scheme,
"n_training_samples": int(n_train),
"features": list(features),
"metrics": {
"r2_mean": round(metrics["block_r2_mean"], 3),
"r2_sd": round(metrics["block_r2_sd"], 3),
"mae": round(metrics["mae"], 4),
"bias": round(metrics["bias"], 4),
},
"residuals": {
"spatial_autocorrelation": round(diagnostics["morans_i"], 3),
},
"area_of_applicability": {
"share_inside": round(applicability_share, 3),
"note": "outside this share the model is extrapolating and the "
"metrics above do not apply",
},
"caveats": [
"metrics are estimated at the training locations",
"error grows with distance from the training data",
f"the model cannot predict outside the training target range",
],
}
text = json.dumps(report, indent=2)
print(text)
if path:
with open(path, "w") as handle:
handle.write(text)
return report
# the applicability share: map cells whose scaled covariates are as close to the training
# data as 95% of training cells are to each other (k = 5 nearest neighbours)
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler().fit(X)
tree = cKDTree(scaler.transform(X))
threshold = np.quantile(tree.query(scaler.transform(X), k=6)[0][:, 1:].mean(axis=1), 0.95)
share = float((tree.query(scaler.transform(X_grid), k=5)[0].mean(axis=1) <= threshold).mean())
evaluation_report(metrics, diagnostics, share, "spatial block CV, 10 blocks, 5 folds, 5 layouts",
metrics["deployment"], FEATURES, len(y), path="evaluation.json"){
"deployment_task": "prediction within the sampled extent",
"validation_scheme": "spatial block CV, 10 blocks, 5 folds, 5 layouts",
"n_training_samples": 3000,
"features": [
"elevation",
"slope",
"northness",
"eastness",
"tpi"
],
"metrics": {
"r2_mean": 0.245,
"r2_sd": 0.008,
"mae": 0.0866,
"bias": 0.007
},
"residuals": {
"spatial_autocorrelation": 0.554
},
"area_of_applicability": {
"share_inside": 0.936,
"note": "outside this share the model is extrapolating and the metrics above do not apply"
},
"caveats": [
"metrics are estimated at the training locations",
"error grows with distance from the training data",
"the model cannot predict outside the training target range"
]
}
The caveats list is the part that survives summarisation. When someone reduces the model to a single number in a slide, the caveats are what they should have carried with it — and writing them down at least makes their omission a choice.
Why R² is weak for spatial targets
Why R² is a weak metric here
For a spatially structured variable, predicting the mean is a very weak baseline, so
Report MAE or RMSE in the target's units, with the target's standard deviation for context.
Why the spread matters as much as the mean
Five-block CV's standard deviation was 0.106 across ten samples of the same landscape with an identical model — a range of 0.32 R².
Two papers reporting −0.06 and +0.26 for the same method in the same place are not disagreeing about anything except which cells they happened to sample, and so which blocks the clustering produced.
Repeating and reporting the spread is what makes the number reproducible. It also, usefully, tells the reader how much to trust it.
Why residual autocorrelation matters
If the residuals are spatially structured — nearby errors similar — the model has missed a spatial process. Something varies smoothly across the area that the covariates do not capture.
That is actionable in a way an error number is not: it points at a missing covariate, or at a case for adding a spatial term such as kriging the residuals.
Residuals that look like noise mean the covariates have captured the structure, which is the goal.
Why applicability belongs in the metric report
The metric was estimated at the training locations. Wherever the map's covariates fall outside the training range, the metric does not describe the error — it describes the error somewhere else.
A model with an RMSE of 0.1 covering 98% of its map inside the area of applicability is a usable product. The same RMSE covering 55% is a research result with a large asterisk.
Reporting both turns a number into a statement about a map.
Bias, residual autocorrelation and written caveats
- State the scheme and the deployment task with every metric.
- Report the spread over repeats, not a point estimate.
- Use the target's units; give its standard deviation as the baseline.
- Report bias separately from RMSE.
- Test residual autocorrelation — structure means a missing covariate.
- Check error against covariate values for heteroscedasticity.
- Report the applicability share.
- Write the caveats down, because they are the first thing dropped in a summary.
Internal links
- Spatial leakage explained: why random cross-validation lies — why the scheme must be stated
- How to do spatial cross-validation in Python — producing the numbers
- Extrapolation in space explained: the area of applicability — the applicability share
- My model scores 0.95 in testing and fails in the field — what a bare metric hides
- Spatial autocorrelation explained — the residual test
- How to calculate Moran's I in Python — computing it properly
- My prediction raster is striped, blocky or full of NoData — what the map shows that metrics do not
- How to cross-validate an interpolated surface — the same argument for interpolation
FAQ
What should I report for a spatial model?
The validation scheme, the deployment task, a metric in the target's units with its spread, residual diagnostics, and the share of the map inside the area of applicability.
Should I use R² or RMSE?
RMSE or MAE, in the target's units, with the target's standard deviation as context. R² compares against predicting the mean, which is a weak baseline for spatial data.
Why report the spread of the score?
Because five-block CV over ten samples of one landscape gave a standard deviation of 0.106. A single number to three decimals is false precision.
What does spatially autocorrelated residual mean?
That the model has missed a spatial process. Something varies smoothly across the area that the covariates do not capture — which points at a missing covariate.
How do I check for heteroscedasticity?
Correlate each covariate with the absolute residual. A strong correlation means the model is accurate over part of the range and not the rest.
Why report the area of applicability with the metric?
Because the metric was estimated at the training locations. Outside the applicability area it describes the error somewhere else entirely.
What if my model has a good RMSE and a bias?
Report both. A systematic offset is a different problem from scatter, and RMSE conflates them.