4. Spatial leakage in spatial simulation

Spatial leakage is the failure mode that most often survives a clean-looking validation report. When observations are spatially dependent, a random train/test split does not produce an independent test set — it produces a test set whose members are near-copies of training samples. The model then reports high accuracy on data it has effectively already seen, and the error surfaces only after deployment, when predictions are made at genuinely new locations. This chapter treats three coupled mechanisms: boundary effects that contaminate features near the edge of a study area, spatial autocorrelation that makes nearby samples redundant, and cross-validation designs that must be built around both.

Boundary effects: contamination at the edge of the domain

A boundary effect arises whenever a feature computed at location s depends on data outside the study area, or on data whose availability changes near the edge. Three common variants matter in simulation work:

• Window truncation. A moving-window statistic (local mean, focal variance, kernel density) computed near the edge uses fewer neighbours than the same statistic in the interior. The feature distribution is therefore non-stationary in a way that is purely an artefact of the domain boundary.
• Buffer leakage. If the study area is a crop from a larger raster, pixels just outside the crop still influence interpolated surfaces and neighbourhood features. Cropping without a buffer silently imports information from outside the declared domain.
• Edge-adjacent splitting. If train and test polygons share a boundary, features on one side are computed from data on the other. The split is spatially leaky even though the polygons are disjoint.

The practical remedy is to define an explicit analysis domain and an explicit buffer. Compute all neighbourhood features on the buffered extent, then discard the buffer before splitting. For a kernel with bandwidth h, a buffer of at least h (and preferably 2h for compactly supported kernels) is a defensible default. Record the buffer width alongside the model so the decision is auditable.

Spatial autocorrelation: why random splits fail

Spatial autocorrelation is the tendency for observations close in space to be more similar than observations far apart. It is usually quantified with Moran's I:

I = (n / W) * (Σ_i Σ_j w_ij (x_i - x̄)(x_j - x̄)) / Σ_i (x_i - x̄)²

where w_ij is a spatial weights matrix (row-standardised contiguity or inverse-distance weights are typical) and W is the sum of all weights. Values near +1 indicate strong positive autocorrelation, near 0 indicate randomness, and negative values indicate dispersion. For a simulated field with a Gaussian variogram of range r, Moran's I at lag distance d typically remains above 0.3 out to roughly d ≈ r/2, which is the scale at which random splitting starts to leak.

The consequence for validation is direct: if the autocorrelation range is 40 km and a random 20% test split is drawn, most test points will have a training neighbour within 10 km. The test set measures interpolation, not extrapolation. The reported RMSE is optimistic by an amount that grows with the autocorrelation range and shrinks only when the test points are pushed beyond that range.

Cross-validation designs that respect spatial structure

Three designs cover most simulation workflows, and they answer different questions:

Block size is the critical parameter. A block smaller than the autocorrelation range leaves training neighbours inside the buffer around each test block, so leakage persists. A practical rule is to set block edge length to at least the autocorrelation range r, and to verify by computing, for every test point, the distance to the nearest training point. If the 10th percentile of that distance is below r/2, increase the block size or add an explicit exclusion buffer between folds.

The following figure shows the three designs on the same simulated field, with the autocorrelation range marked and the leakage path made explicit.

Diagnosing leakage before trusting a score

Two diagnostics should accompany every spatial validation run. First, compute the nearest-training-neighbour distance for each test point and plot its distribution against the autocorrelation range. A test set whose median nearest-neighbour distance is below r/2 is not an independent test set. Second, compute Moran's I on the residuals of the held-out predictions. If residuals are strongly autocorrelated, the model has captured local structure but not the spatial process, and the validation score reflects interpolation rather than generalisation.

A reproducible check in Python, using a simulated field with known range, looks like this:

import numpy as np
from scipy.spatial import cKDTree

# coords: (n, 2) array in metres; y: simulated target
# r: autocorrelation range from the variogram fit (metres)
train_idx, test_idx = split_indices  # from chosen CV design
tree = cKDTree(coords[train_idx])
d_nn, _ = tree.query(coords[test_idx], k=1)

print("median NN distance:", np.median(d_nn))
print("10th percentile:", np.percentile(d_nn, 10))
print("leakage flag:", np.median(d_nn) < r / 2)

If the leakage flag is true, the design must be revised before any score is reported. Increasing block size, adding an exclusion buffer, or switching to leave-one-region-out are the three levers; the choice depends on whether the deployment question is interpolation within the domain or transfer beyond it.

Putting the pieces together

Boundary effects, autocorrelation, and cross-validation design are not independent problems. A buffer that removes boundary contamination also changes the effective sample size available to each fold; a block size chosen to exceed the autocorrelation range reduces the number of independent test blocks and widens confidence intervals. The engineering judgement is to make these trade-offs explicit rather than implicit. In practice, a defensible workflow is: fit the variogram and record r; define the analysis domain with a buffer of at least r; choose block CV with block edge ≥ r for interpolation claims, or leave-one-region-out for transfer claims; run the nearest-neighbour diagnostic; and only then report the score, with all parameters attached. When the diagnostic fails, the correct response is to redesign the split, not to reinterpret the number.