Buffered Leave-One-Out Cross-Validation

TL;DR: For each observation, remove the point itself and every training point within the variogram range using a scipy.spatial.cKDTree radius query, fit or interpolate on what remains, predict the held-out value, and collect the residual. Aggregating those residuals as sqrt(mean(res**2)) gives a buffered leave-one-out (BLOO) RMSE that is not inflated by autocorrelated neighbours sitting a few metres from each test point.

Why This Matters #

Ordinary leave-one-out cross-validation removes a single point and predicts it from everything else — including its immediate neighbours. When data are spatially autocorrelated, those neighbours are near-duplicates of the held-out value, so the score is optimistic in exactly the way random K-fold is. Buffered LOO closes the loophole by excising a disc of radius equal to the variogram range around each test point, so the model must predict across a gap it genuinely has not observed. This page extends the ideas in the Cross-Validation Strategies guide and belongs to the wider Python Workflows for Spatial Modeling & Regression stack. The buffer radius comes straight from a fitted variogram, so this technique is tightly coupled to the Variogram Modeling & Semivariance Analysis pillar.

BLOO is the most conservative common validation for point interpolation. Where spatial block CV gives a handful of regional estimates, BLOO produces a residual at every location, letting you map where the model is weak, not just report a single number.

Why the Buffer Radius Is the Range #

The variogram $\gamma(h)$ rises from the nugget toward the sill as separation $h$ grows. The range $a$ is the distance at which

$$\gamma(a) \approx \text{sill}, \qquad \text{Corr}(Z(s), Z(s+h)) \approx 0 \ \text{for}\ h \ge a.$$

Two observations closer than $a$ are correlated, so leaving a neighbour inside the buffer leaks information about the held-out value. Excluding everything within $a$ makes the nearest surviving training point uncorrelated with the target, which is the condition for an honest test.

Formally, the buffered residual at location $s_i$ is

$$r_i = Z(s_i) - \hat{Z}_{-B(s_i,a)}(s_i),$$

where $\hat{Z}_{-B(s_i,a)}$ is the model fitted with every point in the disc $B(s_i, a)$ of radius $a$ removed. The aggregate error is $\text{RMSE} = \sqrt{\tfrac{1}{n}\sum_i r_i^2}$. Because the training set excludes the correlated neighbourhood, $r_i$ measures prediction skill at separation $a$ rather than skill at interpolating between near-duplicates.

Environment #

pip install \
  geopandas==0.14.4 \
  scikit-gstat==1.0.18 \
  scipy==1.13.0 \
  numpy==1.26.4 \
  scikit-learn==1.4.2

import numpy as np
import geopandas as gpd
from scipy.spatial import cKDTree
from skgstat import Variogram
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error

Assume a projected GeoDataFrame gdf with target y and predictors feat_cols. Coordinates must be metric so the buffer radius is in metres.

Step-by-Step Implementation #

Step 1 — Estimate the Buffer Radius from a Variogram #

Fit an empirical variogram to the target and read off the range. That range becomes the exclusion radius.

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

vgram = Variogram(coords, y, model="spherical", n_lags=15, normalize=False)
buffer_radius = float(vgram.parameters[0])   # [range, sill, nugget] for spherical
print(f"Variogram range (buffer radius): {buffer_radius:,.0f} m")

If the variogram is noisy, round the range up rather than down; a slightly generous buffer errs toward honesty.

Step 2 — Build a Spatial Index #

A cKDTree answers “which points lie within radius r of this location” in logarithmic time, which keeps the per-point buffer step cheap.

tree = cKDTree(coords)

Step 3 — Exclude the Buffer and Predict Each Point #

Loop over observations. For each, query the tree for all indices inside buffer_radius, drop them from the training set, fit on the remainder, and predict the single held-out point.

def bloo_cv(coords, X, y, tree, buffer_radius, seed=0):
    """Buffered leave-one-out CV. Returns per-point residuals and RMSE."""
    n = len(y)
    residuals = np.full(n, np.nan, dtype="float64")

    for i in range(n):
        # Indices to exclude: the point itself + everything within the buffer
        excluded = set(tree.query_ball_point(coords[i], r=buffer_radius))
        excluded.add(i)
        train_mask = np.ones(n, dtype=bool)
        train_mask[list(excluded)] = False

        if train_mask.sum() < 10:
            continue  # too few training points to fit reliably

        model = RandomForestRegressor(
            n_estimators=200, random_state=seed, n_jobs=-1,
        )
        model.fit(X[train_mask], y[train_mask])
        pred = model.predict(X[i].reshape(1, -1))[0]
        residuals[i] = y[i] - pred

    valid = ~np.isnan(residuals)
    rmse = np.sqrt(np.mean(residuals[valid] ** 2))
    return residuals, rmse

X = gdf[feat_cols].to_numpy()
residuals, bloo_rmse = bloo_cv(coords, X, y, tree, buffer_radius)
print(f"Buffered LOO RMSE: {bloo_rmse:.3f}  (n={np.sum(~np.isnan(residuals))})")

