Choosing a Spatial Cross-Validation Strategy

TL;DR: Choose the split from the deployment task, not from the data. In-fill at existing density takes KFold(shuffle=True); whole-surface mapping takes GroupKFold over blocks whose side is about twice the deployment gap and at least the residual variogram range; transfer to a new region takes strata; few clustered samples take a buffered leave-one-out. On one dataset those four report 0.78, 0.51, 0.46 and 0.33.

Why This Matters #

Every entry in the Cross-Validation Strategies family produces a defensible number from the same data and the same model, and those numbers can differ by more than a factor of two. That spread is not a sign that one design is right and the others are wrong. Each simulates a different distance between the point being predicted and the nearest point that trained the model, and each is honest about the deployment it corresponds to. The mistake is not picking the wrong splitter; it is picking one without ever stating what the model will be asked to do.

So the whole family reduces to a single question asked before any code runs: how far from its training data will this model have to predict? Call that the deployment gap. The validation design must reproduce it. A model validated on a harder split than it faces is under-sold — you will discard something that would have worked. A model validated on an easier split is dangerous, because the number you published will not survive contact with the map. Everything below is machinery for measuring the deployment gap, measuring the gap each candidate split creates, and making the two agree.

Environment and Version Pinning #

Nothing here needs a dedicated cross-validation package: scikit-learn supplies the splitters, scipy the neighbour queries, and scikit-gstat the residual variogram that sizes the blocks.

pip install "geopandas==1.0.1" "scikit-learn==1.5.2" "scikit-gstat==1.0.18" \
            "numpy==1.26.4" "pandas==2.2.3" "scipy==1.14.1" "shapely==2.0.6"

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from scipy.spatial import cKDTree
from skgstat import Variogram
from sklearn.cluster import KMeans
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import KFold, GroupKFold, cross_val_predict

Step-by-Step Implementation #

1. Build a survey with realistic clustering #

Real soil and ecological surveys are not uniform random samples. They are clumps around access points, and that clumping is what makes the choice of split consequential.

rng = np.random.default_rng(2026)
EXTENT = 40_000.0                     # 40 km square, metric CRS
N_CENTRES, PER_CENTRE = 40, 20        # 800 samples in 40 loose clumps

centres = rng.uniform(2_000, EXTENT - 2_000, size=(N_CENTRES, 2))
pts = np.repeat(centres, PER_CENTRE, axis=0) + rng.normal(
    0, 900, size=(N_CENTRES * PER_CENTRE, 2))
pts = np.clip(pts, 0.0, EXTENT)

elev = 300 + 0.004 * pts[:, 0] + 40 * np.sin(pts[:, 1] / 6_000)
ndvi = 0.35 + 0.25 * np.cos(pts[:, 0] / 9_000) + 0.10 * rng.normal(size=len(pts))
field = (1.1 * np.sin(pts[:, 0] / 5_200) * np.cos(pts[:, 1] / 4_800)
         + 0.6 * np.sin((pts[:, 0] + pts[:, 1]) / 7_500))
soc = (2.4 + 0.0025 * (elev - 300) + 1.6 * ndvi + field
       + rng.normal(0, 0.25, len(pts)))

gdf = gpd.GeoDataFrame(
    {"soc": soc, "elev": elev, "ndvi": ndvi, "x": pts[:, 0], "y": pts[:, 1]},
    geometry=[Point(a, b) for a, b in pts], crs="EPSG:32633")

X = gdf[["elev", "ndvi", "x", "y"]].to_numpy()
y = gdf["soc"].to_numpy()
XY = gdf[["x", "y"]].to_numpy()
print(f"{len(gdf)} samples, sd(soc) = {y.std(ddof=1):.3f}")

800 samples, sd(soc) = 0.900

2. Measure the deployment gap #

If the model is going to be rasterised across the study area, the deployment gap is the distance from a prediction cell to the nearest sample. Compute it directly — it takes four lines and it settles the argument.

step = 100.0
gx, gy = np.meshgrid(np.arange(step / 2, EXTENT, step),
                     np.arange(step / 2, EXTENT, step))
grid = np.column_stack([gx.ravel(), gy.ravel()])
d_deploy, _ = cKDTree(XY).query(grid, k=1)
print(f"deployment gap: median {np.median(d_deploy) / 1000:.2f} km, "
      f"90th pct {np.quantile(d_deploy, 0.90) / 1000:.2f} km")

deployment gap: median 2.58 km, 90th pct 5.31 km

Half the map lies more than 2.58 km from any sample, and a tenth lies more than 5.31 km away. No split that leaves training points a few hundred metres from every test point is simulating this task.

3. Read the residual variogram range #

The block size has a second constraint: folds should not exchange correlated information. That distance is the range of the variogram of the model’s residuals, not of the response. Fitting and interpreting it is covered in Estimating Nugget, Sill & Range Parameters.

base = RandomForestRegressor(n_estimators=400, min_samples_leaf=3,
                             random_state=0, n_jobs=-1)
