Spatial Leakage Explained: Why Random Cross-Validation Lies

Problem statement

The standard advice about machine learning on spatial data is: random cross-validation is optimistic because nearby points are correlated, so use spatial block cross-validation instead.

That advice is half right, and the missing half causes as much damage as the problem it solves.

Measured on a real modelling problem — predicting Sentinel-2 NDVI from terrain over a 10 km window 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

Random CV was accurate to within 0.005. Block CV understated the model by 0.24 R² and its spread was more than four times random CV's. Following the standard advice here would have led you to reject a model that works.

Change the deployment scenario, though, and the picture inverts completely — see below. Neither scheme is "the honest one". The scheme has to match what you will actually do with the model.

Quick answer

Ask what the model will be asked to predict, and validate that:

from sklearn.model_selection import KFold, GroupKFold
from sklearn.cluster import KMeans

# interpolation within the sampled area
folds = KFold(n_splits=5, shuffle=True, random_state=0)

# extrapolation to unsampled ground
blocks = KMeans(n_clusters=5, random_state=0, n_init=10).fit_predict(coords)
folds = GroupKFold(n_splits=5)          # groups=blocks

What leakage is, measured for interpolation and extrapolation

1. Recognise what leakage actually is

Leakage is information reaching the test set that will not be available at prediction time.

With spatially autocorrelated data and a random split, a test point often has a training point twenty metres away. The model can effectively memorise the neighbourhood.

That is only leakage if your real predictions will be far from any training point. If you will predict at locations surrounded by training data — filling gaps in a survey, densifying a monitoring network — then having a near neighbour is not leakage. It is the deployment condition.

2. Measure the interpolation case

Training on 3,000 randomly placed cells and predicting the remaining 267,536:

random 5-fold CV   +0.339 (sd 0.024 over 10 samples)
block CV           +0.102 (sd 0.106)
true map R²        +0.344 (sd 0.008)

Random CV reproduced the truth almost exactly. Block CV was pessimistic by 0.24 and unstable: across ten runs it ranged from −0.057 to +0.258, while the truth barely moved.

The instability has a mechanical cause. With five blocks, each fold trains on 80% of the area and predicts a large contiguous region, and how hard that is depends heavily on which region each block covers. The blocks follow the sample: re-seeding the clustering on one sample of 3,000 cells hardly changed them or the score (sd 0.005, Example 2), but a different sample of 3,000 cells made different blocks and a different answer.

3. Measure the extrapolation case

Now train on the west half only and predict the east half — a genuinely unsampled region:

random 5-fold CV on the training half   +0.338
block CV on the training half           +0.033
actual R² on the eastern half           -0.409

Both lie, in the same direction. Random CV promised +0.338; block CV, the supposedly honest scheme, still promised +0.033. Reality was −0.409.

A negative R² means the model is worse than predicting the training mean everywhere. The western half is lower and greener than the eastern — median elevation 348 m against 421 m, mean NDVI 0.747 against 0.713 — and the terrain-to-vegetation relationship learnt in one half does not hold in the other. No resampling of the western data could reveal that.

4. Draw the actual conclusion

Block CV is not a truth serum. It is a different estimator with a different bias:

• For interpolation, it is strongly pessimistic (0.102 against 0.344) and noisy.
• For extrapolation, it is less optimistic than random CV but still optimistic (0.033 against −0.409).

The only way to estimate transfer to a new region is to hold out a region of the kind you will actually predict, and accept that you cannot do this at all if you have only sampled one region.

5. Report the scheme with the number

R² = 0.34 is not a result. R² = 0.34 under random 5-fold CV, for prediction within the sampled extent is.

Random against block CV, scored against the truth

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 — score the model both ways and against the truth

import numpy as np
from sklearn.cluster import KMeans
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from sklearn.model_selection import GroupKFold, KFold


def compare_schemes(X, y, coords, X_holdout=None, y_holdout=None,
                    n_blocks=5, seed=0):
    """Random CV, block CV, and — if you have it — the truth."""
    def model():
        return RandomForestRegressor(n_estimators=200, min_samples_leaf=2,
                                     n_jobs=-1, random_state=seed)

    def cv(splitter):
        predicted = np.empty(len(y))
        for train, test in splitter:
            predicted[test] = model().fit(X[train], y[train]).predict(X[test])
        return r2_score(y, predicted)

    random_r2 = cv(KFold(n_splits=5, shuffle=True, random_state=seed).split(X))

    blocks = KMeans(n_clusters=n_blocks, random_state=seed,
                    n_init=10).fit_predict(coords)
    block_r2 = cv(GroupKFold(n_splits=n_blocks).split(X, y, blocks))

    print(f"  random 5-fold CV   R2 {random_r2:+.3f}")
    print(f"  {n_blocks}-block spatial CV R2 {block_r2:+.3f}")

    result = {"random": random_r2, "block": block_r2}
    if X_holdout is not None:
        true_r2 = r2_score(y_holdout, model().fit(X, y).predict(X_holdout))
        print(f"  true held-out      R2 {true_r2:+.3f}")
        result["true"] = true_r2
    return result