Step 4 — Compare Against Naive LOO #

Run ordinary leave-one-out (buffer radius zero, only the point itself removed) so you can see the optimism the buffer removes.

def naive_loo_cv(X, y, seed=0):
    n = len(y)
    residuals = np.empty(n)
    for i in range(n):
        mask = np.ones(n, dtype=bool)
        mask[i] = False
        model = RandomForestRegressor(n_estimators=200, random_state=seed, n_jobs=-1)
        model.fit(X[mask], y[mask])
        residuals[i] = y[i] - model.predict(X[i].reshape(1, -1))[0]
    return np.sqrt(np.mean(residuals ** 2))

naive_rmse = naive_loo_cv(X, y)
print(f"Naive LOO RMSE:    {naive_rmse:.3f}")
print(f"Buffered LOO RMSE: {bloo_rmse:.3f}")
print(f"Optimism ratio:    {bloo_rmse / naive_rmse:.2f}x")

Step 5 — Map the Residuals #

Because BLOO yields a residual per location, attach them back to the GeoDataFrame and inspect their spatial pattern.

gdf = gdf.assign(bloo_residual=residuals)
gdf["abs_residual"] = gdf["bloo_residual"].abs()
print(gdf[["abs_residual"]].describe())
# gdf.plot(column="abs_residual", legend=True)  # reveals where the model is weak

Interpreting the Output #

The BLOO RMSE is your headline accuracy for predicting at an unsampled location. Expect it to exceed the naive LOO RMSE; the ratio quantifies how much ordinary LOO flatters the model. A ratio near 1.0 implies weak short-range autocorrelation (the buffer removed little useful signal), while a ratio of 1.5x–2.5x is typical for densely sampled, strongly autocorrelated fields. If the mapped abs_residual clusters in particular regions, the model is locally biased there — often a sign of an unmodelled trend or of anisotropy that a single isotropic range fails to capture.

Watch the valid count n. Points near the study-area edge lose most of their neighbourhood to the buffer and may be skipped by the < 10 guard; a large skipped fraction means the buffer is large relative to sampling density and the estimate rests on interior points only.

It also helps to summarise the residual distribution, not just its RMSE. The mean residual should sit near zero — a systematic offset signals bias, typically an unmodelled large-scale trend that the local training points cannot recover once the buffer removes the neighbourhood. The interquartile range and the presence of a heavy tail tell you whether a few hard locations dominate the error or whether the model is uniformly mediocre. Reporting both the RMSE and a mapped residual surface turns BLOO from a single accuracy score into a diagnostic of where and why the model fails, which is far more actionable when deciding whether to add covariates, densify sampling, or switch interpolator.

Critical Best Practices #

Take the Radius from a Directional Variogram Under Anisotropy #

A single isotropic range under-buffers along the direction of strongest correlation. If a directional variogram shows geometric anisotropy, buffer with an ellipse using the major-axis range, or take the maximum directional range as a conservative isotropic radius.

Guard Against Empty Neighbourhoods #

Edge points and sparsely sampled regions can lose almost all training data to the buffer. Enforce a minimum training count (the < 10 check above) and record which points were skipped, so the reported RMSE is not silently based on the easy interior only.

Match the Estimator to Deployment #

BLOO validates whatever model you actually deploy. If production uses kriging, refit kriging inside the loop, not a random forest. The buffer logic is estimator-agnostic; only the fit/predict lines change.

Budget the Compute Cost #

BLOO performs n model fits, which is expensive for large n. For big datasets, either subsample the held-out points to a representative few thousand and refit only for those, or fall back to the cheaper Spatial Block Cross-Validation in Python approach. Keep the cKDTree query — it is not the bottleneck; the refits are. The loop is embarrassingly parallel, so joblib.Parallel across held-out indices scales almost linearly with cores; just build the tree once and share it read-only across workers.

Do Not Let the Buffer Erase the Trend Signal #

A large buffer removes not only correlated noise but also the local data a purely spatial interpolator needs to anchor its prediction. If the target has a strong deterministic trend, keep trend covariates in the feature matrix (elevation, distance-to-coast, land cover) so the model still has information to predict from after the neighbourhood is excised. This is the same rationale behind regression-kriging: let covariates carry the large-scale structure while the spatial term handles the residual. Without covariates, an excessively large buffer can make interior predictions collapse toward the global mean and inflate RMSE for reasons unrelated to true skill.

Troubleshooting #

Next Steps #

BLOO and blocked CV answer complementary questions; pair this with Spatial Block Cross-Validation in Python for a faster regional estimate, and review Spatial K-Fold Cross-Validation Setup for the fold-construction fundamentals both rely on.

Related:

• Spatial Block Cross-Validation in Python — a faster, region-level alternative to per-point buffering
• Spatial K-Fold Cross-Validation Setup — foundational spatial fold mechanics
• Variogram Modeling & Semivariance Analysis — where the buffer radius (the range) comes from

← Back to Cross-Validation Strategies