Geostatistical Modeling Spatial Statistics with Python

Spatial K-Fold Cross-Validation Setup in Python

TL;DR: Call KMeans(n_clusters=n_splits).fit_predict(coords) on your projected X/Y coordinates to create spatial groups, then pass those groups to sklearn.model_selection.GroupKFold. That single pattern replaces random KFold and eliminates the geographic information leakage that causes inflated R² scores on spatial datasets.

Why This Matters

Standard random K-fold violates the independence assumption the moment points are spatially autocorrelated — and almost every environmental, urban, or earth-science dataset is. When a random split scatters nearby points across train and test folds, the model exploits local spatial patterns it should never have seen, routinely inflating R² by 15–40% and underreporting RMSE by a comparable margin.

Spatial K-fold is the practical fix. It groups geographically proximate points together, then assigns entire groups to either the training or the test partition, maintaining a spatial gap between the two. This page shows the full setup in Python, including fold-distance validation so you can confirm separation before trusting a single metric.

This task sits within the cross-validation strategies workflow, which in turn is part of the broader Python Workflows for Spatial Modeling & Regression pipeline.

Environment and Version Pinning

text
pip install scikit-learn>=1.4 geopandas>=0.14 numpy>=1.26 scipy>=1.12
python
import numpy as np                                     # 1.26+
import geopandas as gpd                                # 0.14+
from sklearn.cluster import KMeans                     # scikit-learn 1.4+
from sklearn.model_selection import GroupKFold         # scikit-learn 1.4+
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_absolute_error, r2_score
from scipy.spatial.distance import cdist              # scipy 1.12+

All coordinate operations require a projected metric CRS (e.g., UTM or a local equal-area projection). Do not pass geographic lat/lon — Euclidean distances in degree units are geometrically meaningless for clustering. Reproject with GeoDataFrame.to_crs() before extracting coordinates.

Step-by-Step Implementation

Step 1 — Extract projected coordinates

python
def extract_coords(data):
    """
    Return an (n, 2) float array of projected X/Y coordinates.

    Parameters
    ----------
    data : np.ndarray or gpd.GeoDataFrame
        Either a shape-(n, 2) array of [X, Y] columns, or a GeoDataFrame
        whose geometry is already in a metric CRS.
    """
    if isinstance(data, gpd.GeoDataFrame):
        return np.column_stack([data.geometry.x.values,
                                data.geometry.y.values])
    coords = np.asarray(data)
    if coords.ndim != 2 or coords.shape[1] != 2:
        raise ValueError("Expected shape (n, 2) for coordinate array.")
    return coords

Step 2 — Cluster coordinates into spatial groups

KMeans assigns every point to one of n_splits compact geographic clusters. Every resulting group becomes a fold, so points within the same spatial group always stay on the same side of the train/test boundary.

python
def make_spatial_groups(coords, n_splits=5, random_state=42):
    """
    Cluster projected coordinates into n_splits spatially contiguous groups.

    Returns
    -------
    groups : np.ndarray, shape (n,)
        Integer cluster label for each point (0 … n_splits-1).
    """
    if len(coords) < n_splits:
        raise ValueError(
            f"Need at least {n_splits} samples for {n_splits} folds, "
            f"got {len(coords)}."
        )
    km = KMeans(n_clusters=n_splits, random_state=random_state, n_init="auto")
    return km.fit_predict(coords)

Step 3 — Generate fold indices with GroupKFold

GroupKFold guarantees that no group appears in both train and test within a single fold. Combined with spatial groups, this enforces geographic separation automatically.

