Metadata-Version: 2.4
Name: spatialcv
Version: 0.1.0
Summary: Honest cross-validation for spatial machine learning: scikit-learn compatible spatial block, cluster, buffered, NNDM and kNNDM splitters.
Author-email: Yusuf Eminoğlu <yusufeminoglu@gmail.com>
License-Expression: MIT
Project-URL: Homepage, https://github.com/YusufEminoglu/spatialcv_sdk
Project-URL: Repository, https://github.com/YusufEminoglu/spatialcv_sdk
Project-URL: Issues, https://github.com/YusufEminoglu/spatialcv_sdk/issues
Project-URL: Documentation, https://github.com/YusufEminoglu/spatialcv_sdk#readme
Project-URL: Changelog, https://github.com/YusufEminoglu/spatialcv_sdk/blob/main/CHANGELOG.md
Keywords: spatial-cross-validation,cross-validation,spatial-autocorrelation,geoai,machine-learning,scikit-learn,block-cross-validation,nndm,knndm,variogram,gis,remote-sensing,geostatistics,model-evaluation
Classifier: Development Status :: 4 - Beta
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Operating System :: OS Independent
Classifier: Topic :: Scientific/Engineering :: GIS
Classifier: Topic :: Scientific/Engineering :: Artificial Intelligence
Classifier: Topic :: Scientific/Engineering :: Information Analysis
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Developers
Classifier: Typing :: Typed
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.21.0
Requires-Dist: scikit-learn>=1.0.0
Provides-Extra: plot
Requires-Dist: matplotlib>=3.5.0; extra == "plot"
Provides-Extra: geo
Requires-Dist: geopandas>=0.12.0; extra == "geo"
Requires-Dist: matplotlib>=3.5.0; extra == "geo"
Provides-Extra: dev
Requires-Dist: pytest>=7.0.0; extra == "dev"
Requires-Dist: pytest-cov>=4.0.0; extra == "dev"
Requires-Dist: ruff>=0.1.0; extra == "dev"
Requires-Dist: mypy>=1.0.0; extra == "dev"
Requires-Dist: pandas>=1.3.0; extra == "dev"
Dynamic: license-file

# 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.

```bash
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:

| Splitter | Idea | Test sets | Key parameters | Reference |
|---|---|---|---|---|
| `SpatialBlockKFold` | square or hexagonal blocks are held out whole; random, balanced or systematic block-to-fold assignment | partition | `block_size`, `shape`, `assignment`, `balance` | Roberts et al. 2017; Valavi et al. 2019 |
| `SpatialClusterKFold` | k-means on coordinates (on the sphere for lon/lat); clusters are held out | partition | `n_clusters`, `balance` | Brenning 2012 |
| `EnvironmentalBlockKFold` | k-means in standardised *feature* space; tests extrapolation to new conditions | partition | `features`, `standardize` | Roberts et al. 2017; Valavi et al. 2019 |
| `SpatialKFoldWithBuffer` | any K-fold (random by default) plus a dead zone: training points within `buffer` of a test point are dropped | partition | `buffer`, `cv` | Pohjankukka et al. 2017 |
| `BufferedLeaveOneOut` | leave-one-out with a circular dead zone around each test point | single samples | `buffer` | Le Rest et al. 2014 |
| `NNDMLeaveOneOut` | LOO whose training sets are thinned until the CV nearest-neighbour distance distribution matches prediction-to-sample distances | single samples | `prediction_points`, `phi`, `min_train` | Milà et al. 2022 |
| `KNNDMKFold` | k-fold version of NNDM: clusters are merged into folds to minimise the Wasserstein distance W between those distributions | partition | `prediction_points`, `n_splits`, `maxp`, `clustering` | Linnenbrink et al. 2024 |

Helpers:

| Function | Purpose |
|---|---|
| `estimate_autocorrelation_range` | empirical semivariogram + fitted spherical / exponential model (pure numpy); gives a block size or buffer |
| `fold_diagnostics` | nearest-neighbour distance distributions: prediction→train (`Gij`), CV test→train (`Gjstar`), LOO (`Gj`) and the W statistic, for *any* splitter |
| `plot_folds` | map of folds, or of train / test / excluded points of one split (matplotlib) |
| `nndm`, `knndm` | the NNDM / kNNDM algorithms as functions returning folds and distances |
| `make_spatial_regression`, `make_prediction_grid` | synthetic autocorrelated data for experiments |

## Quick start

```python
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):

```text
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`](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:

```text
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.

```python
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).

```python
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

```python
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`](CITATION.cff). Please also cite the papers of the methods you use.

## License

MIT © Yusuf Eminoğlu
