Geostatistical Modeling Spatial Statistics with Python

Spatial Block Cross-Validation in Python

TL;DR: Assign each observation to a spatial block — either a regular grid tile from (x // block_size, y // block_size) or a KMeans grouping on the coordinate array — then feed those block labels as groups to sklearn.model_selection.GroupKFold. Because whole blocks move together into train or test, no held-out point has a training neighbour a few metres away, and the resulting RMSE is an honest estimate of performance on new terrain, typically well above the optimistic random-CV number.

Why This Matters

Standard random K-fold cross-validation silently lies about spatial models. Because spatial autocorrelation makes nearby observations similar, random folds place a test point’s near-neighbours in the training set, so the model interpolates between near-duplicates and reports a flatteringly low error. Deploy that model on a genuinely unsampled region and the true error can be several times larger. Spatial block cross-validation fixes this by holding out contiguous geographic blocks, forcing the model to predict across a gap. This page is a practical companion to the Cross-Validation Strategies guide within the broader Python Workflows for Spatial Modeling & Regression stack.

The single most useful output here is the gap between random and blocked RMSE. A large gap is a red flag that your reported accuracy is inflated by leakage; a small gap suggests the model genuinely generalises across space. Measuring that gap should be a standard step before any spatial model ships.

Environment

bash
pip install \
  geopandas==0.14.4 \
  scikit-learn==1.4.2 \
  numpy==1.26.4 \
  pandas==2.2.1
python
import numpy as np
import pandas as pd
import geopandas as gpd
from sklearn.cluster import KMeans
from sklearn.model_selection import GroupKFold, KFold
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error

Assume a projected (metric) GeoDataFrame gdf with a target column y and feature columns in feat_cols. A metric CRS is mandatory so block sizes are in metres, not degrees.

Step-by-Step Implementation

Step 1 — Extract Coordinates and Build Grid Blocks

The simplest blocking scheme tiles the coordinate plane. Integer-divide each coordinate by a chosen block edge length and combine the two indices into a block id.

python
coords = np.column_stack([gdf.geometry.x.values, gdf.geometry.y.values])

def grid_blocks(coords: np.ndarray, block_size: float) -> np.ndarray:
    """Assign each point to a square grid tile of edge `block_size` (map units)."""
    ix = np.floor((coords[:, 0] - coords[:, 0].min()) / block_size).astype(int)
    iy = np.floor((coords[:, 1] - coords[:, 1].min()) / block_size).astype(int)
    # Combine tile indices into a single label
    return ix * (iy.max() + 1) + iy

block_id = grid_blocks(coords, block_size=25_000)   # 25 km tiles
print(f"{len(np.unique(block_id))} occupied blocks")

Choose block_size at or above the variogram range so autocorrelation does not straddle the block boundary.

Step 2 — Alternative: K-Means Blocks on Coordinates

When sampling is spatially uneven, grid tiles hold wildly different point counts. Clustering the coordinates with k-means yields compact blocks of similar size.

python
def kmeans_blocks(coords: np.ndarray, n_blocks: int, seed: int = 42) -> np.ndarray:
    """Cluster coordinates into `n_blocks` compact spatial blocks."""
    km = KMeans(n_clusters=n_blocks, random_state=seed, n_init=10)
    return km.fit_predict(coords)

block_id = kmeans_blocks(coords, n_blocks=20)

Each k-means block is a contiguous Voronoi region, so held-out blocks are still spatially coherent, just adaptively sized.

Step 3 — Fold Blocks into GroupKFold

GroupKFold guarantees that all rows sharing a group label land in the same fold. Passing the block id as the group makes whole blocks the unit of held-out data.

python
def blocked_cv_rmse(X, y, groups, n_splits=5, seed=0):
    """Return per-fold RMSE using spatial blocks as GroupKFold groups."""
    gkf = GroupKFold(n_splits=n_splits)
    rmses = []
    for train_idx, test_idx in gkf.split(X, y, groups=groups):
        model = RandomForestRegressor(n_estimators=300, random_state=seed, n_jobs=-1)
        model.fit(X[train_idx], y[train_idx])
        pred = model.predict(X[test_idx])
        rmses.append(np.sqrt(mean_squared_error(y[test_idx], pred)))
    return np.array(rmses)

X = gdf[feat_cols].to_numpy()
y = gdf["y"].to_numpy()
blocked_rmse = blocked_cv_rmse(X, y, groups=block_id, n_splits=5)
print(f"Blocked RMSE: {blocked_rmse.mean():.3f} +/- {blocked_rmse.std():.3f}")

Step 4 — Run Random K-Fold for Comparison

Score the identical model with ordinary random folds so the only thing that changes is how points are partitioned.

python
def random_cv_rmse(X, y, n_splits=5, seed=0):
    """Return per-fold RMSE using random (non-spatial) K-fold."""
    kf = KFold(n_splits=n_splits, shuffle=True, random_state=seed)
    rmses = []
    for train_idx, test_idx in kf.split(X):
        model = RandomForestRegressor(n_estimators=300, random_state=seed, n_jobs=-1)
        model.fit(X[train_idx], y[train_idx])
        pred = model.predict(X[test_idx])
        rmses.append(np.sqrt(mean_squared_error(y[test_idx], pred)))
    return np.array(rmses)

random_rmse = random_cv_rmse(X, y, n_splits=5)
print(f"Random RMSE:  {random_rmse.mean():.3f} +/- {random_rmse.std():.3f}")

Step 5 — Quantify the Optimism Gap

Report both numbers and their ratio. The ratio is the multiplier by which random CV understates real-world error.

python
gap = blocked_rmse.mean() / random_rmse.mean()
print(f"Random RMSE:  {random_rmse.mean():.3f}")
print(f"Blocked RMSE: {blocked_rmse.mean():.3f}")
print(f"Optimism ratio (blocked / random): {gap:.2f}x")

Interpreting the Output

A typical result on autocorrelated data looks like a random RMSE of, say, 4.1 against a blocked RMSE of 6.8 — an optimism ratio near 1.7x. The blocked figure is the number to trust and to publish, because it reflects prediction into space the model has not seen. A ratio close to 1.0 means either the target is weakly autocorrelated or your blocks are smaller than the correlation range and leakage persists; enlarge the blocks and re-check. A ratio above 2x is common for dense, strongly clustered sampling and is a strong signal that random CV would have badly misled deployment expectations.

Also inspect the spread across folds. High between-fold variance under blocking is expected and honest: some regions are genuinely harder to predict than others, and blocking surfaces that heterogeneity instead of averaging it away.

Random folds versus spatial block folds Left panel shows points coloured by fold scattered randomly so every test point sits beside training points. Right panel shows points grouped into contiguous square blocks, each block a single fold, so held-out points are separated from training points by a spatial gap. Random K-Fold Spatial Block CV Shade = fold assignment. Blocks keep each fold spatially contiguous.

Critical Best Practices

Size Blocks to the Variogram Range

The block edge must equal or exceed the range at which autocorrelation vanishes. Estimate it first with an empirical variogram (see the Variogram Modeling & Semivariance Analysis pillar) and set block_size to that range or larger. Blocks smaller than the range leave correlated pairs straddling the boundary and quietly restore the optimism you meant to remove.

Keep Enough Blocks for Stable Folds

Five folds need many more than five blocks, or a single block dominates a fold and RMSE becomes noisy. Aim for at least four to eight blocks per fold so each held-out fold samples several regions. With k-means blocking, set n_blocks to roughly n_splits * 5 as a starting point.

Never Fit Preprocessing on the Full Dataset

Scalers, imputers, and feature selectors must be fit inside the fold on training blocks only. Fitting them on all data before splitting leaks held-out statistics across the block boundary and reinflates the score. Wrap preprocessing and model in a Pipeline and pass that to the CV loop.

Match the Blocking Scheme to Sampling Geometry

Uniform sampling suits grid blocks; clustered sampling suits k-means blocks. A grid over clustered data produces empty and overcrowded tiles that unbalance folds. Plot the block assignment before trusting the scores.

Troubleshooting

Symptom Likely cause Fix
Blocked RMSE barely above random Blocks smaller than variogram range Increase block_size to at least the range
Wildly unequal fold sizes Grid tiles over clustered sampling Switch to kmeans_blocks for balanced groups
GroupKFold raises on n_splits > n_groups Fewer occupied blocks than folds Reduce n_splits or create more blocks
Scores unstable across reruns Too few blocks per fold Raise the block count; fix random_state
Blocked RMSE lower than random Preprocessing fit on full data (leakage) Fit all transforms inside the fold via a Pipeline
Block ids collide across tiles Index-combination overflow Use pandas.factorize on (ix, iy) tuples instead of arithmetic

Next Steps

Spatial block CV is one of a family; for the point-level variant that removes a buffer around each held-out observation, see Buffered Leave-One-Out Cross-Validation, and for the general fold-construction mechanics see Spatial K-Fold Cross-Validation Setup.

Related:

← Back to Cross-Validation Strategies

Related

Buffered Leave-One-Out Cross-Validation Spatial K-Fold Cross-Validation Setup in Python