python
def spatial_kfold_split(coords_or_gdf, n_splits=5, random_state=42):
    """
    Build spatially separated K-fold indices.

    Parameters
    ----------
    coords_or_gdf : np.ndarray or gpd.GeoDataFrame
        Projected coordinates, shape (n, 2), or a GeoDataFrame.
    n_splits : int
        Number of spatial folds.
    random_state : int
        Seed for KMeans reproducibility.

    Returns
    -------
    folds : list of (train_idx, test_idx) tuples
    groups : np.ndarray
        Cluster labels — store these alongside your data for
        reproducible reporting.
    """
    coords = extract_coords(coords_or_gdf)
    groups = make_spatial_groups(coords, n_splits=n_splits,
                                  random_state=random_state)
    gkf = GroupKFold(n_splits=n_splits)
    dummy_y = np.zeros(len(coords))
    folds = list(gkf.split(X=coords, y=dummy_y, groups=groups))
    return folds, groups

Step 4 — Validate inter-fold distances

Generating folds is not enough — you need to confirm that the spatial gap between training and test points is at least as wide as the autocorrelation range of your target variable. Estimate that range from a variogram before choosing a threshold.

python
def validate_fold_separation(coords, folds, min_distance=1000.0):
    """
    Print a pass/fail summary for each fold's minimum train-test distance.

    Parameters
    ----------
    coords : np.ndarray, shape (n, 2)
        Projected coordinates (same units as min_distance).
    folds : list of (train_idx, test_idx)
    min_distance : float
        Minimum acceptable separation in CRS units (metres for UTM).
    """
    all_ok = True
    for i, (train_idx, test_idx) in enumerate(folds):
        dm = cdist(coords[test_idx], coords[train_idx], metric="euclidean")
        min_d = dm.min()
        status = "OK" if min_d >= min_distance else "WARNING"
        if status == "WARNING":
            all_ok = False
        print(f"  Fold {i+1}: min train-test distance = {min_d:,.1f}  [{status}]")
    if not all_ok:
        print("\n  At least one fold is below the distance threshold. "
              "Increase n_splits or switch to grid-based blocking.")

Step 5 — Fit, score, and collect per-fold metrics

python
if __name__ == "__main__":
    # --- Synthetic demo: 600 points in a 100 km × 100 km box (metres) ---
    rng = np.random.default_rng(42)
    n = 600
    coords = rng.uniform(0, 100_000, (n, 2))       # metres
    # Target: a gentle spatial trend plus noise
    y = 0.0003 * coords[:, 0] + 0.0005 * coords[:, 1] + rng.normal(0, 5, n)

    folds, groups = spatial_kfold_split(coords, n_splits=5)

    print("Fold separation audit (threshold = 1 000 m):")
    validate_fold_separation(coords, folds, min_distance=1_000.0)
    print()

    results = []
    for i, (train_idx, test_idx) in enumerate(folds):
        model = RandomForestRegressor(n_estimators=200, random_state=42)
        model.fit(coords[train_idx], y[train_idx])
        preds = model.predict(coords[test_idx])
        mae = mean_absolute_error(y[test_idx], preds)
        r2  = r2_score(y[test_idx], preds)
        results.append({"fold": i + 1, "n_test": len(test_idx),
                        "MAE": mae, "R2": r2})
        print(f"  Fold {i+1} | n_test={len(test_idx):4d} | "
              f"MAE={mae:.3f} | R²={r2:.3f}")

    maes = [r["MAE"] for r in results]
    r2s  = [r["R2"]  for r in results]
    print(f"\n  Mean MAE = {np.mean(maes):.3f} ± {np.std(maes):.3f}")
    print(f"  Mean R²  = {np.mean(r2s):.3f} ± {np.std(r2s):.3f}")

Interpreting the Output

Per-fold MAE and R² are your primary diagnostic. Look for two things:

  • High variance across folds (e.g., R² ranging from 0.4 to 0.9) signals that cluster sizes are uneven or that the spatial process has strong local anisotropy. Try increasing n_splits or switching to a grid-based blocking strategy.
  • R² near 1.0 on every fold when the underlying process should not be that predictable is a residual-leakage warning. Re-run validate_fold_separation with a tighter threshold, or compare against a random-KFold baseline — if spatial CV gives dramatically lower R², the random CV was leaking.