resid = y - cross_val_predict(base, X, y,
                              cv=KFold(10, shuffle=True, random_state=0))
V = Variogram(XY, resid, model="spherical", n_lags=18, maxlag=15_000)
print(V)
range_m = float(V.parameters[0])

spherical Variogram
-------------------
Estimator:         matheron
Effective Range:   4812.44
Sill:              0.4131
Nugget:            0.0725

4. Define the four candidate splits #

Each design is a few lines. GroupKFold does the work for both blocking and stratification; only the grouping variable changes.

def block_groups(xy, block_m):
    """Square blocks of side block_m; returns one integer group per sample."""
    ij = np.floor(xy / block_m).astype(int)
    _, groups = np.unique(ij, axis=0, return_inverse=True)
    return groups


def buffered_loo(xy, buffer_m):
    """Hold out one sample and every training sample within buffer_m of it."""
    tree = cKDTree(xy)
    all_idx = np.arange(len(xy))
    for i in all_idx:
        excluded = np.asarray(tree.query_ball_point(xy[i], r=buffer_m))
        train = np.setdiff1d(all_idx, excluded, assume_unique=True)
        yield train, np.array([i])


folds_random = list(KFold(n_splits=10, shuffle=True, random_state=0).split(X))
folds_block = list(GroupKFold(n_splits=10).split(
    X, y, groups=block_groups(XY, 5_000.0)))
folds_bloo = list(buffered_loo(XY, 5_000.0))

strata = KMeans(n_clusters=5, n_init=10, random_state=0).fit_predict(
    gdf[["elev", "ndvi"]].to_numpy())
folds_env = list(GroupKFold(n_splits=5).split(X, y, groups=strata))

The stratification groups on the covariates rather than on coordinates, which is what makes it a proxy for transfer: each held-out fold is a slice of covariate space the model never saw, wherever it happens to sit on the map.

5. Score every split through one harness #

The only way the comparison means anything is if the model, features and hyper-parameters are byte-identical across designs. Record the prediction distance at the same time as the accuracy.

def evaluate(folds, name):
    pred = np.full(len(y), np.nan)
    gaps = []
    for tr, te in folds:
        model = RandomForestRegressor(n_estimators=400, min_samples_leaf=3,
                                      random_state=0, n_jobs=-1)
        model.fit(X[tr], y[tr])
        pred[te] = model.predict(X[te])
        d, _ = cKDTree(XY[tr]).query(XY[te], k=1)
        gaps.append(d)
    gaps = np.concatenate(gaps)
    return {"strategy": name,
            "R2": r2_score(y, pred),
            "RMSE": float(np.sqrt(mean_squared_error(y, pred))),
            "median_gap_km": float(np.median(gaps)) / 1000.0}


results = pd.DataFrame([
    evaluate(folds_random, "random 10-fold"),
    evaluate(folds_block,  "spatial block, 5 km"),
    evaluate(folds_bloo,   "buffered LOO, 5 km"),
    evaluate(folds_env,    "environmental strata (k=5)"),
])
print(results.to_string(index=False, float_format=lambda v: f"{v:.3f}"))

                  strategy    R2  RMSE  median_gap_km
            random 10-fold 0.780 0.422          0.121
       spatial block, 5 km 0.510 0.630          2.412
        buffered LOO, 5 km 0.460 0.661          5.284
environmental strata (k=5) 0.330 0.736         11.043

6. Sweep the block size to see what it controls #

Blocking is the only design with a free dial, and the dial has a simple meaning.

for b in (1_000, 2_000, 5_000, 10_000):
    folds = list(GroupKFold(n_splits=10).split(
        X, y, groups=block_groups(XY, float(b))))
    r = evaluate(folds, f"{b / 1000:g} km blocks")
    print(f"{r['strategy']:>12}: R2 = {r['R2']:.3f}  "
          f"median gap = {r['median_gap_km']:.2f} km")

 1 km blocks: R2 = 0.731  median gap = 0.61 km
 2 km blocks: R2 = 0.680  median gap = 1.02 km
 5 km blocks: R2 = 0.510  median gap = 2.41 km
10 km blocks: R2 = 0.402  median gap = 4.83 km

The median gap is almost exactly half the block side in every row. That is the rule worth memorising: set the block side to twice the deployment gap. Here the deployment gap of 2.58 km asks for about 5.2 km, and the residual range of 4.81 km says a 5 km block is also large enough for the folds to be near-independent. Both constraints land on the same number, which is the comfortable case; the detailed mechanics are in Spatial Block Cross-Validation in Python.

Interpreting the Output #

Read each row as a pair: an accuracy and the distance that accuracy was earned over. Random 10-fold reports R2 = 0.780 at a median gap of 121 metres, which is a fair account of predicting a missing sample inside an existing clump and a wild overstatement of anything else. Spatial blocking at 5 km reports 0.510 at 2.41 km, sitting inside the 2–3 km band the deployment gap occupies. The buffered leave-one-out design reports 0.460 at 5.28 km, which is close to the 90th percentile of the deployment gap rather than its median — it is answering about the sparse tail of the map. Environmental stratification reports 0.330 at 11.04 km, a proper transfer test.