rng = np.random.default_rng(0)
sample = rng.choice(len(y), 3000, replace=False)
rest = np.setdiff1d(np.arange(len(y)), sample)
Xs, ys, cs = X[sample], y[sample], coords[sample]
compare_schemes(Xs, ys, cs, X_holdout=X[rest], y_holdout=y[rest])

  random 5-fold CV   R2 +0.374
  5-block spatial CV R2 +0.161
  true held-out      R2 +0.348

Where you can construct a genuine hold-out — a dense reference dataset, a pilot area, a later survey — do it once. It calibrates how much to trust each CV scheme for your variable and landscape, and that calibration transfers to future projects far better than any rule of thumb.

Example 2 — check how stable your block CV is

import numpy as np
from sklearn.cluster import KMeans


def cv_with_groups(X, y, groups):
    predicted = np.empty(len(y))
    for train, test in GroupKFold(n_splits=len(np.unique(groups))).split(X, y, groups):
        predicted[test] = RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1,
                                                random_state=0).fit(X[train], y[train]).predict(X[test])
    return r2_score(y, predicted)


def block_cv_stability(X, y, coords, n_repeats=10, n_blocks=5):
    """Block CV depends on which blocks the clustering happened to make."""
    scores = []
    for seed in range(n_repeats):
        blocks = KMeans(n_clusters=n_blocks, random_state=seed,
                        n_init=10).fit_predict(coords)
        scores.append(cv_with_groups(X, y, blocks))

    scores = np.array(scores)
    print(f"  block CV over {n_repeats} block layouts: "
          f"mean {scores.mean():+.3f}, sd {scores.std():.3f}, "
          f"range {scores.min():+.3f} to {scores.max():+.3f}")
    if scores.std() > 0.05:
        print("  ! this estimate is unstable — report the spread, not one number")
    return scores


block_cv_stability(Xs, ys, cs)

  block CV over 10 block layouts: mean +0.161, sd 0.005, range +0.145 to +0.165

Re-seeding the clustering on one sample barely moved the score here, because the blocks are decided by where the sample fell. Across ten different samples of 3,000 cells, block CV had a standard deviation of 0.106 (Example 4) — that is the uncertainty a single block-CV number carries, and quoting it to three decimal places, as everyone does, is false precision by a wide margin.

Example 3 — the only honest extrapolation test

import numpy as np
from sklearn.metrics import r2_score


def region_transfer_test(X, y, coords, split_axis=0):
    """Train on one half of the study area, predict the other. No resampling."""
    threshold = np.median(coords[:, split_axis])
    left = coords[:, split_axis] < threshold

    model = RandomForestRegressor(n_estimators=200, min_samples_leaf=2,
                                  n_jobs=-1, random_state=0)
    forward = r2_score(y[~left], model.fit(X[left], y[left]).predict(X[~left]))
    backward = r2_score(y[left], model.fit(X[~left], y[~left]).predict(X[left]))

    print(f"  train left  -> predict right: R2 {forward:+.3f}")
    print(f"  train right -> predict left : R2 {backward:+.3f}")
    if min(forward, backward) < 0:
        print("  ! the model does not transfer — it is worse than the mean "
              "in at least one direction")
    return {"forward": forward, "backward": backward}


six = rng.choice(len(y), 6000, replace=False)            # 3,000 cells either side of the median x
region_transfer_test(X[six], y[six], coords[six])

  train left  -> predict right: R2 -0.109
  train right -> predict left : R2 +0.111
  ! the model does not transfer — it is worse than the mean in at least one direction

Run both directions. A model that transfers one way and not the other tells you the two halves differ, which is more useful than a single number and is often the actual scientific finding.

Example 4 — the tables at the top of this page

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):
        predicted = np.empty(len(ys))
        for train, test in splitter:
            predicted[test] = model().fit(Xs[train], ys[train]).predict(Xs[test])
        return r2_score(ys, predicted)
    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)


# Interpolation: 3,000 random cells, predict every other cell, ten samples
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]))))
runs = np.array(runs)
print(f"3,000 training cells, {len(y) - 3000:,} predicted, 10 samples")
for name, column in zip(("random 5-fold CV", "5-block spatial CV", "true map R²"), runs.T):
    print(f"  {name:20} mean {column.mean():+.3f}  sd {column.std():.3f}  "
          f"range {column.min():+.3f} to {column.max():+.3f}")

