TL;DR Answer
Set the block edge length in spatial block cross-validation to at least the range of spatial autocorrelation in your target variable. Estimate that range from an empirical variogram (skgstat) or a Moran’s correlogram (esda), round it up, and tile the study area into square blocks of that size. Assign whole blocks — not individual points — to folds, so every test block is separated from training data by one autocorrelation range and no signal leaks across the split. This is a core technique in Spatial Cross-Validation Strategies for producing honest error estimates on geospatial data.
import numpy as np
from skgstat import Variogram
# coords: (n, 2) projected metres; z: target values
V = Variogram(coords, z, model="spherical", n_lags=20)
autocorr_range = V.parameters[0] # range in metres
block_size = float(np.ceil(autocorr_range / 100) * 100)
print(f"range = {autocorr_range:.0f} m -> block_size = {block_size:.0f} m")Part of: Spatial Cross-Validation Strategies
Why Block Size Is the Whole Ballgame
Spatial block cross-validation tiles the study area into a grid and assigns each square block to a fold, so that geographically contiguous observations are held out together. It is the grid-based cousin of the K-means approach in Implementing SpatialKFold in Python, and it fixes the same problem: random folds leak information between nearby training and test points. But blocking only works if the blocks are the right size. Get the size wrong and you either keep the leakage you were trying to remove, or you throw away so much statistical power that your error bars become meaningless.
The governing quantity is the range of spatial autocorrelation — the distance beyond which two observations are effectively independent. Inside the range, values covary because they share environmental drivers the model never sees directly: soil parent material, microclimate, hydrology, land-use history. If a block is narrower than the range, a point near one edge of a test block is still correlated with a point just across the boundary in a training block. Because those two blocks live in different folds, the model gets to peek at the answer. The cross-validation score looks great and collapses in production.
If a block is far larger than the range, you have solved leakage but created a new problem. A large block covers a big chunk of the map, so the whole study area contains only a handful of blocks. Each fold now rests on one or two contiguous regions, and any regional peculiarity — a wetter basin, a different sensor swath — dominates that fold’s score. Fold-to-fold variance explodes, and you may not even have enough blocks to fill the number of folds you wanted. The right block size is the smallest one that reliably exceeds the range.
Estimating the Autocorrelation Range
There are two production-grade ways to measure the range, and it is worth computing both and reconciling them.
Empirical variogram with skgstat
A variogram plots semivariance — half the mean squared difference between value pairs — against the lag distance separating them. Semivariance rises with distance and then flattens; the lag at which it levels off is the range, and the plateau height is the sill. Fit a skgstat.Variogram on projected metric coordinates and read the fitted range parameter directly.
import numpy as np
import geopandas as gpd
from skgstat import Variogram
gdf = gpd.read_file("soil_samples.geojson").to_crs("EPSG:32633")
coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
z = gdf["organic_carbon"].to_numpy()
V = Variogram(coords, z, model="spherical", n_lags=25, maxlag="median")
autocorr_range = V.parameters[0] # (range, sill, nugget)
print(f"Estimated range: {autocorr_range:,.0f} m")Choose the model to match the shape of the empirical points — "spherical" and "exponential" are the usual candidates. For an exponential model, the fitted parameter is the effective range only after scaling (roughly three times the nominal range parameter), so read the skgstat documentation for your chosen model rather than assuming the first parameter is the practical range.
Moran’s correlogram as a cross-check
A Moran’s correlogram computes global Moran’s I at a sequence of distance bands and reports where it decays to near zero — another estimate of the range. This ties directly into diagnosing dependence in the first place, covered in Handling Spatial Autocorrelation.
import numpy as np
from libpysal.weights import DistanceBand
from esda.moran import Moran
bands = np.linspace(500, 8000, 16) # metres
for d in bands:
w = DistanceBand(coords, threshold=d, binary=True, silence_warnings=True)
mi = Moran(z, w, permutations=99)
print(f"{d:6.0f} m I={mi.I:+.3f} p={mi.p_sim:.3f}")The distance at which Moran’s I first becomes statistically indistinguishable from zero is your correlogram-based range. If the variogram and correlogram disagree, take the larger of the two — it is the conservative choice for blocking.
From Range to Block Size
The rule is simple: the block edge length must be at least one range. A test block is then surrounded by a training-data buffer no thinner than the range, because the nearest training point lies in an adjacent block whose near edge is a full block width away from the test block’s far edge. In practice, round the range up to a clean number and, if you can afford the blocks, add a small margin.
import numpy as np
def range_to_block_size(autocorr_range: float, margin: float = 1.0,
round_to: float = 100.0) -> float:
"""Translate an autocorrelation range (m) into a square block edge (m).
Parameters
----------
autocorr_range : float
Range of spatial autocorrelation in metres (must be > 0).
margin : float, default 1.0
Safety multiplier on the range. 1.0 = exactly the range;
1.25 adds a 25% buffer at the cost of fewer blocks.
round_to : float, default 100.0
Round the result up to this multiple for clean grid coordinates.
"""
if autocorr_range <= 0:
raise ValueError("autocorr_range must be positive.")
if margin < 1.0:
raise ValueError("margin must be >= 1.0 to avoid leakage.")
raw = autocorr_range * margin
return float(np.ceil(raw / round_to) * round_to)Before committing, sanity-check that the block size still leaves you enough blocks. Divide the width and height of your data’s bounding box by the block size: the product is the maximum number of blocks. If that number is under about 25–30, the folds will be noisy and you should either accept a smaller margin or gather more spatial extent.
Building Blocked Fold Labels
The production function below takes projected coordinates and a block size in metres, tiles the extent into square blocks, and assigns whole blocks to folds. Blocks are shuffled deterministically and dealt round-robin into n_splits folds so that no fold gets a single contiguous strip of the map.
import logging
import numpy as np
logger = logging.getLogger(__name__)
def make_block_folds(
coords: np.ndarray,
block_size: float,
n_splits: int = 5,
random_state: int | None = None,
) -> np.ndarray:
"""Assign each observation to a spatial-block fold.
Parameters
----------
coords : np.ndarray of shape (n_samples, 2)
Projected metric coordinates (easting, northing) in metres.
block_size : float
Square block edge length in metres. Must exceed the range of
spatial autocorrelation to prevent train/test leakage.
n_splits : int, default 5
Number of folds to deal the blocks into.
random_state : int or None, default None
Seed controlling the block-to-fold shuffle.
Returns
-------
np.ndarray of shape (n_samples,)
Integer fold label in [0, n_splits) for every observation.
"""
if coords.ndim != 2 or coords.shape[1] != 2:
raise ValueError("coords must have shape (n_samples, 2).")
if block_size <= 0:
raise ValueError("block_size must be positive.")
if coords.shape[0] < n_splits:
raise ValueError("Need at least n_splits observations.")
origin = coords.min(axis=0)
block_ix = np.floor((coords - origin) / block_size).astype(np.int64)
# Unique block id per (col, row) cell
_, block_of_point = np.unique(block_ix, axis=0, return_inverse=True)
n_blocks = block_of_point.max() + 1
if n_blocks < n_splits:
raise ValueError(
f"Only {n_blocks} blocks for {n_splits} folds — "
"reduce block_size or n_splits."
)
rng = np.random.default_rng(random_state)
shuffled = rng.permutation(n_blocks)
block_fold = np.empty(n_blocks, dtype=np.int64)
block_fold[shuffled] = np.arange(n_blocks) % n_splits
fold_labels = block_fold[block_of_point]
for f in range(n_splits):
logger.debug("Fold %d: %d points", f, int((fold_labels == f).sum()))
return fold_labelsTo plug this into scikit-learn, feed the labels to PredefinedSplit, which yields standard (train_idx, test_idx) pairs usable by cross_val_score, GridSearchCV, and Pipeline:
from sklearn.model_selection import PredefinedSplit
fold_labels = make_block_folds(coords, block_size=block_size,
n_splits=5, random_state=42)
cv = PredefinedSplit(test_fold=fold_labels)How Block Size Controls Leakage
The diagram contrasts a block that is too small against one sized to the range. When the block edge is shorter than the range, the correlation shadow of a test point reaches into the neighbouring training block; sizing the block to the range pushes the nearest training point a full range away.
Walkthrough: Estimate the Range, Then Set the Block Size
Putting the pieces together on a single dataset:
import numpy as np
import geopandas as gpd
from skgstat import Variogram
# 1. Load points in a projected metric CRS
gdf = gpd.read_file("soil_samples.geojson").to_crs("EPSG:32633")
coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
z = gdf["organic_carbon"].to_numpy()
# 2. Estimate the range from the variogram
V = Variogram(coords, z, model="spherical", n_lags=25, maxlag="median")
autocorr_range = V.parameters[0]
print(f"Range: {autocorr_range:,.0f} m") # e.g. 2,300 m
# 3. Translate range -> block size with a modest margin
block_size = range_to_block_size(autocorr_range, margin=1.15, round_to=100)
print(f"Block size: {block_size:,.0f} m") # e.g. 2,700 m
# 4. Confirm block count is adequate
w = coords[:, 0].ptp()
h = coords[:, 1].ptp()
print(f"Approx blocks: {int((w / block_size) * (h / block_size))}")
# 5. Build fold labels
fold_labels = make_block_folds(coords, block_size, n_splits=5, random_state=42)If step 4 reports fewer than about 25 blocks, drop the margin back toward 1.0 or widen the study extent before proceeding. Once the block count is comfortable, the fold labels feed straight into evaluation — pair the resulting scores with Evaluating Model Performance with Spatial Metrics to confirm residuals are no longer spatially structured.
Verification: Minimum Train–Test Distance ≥ Range
The definitive check is geometric. For every fold, the closest training point to any test point must be at least one autocorrelation range away. If it is not, some blocks were undersized or points sit exactly on a shared boundary.
import numpy as np
from scipy.spatial import cKDTree
def audit_block_separation(coords: np.ndarray, fold_labels: np.ndarray,
autocorr_range: float) -> None:
"""Assert every fold's min train-test distance meets the range."""
for f in np.unique(fold_labels):
test = coords[fold_labels == f]
train = coords[fold_labels != f]
d_min = cKDTree(train).query(test, k=1)[0].min()
status = "OK" if d_min >= autocorr_range else "LEAK"
print(f"Fold {f}: min train-test = {d_min:,.0f} m [{status}]")
audit_block_separation(coords, fold_labels, autocorr_range)Every fold should print OK. A LEAK flag means the effective separation fell below the range — usually because test and training blocks touch along an edge, so add a one-block buffer (drop the training blocks immediately bordering each test block) or increase the margin. This distance audit is the same principle used to reduce contamination in Reducing Spatial Leakage in Model Training, applied here at the block boundary rather than the point level.
FAQ
What happens if my block size is smaller than the autocorrelation range?
If a block is narrower than the range of spatial autocorrelation, points in one block still correlate with points in the neighbouring block. Because adjacent blocks land in different folds, training and test samples on either side of a block boundary share the same unmodeled signal, so leakage persists and your cross-validation scores stay optimistically inflated. The block edge length must meet or exceed the range so that any test block is separated from training data by at least one autocorrelation range.
How do I estimate the range of spatial autocorrelation in Python?
Fit an empirical variogram to your target variable using skgstat.Variogram on projected metric coordinates, then read the fitted range parameter — the lag distance at which the variogram flattens to the sill. Alternatively compute a Moran’s correlogram with esda and libpysal and take the distance at which Moran’s I drops to near zero. Both give a distance in metres that you round up to set the block edge length.
Why not just make the blocks as large as possible to be safe?
Oversized blocks trade leakage for variance. Larger blocks mean fewer independent blocks across your study area, so each fold rests on a handful of geographic units and the fold-to-fold score spread widens. With too few blocks you also drop below the number of folds you want and lose the ability to average out regional effects. Aim for the smallest block that comfortably exceeds the range so you keep many blocks while blocking leakage.
Related
- Implementing SpatialKFold in Python
- Leave-One-Region-Out Cross-Validation with scikit-learn
- Evaluating Model Performance with Spatial Metrics
- Handling Spatial Autocorrelation
Part of: Spatial Cross-Validation Strategies — Training Geospatial Predictive Models in Python