What good looks like is a chosen strategy whose median gap lands within a factor of about 1.5 of the measured deployment gap, and an RMSE you would be willing to attach to the map legend. The warning signs are two. First, a median gap below a few hundred metres on clustered data, which almost always means duplicate or near-duplicate locations are straddling the fold boundary. Second, a gap far beyond the deployment gap, which produces a pessimistic figure that will be used against the model in review — the 10 km block row above is already in that territory.

Critical Best Practices #

Size blocks from the residual variogram, not the raw one #

The variogram of soc itself carries the covariate trend and will show a longer range than the residual variogram, because elevation and NDVI vary smoothly across the whole square. What matters for fold independence is the correlation the model cannot already explain, so fit the variogram to out-of-fold residuals. On this dataset the residual effective range is 4.81 km; the raw response range is considerably longer, and using it would have pushed you to blocks large enough to report 0.40 for a task that deserves 0.51.

Never choose the strategy after seeing the scores #

Fixing the validation design after inspecting four candidate numbers makes it a tuned parameter, and the thing you tuned it against was the estimate of generalisation itself. Write the deployment description and the resulting choice into the project notes before running evaluate, and treat the other three rows as diagnostics of the split rather than candidates for the headline. Reporting all four with their median gaps is fine and honest; reporting the most flattering one as the accuracy is not.

Keep everything but the fold assignment identical #

random_state, n_estimators, min_samples_leaf and the feature matrix must be the same across designs, otherwise the spread mixes validation effects with model effects. The harness above constructs a fresh RandomForestRegressor with fixed arguments inside every fold for exactly this reason. If you tune hyper-parameters, the tuning must happen in an inner loop that uses the same spatial split as the outer one, or the leakage you removed outside comes back through the inner folds.

Report the prediction distance beside the accuracy #

median_gap_km costs one cKDTree query per fold and it is the only column that makes the accuracy interpretable to someone else. A reader who sees R2 = 0.51, median gap 2.41 km can decide whether that matches their own use; a reader who sees R2 = 0.51 alone cannot. It also catches silent failures — an unexpectedly small gap under a blocking design usually means the grouping variable was passed positionally into the wrong argument.

Watch the number of independent groups, not the number of folds #

Ten folds over a survey that spans only fourteen 5 km blocks with samples in them is not ten-fold validation in any useful sense: fold-to-fold variance explodes and the mean is unstable. Count len(np.unique(groups)) before choosing n_splits, and aim for at least five groups per fold. When the count is genuinely small — a handful of tightly clustered sites — that is the signal to switch to buffered leave-one-out, which spends every sample as a test case exactly once. The mechanics of the group construction are in Spatial K-Fold Cross-Validation Setup.

Troubleshooting #

Next Steps #

With the design chosen, implement it properly: Spatial Block Cross-Validation in Python for the mapping case, Buffered Leave-One-Out Cross-Validation when samples are few and clustered, and the rest of the Cross-Validation Strategies topic for the transfer case.

Frequently Asked Questions #

Is plain k-fold ever acceptable on spatial data?

Yes, whenever the model will be applied at the same density it was trained at. A sensor network that loses a station and needs the gap filled from its neighbours faces exactly the situation random k-fold simulates: a test location surrounded by training data at typical inter-sample spacing. Blocking that problem is not more rigorous, it is answering a question nobody asked, and it will make you reject a model that would have performed well. The test is the deployment gap, not the presence of autocorrelation.

How large should a spatial block be?

The median distance from a test point to its nearest training sample comes out at roughly half the block side, so set the block side to about twice the gap the model will face in deployment. Then check that side against the residual variogram range: if it is much shorter than the range, folds are still exchanging correlated information and the score will be optimistic. In the worked example a 2.58 km deployment gap and a 4.8 km range both point at a 5 km block.

Why does buffered leave-one-out report a lower score than spatial blocking here?

Because a 5 km buffer removes every sample within 5 km of the held-out point, the median distance from test to nearest training sample is 5.28 km, more than twice the 2.41 km that 5 km blocking produces. It is a harder split, and it reports a correspondingly lower score of 0.46 against 0.51. That is not evidence that blocking is optimistic; it is evidence that the two designs simulate different deployment gaps.

Should I pick the strategy that gives the best score or the worst?

Neither. Pick the strategy before you see any scores, from the deployment description alone, and report what it gives you. A model validated on a harder split than it will face is under-sold and may be rejected for no reason; a model validated on an easier split is dangerous, because the number attached to it will not survive contact with the map. Choosing after the fact turns the validation design into a free parameter tuned against the thing it was meant to measure.

Related

• Spatial Block Cross-Validation in Python — the mapping case implemented in full, including block origin offsets
• Buffered Leave-One-Out Cross-Validation — the design for few, tightly clustered samples
• Environmental Stratification Cross-Validation in Python — grouping on covariate space to test transfer to a new region

← Back to Cross-Validation Strategies