The groups array returned alongside folds is worth storing as a column in your GeoDataFrame. It enables spatial error mapping (e.g., plotting residuals coloured by fold assignment) and makes the validation fully reproducible.

Fold size imbalance is normal with KMeans because clusters adapt to point density. Imbalance above roughly 3:1 between the largest and smallest test set is worth investigating — large test folds deflate per-fold variance estimates, while tiny test folds inflate it.

Critical Best Practices

Match blocking scale to the autocorrelation range

The spatial scale of your folds must exceed the autocorrelation range of the target variable. Use spatial autocorrelation metrics or fit a variogram to estimate the range before choosing n_splits. If the range is 5 km and your folds are only 2 km wide, leakage persists despite spatial grouping.

Always use a projected CRS before clustering

KMeans computes Euclidean distances. Feeding it WGS84 latitude/longitude produces clusters distorted by the non-linear relationship between degrees and metres — a 1° longitude near the equator is ~111 km but only ~55 km at 60° N. Reproject to a locally appropriate metric CRS (e.g., EPSG:32633 for UTM zone 33N) before calling extract_coords.

Handle cluster imbalance proactively

If the largest fold contains more than three times as many test points as the smallest, consider grid-based blocking instead. Construct a regular grid over the bounding box, assign each point to the grid cell it falls in, then use cell IDs as GroupKFold groups. This guarantees roughly uniform geographic coverage at the cost of potentially unequal sample counts per cell.

Preserve temporal ordering in spatiotemporal data

For datasets with both spatial and temporal dimensions, apply spatial blocking first and then sort folds chronologically within each spatial block. Mixing temporal and spatial dimensions in a single clustering step confounds two orthogonal sources of leakage and makes the evaluation uninterpretable.

Store and version fold assignments

Assign the groups array as a column (gdf["spatial_fold"] = groups) and commit it to your data version-control system. Fold geometry changes between runs if random_state or n_splits changes, silently invalidating comparisons. Pinning the fold assignment ensures that model A and model B are always evaluated on exactly the same spatial partitions.

Troubleshooting

Symptom Likely cause Fix
ValueError: n_samples must be >= n_splits Fewer data points than folds Reduce n_splits or increase the dataset size
Min train-test distance = 0.0 Points at identical coordinates (duplicates) Deduplicate or apply a small spatial jitter before clustering
R² identical between random and spatial CV Autocorrelation range shorter than cluster diameter Use a finer variogram to re-estimate range; reduce n_splits
KMeans fails to converge Coordinates not normalised, very large CRS units Standardise coordinates with StandardScaler before clustering, then discard the scaler
Highly unequal fold sizes (> 5:1 ratio) Point clustering creates one dominant dense group Switch to grid-based blocking; choose cell size ≈ autocorrelation range
GroupKFold returns fewer splits than expected One or more clusters has zero points Verify that n_clusters == n_splits and no cluster is empty after fitting

Next Steps

For the broader context — spatial block CV, leave-one-out variants, and how to choose between them — see the cross-validation strategies guide. If your model outputs include kriging variance surfaces rather than point predictions, the uncertainty and variance mapping workflow explains how to interpret and validate those outputs alongside fold-based metrics.

Spatial K-Fold partitioning workflow A left-to-right flow showing: projected coordinates fed into KMeans clustering, which produces group labels, which are passed to GroupKFold to yield per-fold train and test index arrays. Projected coordinates (n × 2, metres) KMeans clustering n_clusters = k Group labels [0, 2, 1, 0 …] shape (n,) GroupKFold .split(X, y, groups) → train/test idx Fold 1: groups {1,2,3,4} train | group {0} test Fold 2: groups {0,2,3,4} train | group {1} test Fold 3 … k follows same pattern Spatial K-Fold: coordinate clustering → grouped splits Must be metric CRS (e.g. UTM / EPSG:326xx)

← Back to Cross-Validation Strategies for Spatial Modeling in Python