Spatial Block Cross-Validation for SDM in Python

Random k-fold cross-validation reports a species distribution model as far better than it is, because spatially adjacent occurrence records leak across the train/test partition and the model is scored on ground it has effectively already seen. This guide isolates the fix: assigning occurrences to spatial blocks — regular grids or a checkerboard — and driving the split through scikit-learn’s GroupKFold so no block ever spans both partitions. It is a focused technique within the parent Model Validation & AUC Metrics workflow, part of the broader Species Distribution Modeling with MaxEnt pipeline. The result is an AUC that reflects transfer to new geography rather than memory of the sampled patches.

When to use spatial blocks and which layout

Any presence-only model built from clustered occurrences needs spatially structured folds; the only question is the block layout. Regular systematic blocks are the general default; a checkerboard interleaves training and test regions to balance environmental coverage when a single contiguous hold-out would leave a fold with no examples of some conditions.

Choose the block edge larger than the range of residual spatial autocorrelation — a variogram on a pilot fit gives it, and 10–50 km is typical for regional forest taxa. Pair spatial folds with spatial thinning of occurrence records first: thinning removes duplication within survey patches, blocking removes leakage between train and test.

Minimal reproducible example

The block-id assignment is pure arithmetic on projected coordinates. Below, points_gdf is a GeoDataFrame of presence and background points in a projected CRS (metres), sampled against the environmental predictor stack; X is the feature matrix and y is the label (1 presence, 0 background).

import numpy as np
import geopandas as gpd
from sklearn.model_selection import GroupKFold
from sklearn.metrics import roc_auc_score


def assign_block_ids(points_gdf, block_size_m=25000):
    """Bucket each point into a square grid cell and return integer ids."""
    xs = points_gdf.geometry.x.to_numpy()
    ys = points_gdf.geometry.y.to_numpy()
    col = np.floor((xs - xs.min()) / block_size_m).astype(int)
    row = np.floor((ys - ys.min()) / block_size_m).astype(int)
    n_cols = col.max() + 1
    return row * n_cols + col, row, col


def checkerboard_ids(row, col):
    """Two interleaved folds: cell colour is (row + col) parity."""
    return (row + col) % 2


def spatial_block_cv(model_factory, X, y, groups):
    """GroupKFold on block ids so no block spans train and test."""
    n_folds = min(5, len(np.unique(groups)))
    gkf = GroupKFold(n_splits=n_folds)
    per_fold_auc = []
    for train_idx, test_idx in gkf.split(X, y, groups=groups):
        model = model_factory()
        model.fit(X[train_idx], y[train_idx])
        scores = model.predict_proba(X[test_idx])[:, 1]
        per_fold_auc.append(roc_auc_score(y[test_idx], scores))
    return np.array(per_fold_auc)


block_id, row, col = assign_block_ids(points_gdf, block_size_m=25000)
fold_auc = spatial_block_cv(model_factory, X, y, groups=block_id)
print(f"Spatial block AUC: {fold_auc.mean():.3f} "
      f"(per fold {fold_auc.round(3).tolist()})")

# Checkerboard alternative (two interleaved folds)
cb_id = checkerboard_ids(row, col)
cb_auc = spatial_block_cv(model_factory, X, y, groups=cb_id)
print(f"Checkerboard AUC: {cb_auc.mean():.3f}")

GroupKFold guarantees every block lands entirely in one fold, which is exactly the property random KFold lacks. model_factory is any zero-argument callable returning an unfitted estimator with fit and predict_proba — an elapid MaxentModel wrapper, a sklearn classifier, or a maxnet fit.

Random versus spatial: the inflation you are removing

Run both partitions on the same data to see the leakage the spatial folds eliminate. The random split reuses GroupKFold machinery only conceptually — a plain KFold shuffles points regardless of geography.

from sklearn.model_selection import KFold


def random_cv(model_factory, X, y, k=5, seed=0):
    kf = KFold(n_splits=k, shuffle=True, random_state=seed)
    aucs = []
    for tr, te in kf.split(X):
        m = model_factory()
        m.fit(X[tr], y[tr])
        aucs.append(roc_auc_score(y[te], m.predict_proba(X[te])[:, 1]))
    return np.array(aucs)


rand_auc = random_cv(model_factory, X, y)
print(f"Random CV AUC:  {rand_auc.mean():.3f}")
print(f"Spatial CV AUC: {fold_auc.mean():.3f}")
print(f"Optimism (inflation): {rand_auc.mean() - fold_auc.mean():+.3f}")

A positive optimism value is the inflation you would otherwise have published. A wide spread across spatial folds is not a defect of the method but a finding: it localizes where a driver is missing from the predictor stack.

Parameter reference

Expected output and verification

Two invariants distinguish a correct spatial split from a broken one: no block may appear in both a train and a test partition, and the spatial AUC should be no higher than the random AUC. If spatial CV scores higher than random, the blocks are too small to break autocorrelation.

import numpy as np
from sklearn.model_selection import GroupKFold


def verify_spatial_cv(X, y, groups, rand_mean, spatial_mean):
    gkf = GroupKFold(n_splits=min(5, len(np.unique(groups))))
    for train_idx, test_idx in gkf.split(X, y, groups=groups):
        shared = set(groups[train_idx]) & set(groups[test_idx])
        assert not shared, f"block(s) {shared} leaked across the split"

    # Spatial CV should not exceed random CV; if it does, blocks too small.
    assert spatial_mean <= rand_mean + 1e-6, (
        "spatial AUC exceeds random AUC — increase block_size_m"
    )
    print(f"OK — no block leakage; optimism removed "
          f"{rand_mean - spatial_mean:+.3f} AUC")

Once the split is verified, feed the pooled out-of-fold scores into the discrimination and calibration metrics the parent workflow defines. Pair AUC with the continuous Boyce index, which is designed for presence-only settings where a true-absence AUC is awkward. Log the block size, layout, seed, and per-fold AUCs beside the suitability layer so later model updates stay comparable.

Common pitfalls

• Block size below the autocorrelation range. Small blocks still place near-duplicate neighbours in adjacent folds, so the AUC stays inflated. Size the block from a variogram of pilot-fit residuals, not by guesswork.
• Running the split in geographic degrees. Grid arithmetic on longitude/latitude produces cells that are physically unequal and distances that are wrong. Reproject to a metre-based CRS before assigning blocks.
• Using plain KFold and calling it spatial. Only GroupKFold (or an explicit block-to-fold map) keeps a block intact. A shuffled KFold on block-tagged data still leaks — the tag has to be the grouping key.
• Reporting the mean and hiding the spread. A mean AUC of 0.85 over folds of {0.95, 0.65} is not an 0.85 model. Report the per-fold range and map it geographically to find where the model collapses.

Frequently asked questions

Related

• Model Validation & AUC Metrics — the parent workflow this folding technique feeds
• Evaluating SDM with the Continuous Boyce Index — a presence-only calibration metric to pair with spatial-fold AUC
• Spatial Thinning of Occurrence Records in Python — remove within-patch duplication before blocking between folds
• MaxEnt Model Training & Tuning — the model whose hyperparameters these folds help select
• Species Distribution Modeling with MaxEnt — the full presence-only modeling pipeline