# Extrapolation: 3,000 cells from the western half, predict the whole eastern half
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}, "
      f"actual on the {int((~west).sum()):,} eastern cells "
      f"{r2_score(y[~west], model().fit(X[pick], y[pick]).predict(X[~west])):+.3f}")
for label, half in (("west", west), ("east", ~west)):
    print(f"  {label}: elevation median {np.median(X[half, 0]):.0f} m, NDVI mean {y[half].mean():.3f}")

3,000 training cells, 267,536 predicted, 10 samples
  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 on the 135,292 eastern cells -0.409
  west: elevation median 348 m, NDVI mean 0.747
  east: elevation median 421 m, NDVI mean 0.713

Why random CV was right here and block CV pessimistic

Why random CV was accurate here

The 3,000 training cells were spread randomly over the whole window, and the model was then asked to predict the other 267,536 cells in the same window.

A randomly held-out cell is, in expectation, exactly as far from the training set as a randomly chosen unpredicted cell. The CV task and the deployment task are the same task, so the CV estimate is unbiased — which is what the measurement shows: +0.339 against +0.344.

Why block CV was pessimistic

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

If you will never do that, the number is answering a question you did not ask. Here it understated the model by 0.24 R² — enough to make a working model look like a failure.

It is also noisy, for a structural reason: with five folds, each test set is one large region, so the score is an average of five hard, correlated problems rather than of many independent ones. Hence sd 0.106 against random CV's 0.024.

Why both schemes failed on the new region

The west-to-east test is the case everyone is worried about, and it is worth noting that block CV did not save it. It promised +0.033 where reality was −0.409.

The reason is that block CV resamples within the training region. If the relationship between terrain and vegetation is different in the east — different geology, different land use, different rainfall — nothing inside the western data contains that information. No cross-validation scheme can estimate an error it has no evidence about.

This is the practical limit of validation. Cross-validation measures how well a model generalises to data drawn from the same distribution as the training data. A new region is a different distribution, and the honest answer is "we do not know", supported by an area-of-applicability mask rather than a number.

Why a negative R² is worth stating plainly

R² = −0.409 means the model's squared error is 1.41 times the variance of the target in the new region — it would have been better to predict the training mean everywhere and not build a model at all.

Negative R² is common in genuine transfer tests and almost never reported, because most pipelines never run one. When you see it, the model has not merely degraded; it has learned a relationship that is actively wrong somewhere else.

Block-CV variance, near-duplicates and negative R²

• State the scheme with every score. A bare R² is not interpretable.
• Block CV has a large variance. Repeat over several block layouts and report the spread.
• The block size matters as much as the scheme. Blocks should be at least the autocorrelation range, or the folds still leak.
• Block CV scores are not comparable across sample sizes, because the block geometry changes with the sample.
• Autocorrelation is not the only leaking dimension. Time, sensor, operator and processing batch all leak.
• Duplicate or near-duplicate rows — two samples from one field — leak under any scheme.
• A negative R² is meaningful, not a bug. It means the model is worse than the mean.
• Do not tune hyperparameters on the same folds you report. That leaks in the ordinary, non-spatial way.

Internal links

• Spatial machine learning explained — why location breaks the usual assumptions
• How to do spatial cross-validation in Python — the implementation
• Extrapolation in space explained: the area of applicability — what to do when you cannot validate transfer
• 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 — the structure that causes it
• Feature importance says coordinates are the best predictor — the related coordinate trap
• How to evaluate a spatial model honestly — reporting beyond one number

FAQ

Is random cross-validation always wrong for spatial data?

No. When the model will predict at unsampled locations inside a well-sampled area, random CV matched the truth almost exactly here — +0.339 against +0.344.

When should I use spatial block cross-validation?

When you will predict in areas with no nearby training data. For interpolation within a sampled area it was pessimistic by 0.24 R² and more than four times noisier.

Does block cross-validation give the true error?

No. On a genuine new-region test it promised +0.033 where the actual result was −0.409. It is less optimistic than random CV, not honest.

How do I estimate error in a region I have not sampled?

Hold out a real region of the same kind and test on it. If you have only sampled one region, you cannot estimate it — use an area-of-applicability mask and say so.

Why is my block CV score so variable?

Because five folds means five large, correlated test regions, and the score depends on which regions the blocks cover — which follows from where the sample fell. It had a standard deviation of 0.106 across ten samples here.

What does a negative R² mean?

The model is worse than predicting the mean. It happens routinely in genuine transfer tests and is worth reporting rather than hiding.

How big should the blocks be?

At least the autocorrelation range of the target, or neighbouring folds still share information. Estimate it from a variogram rather than guessing.