Designing Spatial Cross-Validation for Raster Models

Size the blocks from the data, not from habit: estimate how far the target stays correlated, then make the blocks bigger than that.

import numpy as np

def empirical_variogram(x, y, values, *, bins=np.arange(0, 40_000, 2_000), n=5_000):
    """Semivariance against separation distance for a sample of point pairs."""
    rng = np.random.default_rng(0)
    i, j = rng.integers(0, len(values), n), rng.integers(0, len(values), n)
    d = np.hypot(x[i] - x[j], y[i] - y[j])
    gamma = 0.5 * (values[i] - values[j]) ** 2
    idx = np.digitize(d, bins)
    return bins, np.array([gamma[idx == k].mean() if (idx == k).any() else np.nan
                           for k in range(1, len(bins) + 1)])

The distance at which that curve flattens is the correlation range, and any block smaller than it leaks. This page belongs to validating raster predictions and accuracy assessment in Raster Machine Learning & Model Inference.

Reading the Range off a Variogram

The empirical variogram is the cheapest possible diagnostic for how much structure a target carries, and it answers the block-size question directly rather than by argument.

Semivariance against separation distance Semivariance rises steeply over the first few kilometres and flattens to a sill at around fourteen kilometres. Below the range, two samples are partly redundant; above it they are effectively independent. The block edge for cross-validation should be set at or above that range. The range is where the curve stops rising sill 0 range ≈ 14 km block edge: 15 km or more samples here are partly the same observation samples here are effectively independent 0 km 36 km

For a categorical target, run the same computation on an indicator variable — one for the class of interest, zero otherwise. The resulting range describes how far patches of that class extend and cluster, which is exactly the quantity a block needs to exceed.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Variogram sampling and fold arithmetic
scikit-learn >=1.3 GroupKFold, LeaveOneGroupOut, scoring
geopandas >=0.14 Region polygons and spatial joins
pandas >=2.0 The sample table carrying coordinates and folds 
pip install "numpy>=1.23" "scikit-learn>=1.3" "geopandas>=0.14" "pandas>=2.0"

Complete Working Example

import numpy as np
import pandas as pd
from sklearn.base import clone
from sklearn.metrics import f1_score
from sklearn.model_selection import GroupKFold


def spatial_cv(model, X: np.ndarray, y: np.ndarray, xy: np.ndarray,
               *, block_m: float = 15_000, n_folds: int = 5,
               buffer_m: float = 0.0) -> pd.DataFrame:
    """Blocked cross-validation with an optional discarded buffer between folds."""
    bx = np.floor(xy[:, 0] / block_m).astype(int)
    by = np.floor(xy[:, 1] / block_m).astype(int)
    block = bx * 1_000_003 + by

    rows = []
    for fold, (train, test) in enumerate(
            GroupKFold(n_splits=n_folds).split(X, y, groups=block)):

        if buffer_m > 0:
            # Drop training samples within buffer_m of any test sample's block
            test_blocks = set(block[test].tolist())
            near = np.array([
                (bx[i] * 1_000_003 + by[i]) in test_blocks
                or any((bx[i] + dx) * 1_000_003 + (by[i] + dy) in test_blocks
                       for dx in (-1, 0, 1) for dy in (-1, 0, 1))
                for i in train
            ])
            train = train[~near]

        fitted = clone(model).fit(X[train], y[train])
        pred = fitted.predict(X[test])
        rows.append({
            "fold": fold,
            "n_train": len(train),
            "n_test": len(test),
            "macro_f1": f1_score(y[test], pred, average="macro"),
            "accuracy": (pred == y[test]).mean(),
        })

    out = pd.DataFrame(rows)
    print(out.to_string(index=False))
    print(f"\nmacro F1: {out.macro_f1.mean():.3f} "
          f"(range {out.macro_f1.min():.3f}-{out.macro_f1.max():.3f})")
    return out

The buffer implementation deliberately drops from the training side rather than the test side. Shrinking the training set slightly is a cost paid in model quality; shrinking the test set would change what is being measured and make folds incomparable.

Interpreting the Spread

Two models, one mean, different products Both models average 0.80 macro F1 across five spatial folds. The first scores between 0.78 and 0.82 everywhere and can be deployed with confidence. The second ranges from 0.62 to 0.94, meaning its behaviour depends heavily on where it is applied, and a single mean would hide that entirely. Both average 0.80 macro F1 1.0 0.5 model A — spread 0.04 deploy it model B — spread 0.32 find out what the bad folds have in common

A wide spread is information, not noise. Map the fold scores back onto the study area and look at the low-scoring blocks: they usually share something identifiable — a different ecological zone, a persistent cloud regime, a season with few observations, a region where the labels were collected by a different team. Each of those has a different remedy, and none of them is visible from the mean.

This is also the check that decides whether one model or several is the right answer. If the low fold corresponds to a genuinely different landscape, a separate model for that landscape will usually beat any amount of tuning on a single global one — and the fold scores are the evidence for that decision rather than an intuition about it.

Verification

import numpy as np

# No block appears in both sides of any split
for fold, (train, test) in enumerate(GroupKFold(5).split(X, y, groups=block)):
    assert not (set(block[train]) & set(block[test])), f"fold {fold} leaks blocks"

# Folds are comparable in size and class mix
sizes = [len(test) for _, test in GroupKFold(5).split(X, y, groups=block)]
assert max(sizes) / min(sizes) < 2.0, f"fold sizes very unbalanced: {sizes}"
Map the fold scores, do not just average them The five folds are drawn in their geographic positions and shaded by score. The two low-scoring folds are adjacent and both cover the upland part of the study area, which points at a landscape difference rather than random variation and suggests either more labels there or a separate model. Fold scores in their geographic positions 0.84 0.82 0.86 0.63 — uplands 0.66 — uplands Two adjacent weak folds is a pattern; two scattered weak folds is variance.

The size check is the one that catches a design problem rather than a code problem. Blocks vary in how many samples they contain, and if one block holds a third of the data no grouping scheme will give even folds; the fix is a better sampling design, not a different splitter. The general fold-assignment mechanics are covered in splitting chips into spatially disjoint folds.

Common Errors

The variogram is flat from the first bin

Either the target has no spatial structure at the sampled distances, or the sample is too small. Reduce the bin width and resample; a genuinely flat variogram means a random split would have been fine, which is rare but not impossible.

Buffering removes most of the training data

The buffer is large relative to the block. Keep the buffer to one block ring, and enlarge the blocks rather than the buffer if leakage is still suspected.

Scores vary wildly between runs with the same settings

Blocks are being assigned to folds in geographic order, so the split depends on iteration order. Shuffle the block list with a fixed seed before assigning.

Frequently Asked Questions

Q: How do I find the autocorrelation range in practice? Compute an empirical variogram of the target variable on a sample of point pairs and find the distance at which the semivariance flattens. For categorical targets use an indicator variogram on the class of interest. Round the answer up and use it as the block edge.

Q: Is the mean of the fold scores the number to report? Report both the mean and the range. A model scoring 0.80 on every fold and one averaging 0.80 across folds of 0.62 and 0.94 are very different products, and only the second warns you that performance depends heavily on where it is applied.

Q: Does spatial cross-validation apply to time as well? The same logic applies to any dimension with autocorrelation. For a model that will be applied to a future year, hold out whole years rather than random dates, or the score measures interpolation within a season you have already seen.