⚙️ How-To Spatial machine learning Spatial Machine Learning · step 6 of 14 Updated

How to Do Spatial Cross-Validation in Python

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 a vegetation index from terrain, 3,000 training cells, repeated over ten random samples:

                        mean R²    sd
random 5-fold CV         +0.555   0.019
spatial block CV         +0.213   0.106
true error over the
whole unsampled map      +0.557   0.014

Random CV was accurate to within 0.002. Block CV understated the model by 0.34 R² and was five times more variable.

Change the task to genuine extrapolation — train on the western half, predict the east — and both lie:

random CV      +0.649
block CV       +0.385
actual          -0.771

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
train = region != "east"
Three validation schemes matched to three deployment scenarios: interpolation, unsampled areas, and a new region.
There is no universally honest scheme. Each estimates a different quantity.

Step-by-step solution

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.002 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 range was several kilometres, so five blocks over a 10 km window was, if anything, too few.

3. Use enough blocks, and repeat

Five blocks gave a standard deviation of 0.106 across ten block layouts, ranging from +0.061 to +0.429. A single number from a single layout is false precision.

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.pipeline import make_pipeline

model = make_pipeline(StandardScaler(), RandomForestRegressor())
cross_val_score(model, X, y, cv=folds)

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.55 is not a result. R² = 0.55, random 5-fold, for prediction within the sampled extent is.

Block CV scores ranging from 0.061 to 0.429 across ten block layouts while the true error stays near 0.557.
Block CV has five times the variance of random CV. One layout is not a measurement.

Code examples

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

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

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: a 500 m buffer on 10 km blocks removes a noticeable fraction of the training data, which makes the estimate more pessimistic still.

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
  left -> right: R2 -0.771, MAE 0.1523 (3,000 train, 30,000 test)

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.

Explanation

Why block CV was pessimistic here

Block CV asks a harder question than the deployment does. Each fold trains on 80% of the area and predicts a large contiguous region with no nearby training data.

If the model will actually predict at locations surrounded by training data — filling gaps in a survey — that is not the question. Here it understated the model by 0.34 R², enough to make a working model look like a failure.

It is also structurally noisy: with five folds, the score averages five hard, correlated problems rather than many independent ones. Hence 0.106 standard deviation against random CV's 0.019.

Why neither scheme handles a new region

The west-to-east test is the case everyone worries about, and block CV did not save it: it promised +0.385 where reality was −0.771.

Block CV resamples within the training region. If the relationship between covariates and target differs in the new region — different geology, land use, climate — nothing inside the training data contains that information.

Cross-validation estimates generalisation to data from the same distribution. A new region is a different distribution, and the honest answer is "we do not know", supported by an area-of-applicability mask.

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

KMeans on coordinates produces different blocks for different seeds, and the score depends on which regions each fold happens to hold out.

Measured: ten layouts gave +0.061 to +0.429, a range of 0.37 R² on identical data 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.

Blocks smaller than the autocorrelation range leaking across boundaries against larger buffered blocks.
A block smaller than the range produces a score between random and true block CV, with nothing to say which.

Edge cases or notes

  • 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 block layouts. One layout has a standard deviation of 0.106 here.
  • Use a Pipeline so 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.

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.002 R² while block CV was pessimistic by 0.34.

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?

Ten to twenty, with several repeats. Five blocks gave a standard deviation of 0.106 across layouts 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.385 where a real west-to-east transfer gave −0.771. 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.