How to Do Spatial Cross-Validation in Python
On this page
Problem statement
Spatial block cross-validation is the standard recommendation for spatial machine learning. Applied without thinking about the deployment scenario, it gives a worse answer than the random k-fold it replaced.
Measured on a real problem — predicting Sentinel-2 NDVI from terrain on Snowdon, 3,000 training cells, repeated over ten random samples (Example 4):
mean R² sd
random 5-fold CV +0.339 0.024
spatial block CV +0.102 0.106
true error over the
whole unsampled map +0.344 0.008
Against the truth, random CV missed by 0.005 and block CV by 0.24, with four times the spread from sample to sample. The scheme to use is the one whose question matches the deployment, and the comparison above is how to check it.
Change the task to genuine extrapolation — train on the western half, predict the east — and both lie:
random CV +0.338
block CV +0.033
actual -0.409
Quick answer
Match the validation to the deployment:
from sklearn.cluster import KMeans
from sklearn.model_selection import GroupKFold, KFold
# predicting inside the sampled area
folds = KFold(n_splits=5, shuffle=True, random_state=0).split(X)
# predicting in unsampled areas within the same region
blocks = KMeans(n_clusters=10, random_state=0, n_init=10).fit_predict(coords)
folds = GroupKFold(n_splits=5).split(X, y, blocks)
# predicting in a genuinely new region: hold one out for real
region = np.where(coords[:, 0] < np.median(coords[:, 0]), "west", "east")
train = region != "east"What you get — from 4. Keep everything inside the fold
[0.429 0.403 0.327 0.312 0.377]
The deployment task, block size, repeats, the fold boundary
1. Decide what the model will be asked to do
- Fill gaps inside a surveyed area. Random k-fold matches this and was accurate to 0.005 R².
- Predict unsampled ground in the same region. Block CV is the right shape, and expect it to be pessimistic.
- Transfer to a new region. Neither works. Hold out a real region.
Writing this down before choosing a splitter is the whole method.
2. Build blocks larger than the autocorrelation range
Blocks smaller than the range still leak: a held-out point has training points within its correlation distance, just in a different fold.
Estimate the range from a variogram of the target, and make blocks at least that size. On the measured data the semivariance of NDVI reached 95% of its variance at 2.5 km (Example 4). Ten KMeans blocks over the 10 km window spanned a median 3.65 km, just clear of it; five blocks are larger still.
3. Use enough blocks, and repeat
Five blocks gave a standard deviation of 0.106 across ten samples of 3,000 cells, ranging from −0.057 to +0.258. A single number from a single sample is false precision. On one sample, five different ten-block layouts varied much less (+0.234 to +0.258, Example 1): the blocks follow the sample, so repeat over samples, not just seeds.
Ten to twenty blocks with a 30% hold-out is a better default, and repeating over several random block layouts and reporting the spread is better still.
4. Keep everything inside the fold
Feature scaling, imputation, feature selection and hyperparameter tuning all leak if fitted on the full dataset before splitting.
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import KFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
sample = np.random.default_rng(0).choice(len(y), 3000, replace=False)
Xs, ys, cs = X[sample], y[sample], coords[sample]
folds = KFold(n_splits=5, shuffle=True, random_state=0)
model = make_pipeline(StandardScaler(), RandomForestRegressor(n_estimators=200, n_jobs=-1, random_state=0))
print(cross_val_score(model, Xs, ys, cv=folds).round(3))[0.429 0.403 0.327 0.312 0.377]
A Pipeline refits every step per fold. Scaling outside it uses the test fold's statistics — an ordinary, non-spatial leak that spatial CV does nothing about.
5. Report the scheme with the score
R² = 0.34 is not a result. R² = 0.34, random 5-fold, for prediction within the sampled extent is.
Buffered blocks and an honest transfer test
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, y and coords hold every remaining cell; the examples sample from them.
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"]
print(f"{item.id}: {shape[1]} x {shape[0]} cells of 20 m, {valid.sum():,} usable "
f"({valid.mean():.1%}); NDVI mean {y.mean():.3f}, sd {y.std():.3f}")S2C_30UVD_20250712_0_L2A: 543 x 509 cells of 20 m, 270,536 usable (97.9%); NDVI mean 0.730, sd 0.145
Example 1 — spatial blocks with the spread reported
import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import r2_score
from sklearn.model_selection import GroupKFold
def spatial_cv(model_factory, X, y, coords, n_blocks=10, n_splits=5,
repeats=5):
"""Block CV repeated over several block layouts."""
scores = []
for seed in range(repeats):
blocks = KMeans(n_clusters=n_blocks, random_state=seed,
n_init=10).fit_predict(coords)
predicted = np.empty(len(y))
for train, test in GroupKFold(n_splits=n_splits).split(X, y, blocks):
predicted[test] = model_factory().fit(X[train], y[train]) \
.predict(X[test])
scores.append(r2_score(y, predicted))
scores = np.array(scores)
print(f" {n_blocks} blocks, {n_splits} folds, {repeats} layouts")
print(f" R2 {scores.mean():+.3f} ± {scores.std():.3f} "
f"(range {scores.min():+.3f} to {scores.max():+.3f})")
if scores.std() > 0.05:
print(" ! unstable — report the spread, not a single number")
return scores
forest = lambda: RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1, random_state=0)
spatial_cv(forest, Xs, ys, cs) 10 blocks, 5 folds, 5 layouts
R2 +0.245 ± 0.008 (range +0.234 to +0.258)
Example 2 — blocks with a buffer, so folds do not touch
import numpy as np
from scipy.spatial import cKDTree
def buffered_blocks(coords, n_blocks=10, buffer_m=500, seed=0):
"""Assign blocks, then drop training points within a buffer of the test fold."""
from sklearn.cluster import KMeans
blocks = KMeans(n_clusters=n_blocks, random_state=seed,
n_init=10).fit_predict(coords)
tree = cKDTree(coords)
for block in np.unique(blocks):
test = blocks == block
if not test.any():
continue
near = set()
for idx in tree.query_ball_point(coords[test], r=buffer_m):
near.update(idx)
train = np.ones(len(coords), bool)
train[test] = False
excluded = np.zeros(len(coords), bool)
excluded[list(near)] = True
train &= ~excluded
print(f" block {block}: {int(test.sum()):5,} test, "
f"{int(train.sum()):5,} train, "
f"{int((~test & excluded).sum()):5,} buffered out")
yield train, test
predicted = np.empty(len(ys))
for train, test in buffered_blocks(cs, n_blocks=10, buffer_m=500):
predicted[test] = forest().fit(Xs[train], ys[train]).predict(Xs[test])
print(f" 10 blocks with a 500 m buffer: R2 {r2_score(ys, predicted):+.3f}") block 0: 349 test, 2,567 train, 84 buffered out
block 1: 333 test, 2,499 train, 168 buffered out
block 2: 299 test, 2,627 train, 74 buffered out
block 3: 244 test, 2,696 train, 60 buffered out
block 4: 349 test, 2,538 train, 113 buffered out
block 5: 258 test, 2,632 train, 110 buffered out
block 6: 257 test, 2,636 train, 107 buffered out
block 7: 267 test, 2,668 train, 65 buffered out
block 8: 343 test, 2,534 train, 123 buffered out
block 9: 301 test, 2,595 train, 104 buffered out
10 blocks with a 500 m buffer: R2 +0.203
The buffer is what makes block CV do what it claims. Without it, a training point one metre outside a block boundary is a near neighbour of a test point one metre inside — exactly the leak blocks were meant to prevent.
The cost is real. Here a 500 m buffer around ten blocks removed 60 to 168 training cells per fold, 2–6% of the training data, and the score fell to +0.203, against +0.234 to +0.258 for unbuffered ten-block layouts of the same sample.
Example 3 — the honest transfer test
import numpy as np
from sklearn.metrics import r2_score, mean_absolute_error
def transfer_test(model_factory, X, y, coords, axis=0):
"""Train on one half, predict the other. Both directions."""
threshold = np.median(coords[:, axis])
left = coords[:, axis] < threshold
results = {}
for label, train, test in (("left -> right", left, ~left),
("right -> left", ~left, left)):
model = model_factory().fit(X[train], y[train])
predicted = model.predict(X[test])
results[label] = {
"r2": r2_score(y[test], predicted),
"mae": mean_absolute_error(y[test], predicted),
"n_train": int(train.sum()), "n_test": int(test.sum()),
}
print(f" {label}: R2 {results[label]['r2']:+.3f}, "
f"MAE {results[label]['mae']:.4f} "
f"({int(train.sum()):,} train, {int(test.sum()):,} test)")
worst = min(r["r2"] for r in results.values())
if worst < 0:
print(" ! the model does not transfer — worse than predicting the "
"training mean in at least one direction")
return results
six = np.random.default_rng(1).choice(len(y), 6000, replace=False) # 3,000 cells either side
transfer_test(forest, X[six], y[six], coords[six]) left -> right: R2 +0.083, MAE 0.1025 (2,997 train, 3,003 test)
right -> left: R2 +0.086, MAE 0.0838 (3,003 train, 2,997 test)
On this sample of 6,000 cells both directions transferred, barely: +0.083 and +0.086. Trained on a different 3,000 western cells and applied to every eastern cell (Example 4), the same model scored −0.409. A transfer score is itself unstable, which is one more reason to run it on several samples. A negative R² means the model is worse than predicting the training mean. It is common in genuine transfer tests and almost never reported, because most pipelines never run one.
Example 4 — the tables at the top of this page, and the range
from sklearn.model_selection import GroupKFold
def scheme_scores(Xs, ys, cs, seed):
model = lambda: RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1, random_state=seed)
def cv(splitter):
p = np.empty(len(ys))
for train, test in splitter:
p[test] = model().fit(Xs[train], ys[train]).predict(Xs[test])
return r2_score(ys, p)
blocks = KMeans(n_clusters=5, random_state=seed, n_init=10).fit_predict(cs)
return (cv(KFold(n_splits=5, shuffle=True, random_state=seed).split(Xs)),
cv(GroupKFold(n_splits=5).split(Xs, ys, blocks)), model)
runs = []
for seed in range(10):
pick = np.random.default_rng(seed).choice(len(y), 3000, replace=False)
others = np.setdiff1d(np.arange(len(y)), pick)
random_r2, block_r2, model = scheme_scores(X[pick], y[pick], coords[pick], seed)
runs.append((random_r2, block_r2, r2_score(y[others], model().fit(X[pick], y[pick]).predict(X[others]))))
for name, column in zip(("random 5-fold CV", "5-block spatial CV", "true map R²"), np.array(runs).T):
print(f"{name:20} mean {column.mean():+.3f} sd {column.std():.3f} range {column.min():+.3f} to {column.max():+.3f}")
west = coords[:, 0] < np.median(coords[:, 0])
pick = np.random.default_rng(0).choice(np.nonzero(west)[0], 3000, replace=False)
random_r2, block_r2, model = scheme_scores(X[pick], y[pick], coords[pick], 0)
print(f"trained in the west: random CV {random_r2:+.3f}, block CV {block_r2:+.3f}, actual in the east "
f"{r2_score(y[~west], model().fit(X[pick], y[pick]).predict(X[~west])):+.3f}")
# Empirical variogram of the target, 5,000 cells, 250 m lags
from scipy.spatial.distance import pdist
v = np.random.default_rng(2).choice(len(y), 5000, replace=False)
d, g = pdist(coords[v]), 0.5 * pdist(y[v, None], "sqeuclidean")
edges = np.arange(0, 6001, 250)
gamma = np.array([g[(d >= a) & (d < b)].mean() for a, b in zip(edges[:-1], edges[1:])])
sill = y.var()
reach = edges[1:][np.argmax(gamma >= 0.95 * sill)] if (gamma >= 0.95 * sill).any() else None
print(f"variance {sill:.4f}; semivariance at 250 m {gamma[0]:.4f}, 1 km {gamma[3]:.4f}, 3 km {gamma[11]:.4f}, "
f"6 km {gamma[-1]:.4f}; first reaches 95% of the variance at {reach} m")
blocks10 = KMeans(n_clusters=10, random_state=0, n_init=10).fit_predict(cs)
print(f"10 KMeans blocks on the sample: median block span "
f"{np.median([np.ptp(cs[blocks10 == k], axis=0).max() for k in range(10)]):,.0f} m")random 5-fold CV mean +0.339 sd 0.024 range +0.303 to +0.375
5-block spatial CV mean +0.102 sd 0.106 range -0.057 to +0.258
true map R² mean +0.344 sd 0.008 range +0.335 to +0.360
trained in the west: random CV +0.338, block CV +0.033, actual in the east -0.409
variance 0.0210; semivariance at 250 m 0.0092, 1 km 0.0160, 3 km 0.0213, 6 km 0.0227; first reaches 95% of the variance at 2500 m
10 KMeans blocks on the sample: median block span 3,650 m
Why a block has to outrun the autocorrelation
Why block CV was pessimistic here, and why no scheme can score a new region, is the subject of spatial leakage explained. What matters when building the folds is their size.
Why blocks need to exceed the autocorrelation range
The point of a block is that test points have no training points within their correlation distance.
If blocks are smaller than the range, a test point near a block edge has training points from the adjacent block within range, and the fold leaks. The score then sits between random CV and true block CV, and nothing indicates which.
Estimating the range from a variogram of the target — not of the residuals, not of a covariate — is the defensible way to size them.
Why to repeat over layouts
The score depends on which regions each fold happens to hold out, and the blocks are cut from wherever the sample fell. Re-seeding KMeans on one sample moves them little; drawing a new sample moves them a lot.
Measured: ten samples of 3,000 cells gave block-CV scores from −0.057 to +0.258, a range of 0.32 R² from the same landscape and an identical model.
Reporting one of those numbers to three decimal places is false precision by a wide margin. Repeat, and report the mean and spread.
Pipelines, tuning folds and grouping by cluster
- Match the scheme to the deployment, not to a rule of thumb.
- Blocks must exceed the autocorrelation range, or folds still leak.
- Buffer between folds to remove edge leakage, and accept the extra pessimism.
- Repeat over samples and block layouts. Block CV had a standard deviation of 0.106 across samples here.
- Use a
Pipelineso scaling and imputation are refit per fold. - Do not tune on the folds you report.
- A negative R² is meaningful — worse than the mean.
- Group by cluster, not just by location, when samples come in plots or transects.
Internal links
- Spatial leakage explained: why random cross-validation lies — the measurements behind this
- Spatial machine learning explained — why location breaks the usual rules
- Extrapolation in space explained: the area of applicability — what to do when transfer cannot be validated
- How to evaluate a spatial model honestly — what else to report
- My model scores 0.95 in testing and fails in the field — the symptom
- How to cross-validate an interpolated surface — the same problem in interpolation
- Spatial autocorrelation explained — measuring the range
- The variogram explained — sizing the blocks
- How to build a land-use regression model for NO2 — leave-one-out with 17 stations
FAQ
How do I do spatial cross-validation in Python?
Cluster the coordinates into blocks with KMeans, use GroupKFold with those blocks, and repeat over several block layouts reporting the spread.
Is spatial CV always better than random CV?
No. For predicting inside a sampled area, random CV matched the truth to 0.005 R² while block CV was pessimistic by 0.24.
How big should the blocks be?
Larger than the target's autocorrelation range, estimated from a variogram. Smaller blocks still leak across fold boundaries.
How many blocks should I use?
Enough that each block is larger than the autocorrelation range — about 2.5 km here — with several repeats. Five blocks gave a standard deviation of 0.106 across samples here.
Should I buffer between folds?
Yes if you want block CV to mean what it claims. It costs training data and makes the estimate more pessimistic, which is the point.
Can spatial CV tell me if my model works in a new region?
No. It promised +0.033 where a real west-to-east transfer gave −0.409. Hold out a real region instead.
Where do scaling and imputation go?
Inside a Pipeline, so they are refit per fold. Fitting them on the full dataset leaks in the ordinary, non-spatial way.