spatialcv

Honest cross-validation for machine learning on spatial data.

spatialcv provides scikit-learn compatible splitters that keep spatially autocorrelated neighbours out of each test fold's training data, plus the tools to choose a block size or buffer and to check whether your cross-validation is as hard as the real prediction task.

pip install spatialcv              # numpy + scikit-learn
pip install "spatialcv[plot]"      # + matplotlib for fold maps
pip install "spatialcv[geo]"       # + geopandas (GeoDataFrame input) and matplotlib

Python 3.9+. Every splitter has split(X, y, groups) and get_n_splits(), so it plugs into cross_val_score, cross_validate, cross_val_predict, GridSearchCV, RandomizedSearchCV and pipelines unchanged.

Why spatial cross-validation?

Nearby observations are similar (Tobler's first law). Random K-fold puts a sample's neighbours in the training set, so a flexible model (random forest, gradient boosting, a neural network, anything that sees coordinates or spatially smooth covariates) can "look up" the answer instead of learning a transferable relationship. The CV score then reflects interpolation between close neighbours, while the map is used far from any sample. Large-scale mapping studies have reported random-CV accuracies that collapse under spatial validation (Ploton et al. 2020).

Spatial CV makes the test data as far from the training data as the prediction locations are. The splitters here do this in different ways:

Helpers:

Quick start

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GridSearchCV, KFold, cross_val_score

import spatialcv as scv

coords, X, y = scv.make_spatial_regression(400, random_state=42)
features = np.column_stack([X, coords])
model = RandomForestRegressor(n_estimators=100, random_state=0, n_jobs=-1)

# 1. How far does spatial autocorrelation reach?
vg = scv.estimate_autocorrelation_range(coords, y, model="auto")
print(f"range = {vg.effective_range:.0f} m, model = {vg.model}, nugget = {vg.nugget:.2f}")

# 2. Random vs spatial K-fold
random_cv = KFold(5, shuffle=True, random_state=0)
block_cv = scv.SpatialBlockKFold(5, block_size=vg.effective_range, coords=coords, random_state=0)
for name, cv in [("random KFold", random_cv), ("SpatialBlockKFold", block_cv)]:
    r2 = cross_val_score(model, features, y, cv=cv, scoring="r2")
    print(f"{name:<18} R2 = {r2.mean():.3f} +/- {r2.std():.3f}")

# 3. Hyper-parameter search with spatial folds
search = GridSearchCV(model, {"max_features": [0.3, 1.0]}, cv=block_cv, scoring="r2")
search.fit(features, y)
print("best:", search.best_params_, f"R2 = {search.best_score_:.3f}")

# 4. Are the CV test points as far from training data as the prediction grid?
grid = scv.make_prediction_grid(extent=10_000, n_per_side=40)
for name, cv in [("random KFold", random_cv), ("SpatialBlockKFold", block_cv)]:
    d = scv.fold_diagnostics(cv, features, prediction_points=grid, coords=coords)
    print(
        f"{name:<18} median test->train {np.median(d.Gjstar):6.0f} m | "
        f"prediction->train {np.median(d.Gij):6.0f} m | W = {d.W:.0f}"
    )

Output (scikit-learn 1.9, numpy 2.5):

range = 2760 m, model = spherical, nugget = 0.00
random KFold       R2 = 0.835 +/- 0.041
SpatialBlockKFold  R2 = 0.143 +/- 0.299
best: {'max_features': 0.3} R2 = 0.225
random KFold       median test->train    122 m | prediction->train    694 m | W = 711
SpatialBlockKFold  median test->train    827 m | prediction->train    694 m | W = 201

Random K-fold tests points that are typically 122 m from a training sample, while the map is predicted 694 m (median) from the nearest sample: its R² of 0.84 measures something the model will not be asked to do.

Which CV matches the true map error?

examples/compare_cv.py generates 300 clustered samples of a spatial field together with the true values on a 40 × 40 prediction grid over the whole 10 km × 10 km area, so the real map RMSE is known. Pooled CV predictions (cross_val_predict) give:

autocorrelation range of y: 2177 m (spherical)
true map RMSE: 1.035

splitter                  CV RMSE    error       W
KFold (random)              0.499   -0.537     974
SpatialBlockKFold           1.097   +0.062     482
SpatialClusterKFold         1.173   +0.138     403
SpatialKFoldWithBuffer      1.216   +0.180     552
BufferedLeaveOneOut         1.093   +0.058     760
NNDMLeaveOneOut             0.946   -0.089      71
KNNDMKFold                  1.210   +0.175     341

Random CV halves the true error. The spatial methods land within ±0.2, and NNDM has by far the best-matched distance distribution (lowest W). Spatial CV is not automatically "correct" either: strategies that hold out more than the prediction task requires are pessimistic (Wadoux et al. 2021), which is why matching the distance distributions (NNDM / kNNDM) and checking them with fold_diagnostics is worthwhile.

Coordinates

Splitters look up sample coordinates in this order:

1. coords=: an (n_samples, 2) array, or a GeoSeries / GeoDataFrame (point coordinates, or centroids for other geometries);
2. coord_columns=: two column names (DataFrame) or positions inside X;
3. X itself when it is a GeoDataFrame;
4. X itself when it has exactly two columns.

Units: with metric="euclidean" (default) coordinates are planar, e.g. metres in a projected CRS, and block_size / buffer use the same unit. With metric="haversine" coordinates are (lon, lat) in degrees and all distances (buffer, block_size, phi, variogram range, diagnostics) are great-circle metres. K-means clustering then runs on 3D unit vectors, so clusters are correct across the antimeridian; blocks are laid on a local equirectangular projection.

cv = scv.SpatialKFoldWithBuffer(5, buffer=2_000, coords=lonlat, metric="haversine")
cv = scv.SpatialClusterKFold(5, coord_columns=["x", "y"])       # columns of a DataFrame X
cv = scv.SpatialBlockKFold(5, block_size=5_000)                  # X is a GeoDataFrame

Splitter notes

Blocks. block_size should be at least the autocorrelation range. Hexagons have the same area as squares of side block_size. assignment="systematic" numbers grid cells row by row and sends cell i to fold i mod n_splits (deterministic); balance=True (random assignment) gives each block, largest first, to the fold with the fewest samples. block_index(X) returns the block id of each sample, and every partitioning splitter has assign_folds(X).

Buffers. A training point is dropped when its distance to any test point is <= buffer. Test sets still partition the data, so cross_val_predict works; training sets shrink, so very large buffers raise an error rather than yield empty training sets. Wrap any splitter: SpatialKFoldWithBuffer(buffer=1000, cv=SpatialBlockKFold(...)).

NNDM (Milà et al. 2022) needs the locations where the model will be applied (prediction_points, e.g. a regular sample of a few thousand points of the area of interest). Starting from LOO, it walks through nearest-neighbour distances r in increasing order and removes test point j's nearest neighbour from its training set whenever the LOO distance CDF Gj*(r) would still be at or above the prediction CDF Gij(r), keeping more than min_train of the data and stopping at phi (default: the largest distance). It follows CAST's reference nndm(). Cost is O(n²) memory.

kNNDM (Linnenbrink et al. 2024) clusters the samples (Ward hierarchical clustering, or k-means) into q clusters for a log-spaced grid of q between n_splits and n − 2. Clusters with at least n/k samples become folds on their own; the rest, ordered along the first principal component of the coordinates, are dealt round-robin into the remaining folds. Configurations with a fold larger than maxp are rejected and the one with the smallest Wasserstein distance W between Gj* and Gij wins. If random CV already matches (W < 1e-4) random folds are returned. It follows CAST's reference knndm() (geographical space).

res = scv.knndm(coords, grid, n_splits=5, random_state=0)
res.labels, res.n_clusters, res.W, res.method     # folds, chosen q, W, "hierarchical"
res = scv.nndm(coords, grid, phi=3_000)
res.train_indices[0], res.Gjstar, res.W

Diagnostics and plots

d = scv.fold_diagnostics(cv, X, prediction_points=grid, coords=coords)
print(d)            # medians of Gij / Gjstar / Gj, W and W_random
d.plot()            # ECDFs of the three distance distributions
scv.plot_folds(cv, X)            # colour = test fold
scv.plot_folds(cv, X, fold=0)    # train / test / excluded (buffer) points of split 0
vg.plot()                        # empirical semivariogram and fitted model

fold_diagnostics works with any object that has split, including scikit-learn's own splitters, so you can check a GroupKFold built from administrative units too.

Variogram details

estimate_autocorrelation_range(coords, values, model="spherical" | "exponential" | "auto") computes Matheron's estimator in n_lags equal-width bins up to max_lag (default: a third of the bounding-box diagonal) and fits nugget + partial sill + range by weighted least squares with weights N(h)/h² (Cressie 1985). The returned Variogram has effective_range (the range for the spherical model, 3 × the range parameter for the exponential one, i.e. 95 % of the sill), nugget, partial_sill, sill, the empirical lags/gamma/counts, and is callable: vg(h). Use model residuals as values when possible; above 2000 samples a random subset is used (max_points).

References

• Brenning, A. (2012). Spatial cross-validation and bootstrap for the assessment of prediction rules in remote sensing: the R package sperrorest. IEEE IGARSS 2012, 5372–5375. https://doi.org/10.1109/IGARSS.2012.6352393
• Cressie, N. (1985). Fitting variogram models by weighted least squares. Mathematical Geology 17, 563–586. https://doi.org/10.1007/BF01032109
• Le Rest, K., Pinaud, D., Monestiez, P., Chadoeuf, J., Bretagnolle, V. (2014). Spatial leave-one-out cross-validation for variable selection in the presence of spatial autocorrelation. Global Ecology and Biogeography 23, 811–820. https://doi.org/10.1111/geb.12161
• Linnenbrink, J., Milà, C., Ludwig, M., Meyer, H. (2024). kNNDM CV: k-fold nearest-neighbour distance matching cross-validation for map accuracy estimation. Geoscientific Model Development 17, 5897–5912. https://doi.org/10.5194/gmd-17-5897-2024
• Milà, C., Mateu, J., Pebesma, E., Meyer, H. (2022). Nearest neighbour distance matching Leave-One-Out Cross-Validation for map validation. Methods in Ecology and Evolution 13, 1304–1316. https://doi.org/10.1111/2041-210X.13851
• Ploton, P., et al. (2020). Spatial validation reveals poor predictive performance of large-scale ecological mapping models. Nature Communications 11, 4540. https://doi.org/10.1038/s41467-020-18321-y
• Pohjankukka, J., Pahikkala, T., Nevalainen, P., Heikkonen, J. (2017). Estimating the prediction performance of spatial models via spatial k-fold cross validation. International Journal of Geographical Information Science 31, 2001–2019. https://doi.org/10.1080/13658816.2017.1346255
• Rahimi, A., Recht, B. (2007). Random features for large-scale kernel machines. Advances in Neural Information Processing Systems 20 (used by the synthetic data generator).
• Roberts, D. R., et al. (2017). Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography 40, 913–929. https://doi.org/10.1111/ecog.02881
• Valavi, R., Elith, J., Lahoz-Monfort, J. J., Guillera-Arroita, G. (2019). blockCV: an R package for generating spatially or environmentally separated folds for k-fold cross-validation of species distribution models. Methods in Ecology and Evolution 10, 225–232. https://doi.org/10.1111/2041-210X.13107
• Wadoux, A. M. J.-C., Heuvelink, G. B. M., de Bruin, S., Brus, D. J. (2021). Spatial cross-validation is not the right way to evaluate map accuracy. Ecological Modelling 457, 109692. https://doi.org/10.1016/j.ecolmodel.2021.109692

The NNDM and kNNDM implementations follow the R package CAST (Meyer, Milà, Ludwig, Linnenbrink et al.); block and buffer designs follow blockCV and spatialsample.

Citation

See CITATION.cff. Please also cite the papers of the methods you use.

License

MIT © Yusuf Eminoğlu