Leave-One-Region-Out Cross-Validation with scikit-learn

Use scikit-learn LeaveOneGroupOut for leave-one-region-out cross-validation with administrative or ecological region groups to test true spatial transferability.

TL;DR Answer

When your data has natural spatial groupings — watersheds, ecoregions, provinces — assign each sample a region id and pass it as the groups argument to sklearn.model_selection.LeaveOneGroupOut. Each fold holds out one entire region as the test set and trains on all the others, so every score answers the question that actually matters in deployment: how well does the model predict in a place it has never seen? This is one of the core Spatial Cross-Validation Strategies for producing honest, transferable performance estimates.

from sklearn.model_selection import LeaveOneGroupOut, cross_val_score
from sklearn.ensemble import RandomForestRegressor

logo = LeaveOneGroupOut()
rf = RandomForestRegressor(n_estimators=300, random_state=42, n_jobs=-1)

scores = cross_val_score(rf, X, y, groups=region_id, cv=logo,
                         scoring="r2")
print(f"Regions held out: {logo.get_n_splits(groups=region_id)}")
print(f"Mean transfer R²: {scores.mean():.3f}")

Part of: Spatial Cross-Validation Strategies

Why Region-Held-Out Splitting Simulates Deployment

Most geospatial models are not deployed to predict on scattered points inside their training footprint. They are deployed to a new region — a watershed that was never surveyed, an administrative unit that just requested a map, an ecoregion added to the mandate this year. A random k-fold split, and even a coordinate-clustered split, cannot answer whether the model survives that jump, because they still let the model peek at data from every region during training.

Leave-one-region-out (LORO) cross-validation makes the evaluation match the deployment. By holding out one complete region at a time, the test fold is separated from training data not just by distance but by whatever unmodeled regional structure exists: a different soil parent material, a different rainfall regime, a different sensor calibration, a different mapping crew. If the model only learned associations that hold within a region, its LORO score collapses — and that collapse is exactly the signal you want before you ship, not after.

This complements the leakage controls covered in Reducing Spatial Leakage in Model Training. Leakage reduction fixes the features and the split boundary; region-held-out validation stress-tests the result against the hardest realistic generalization gap. The two techniques are most useful together: reduce leakage during feature engineering, then confirm transferability with LORO.

The trade-off is that LORO answers a harder question than SpatialKFold or block cross-validation. Expect lower, more variable scores. That is not a flaw in the method — it is the method correctly reporting a risk that optimistic splits hide.

Core Principles for Region-Based Grouping

  • Groups come from geography, not from the model. Assign region labels with a spatial join against an authoritative polygon layer. Never derive groups from clustering the features, or you leak target structure into the split definition.
  • Every region is held out exactly once. LeaveOneGroupOut produces one fold per unique group value. There is no n_splits to tune — the number of folds equals the number of regions.
  • No group may straddle the boundary. A correct implementation guarantees a region’s samples are either all in train or all in test for a given fold. Verify this rather than assuming it.
  • Match the CRS before the join. Points and region polygons must share a coordinate system, or sjoin silently returns no matches. Reproject both to a common CRS first.
  • Handle unmatched points explicitly. Points falling outside every polygon produce a null group. Drop them or snap them to the nearest region deliberately — never let them collapse into an accidental NaN fold.
  • Weight aggregation by region size. Regions vary wildly in sample count. A plain mean of fold scores over-weights tiny regions. Report per-region scores and a size-weighted aggregate.

Assigning Group Labels from a Spatial Join

The groups array is the whole game. Build it by joining your sample points to a region polygon layer and lifting the region identifier onto each point.

import geopandas as gpd
import numpy as np

# Sample points with features and target
samples = gpd.read_file("soil_samples.geojson").to_crs("EPSG:5070")

# Authoritative region polygons (e.g. Level-III ecoregions)
regions = gpd.read_file("ecoregions.gpkg").to_crs("EPSG:5070")

# Spatial join: each point inherits the region it falls within
joined = gpd.sjoin(
    samples,
    regions[["region_id", "geometry"]],
    how="left",
    predicate="within",
)

# Fail loudly on points outside every polygon rather than shipping NaN groups
unmatched = joined["region_id"].isna().sum()
if unmatched:
    raise ValueError(
        f"{unmatched} points fell outside all region polygons. "
        "Inspect coverage or snap to nearest region before grouping."
    )

feature_cols = ["ndvi", "elevation", "slope", "aspect", "clay_pct"]
X = joined[feature_cols].to_numpy()
y = joined["organic_carbon"].to_numpy()
region_id = joined["region_id"].to_numpy()

print(f"Samples: {len(X)}  |  Regions: {np.unique(region_id).size}")

Two failure modes dominate here. The first is a CRS mismatch — if samples and regions are in different coordinate systems, predicate="within" matches nothing and every point becomes unmatched. Reprojecting both to the same equal-area CRS (here EPSG:5070) before the join removes that risk. The second is points on shared borders being matched to two polygons; how="left" with within keeps the first match, which is acceptable for cross-validation grouping as long as the assignment is deterministic.

How Leave-One-Region-Out Splits the Data

The diagram below shows four regions being cycled through as the held-out test set. In every fold, one region is fully withheld and the model trains on the union of the rest — never seeing a single sample from the region it will be scored on.

Leave-One-Region-Out: each region held out once Four rows, one per fold. In each fold a different region A, B, C or D is marked as the held-out test set while the remaining three regions form the training set, showing every region is tested exactly once. Four regions, four folds — each region is the test set once Region A Region B Region C Region D Fold 1 TEST train train train Fold 2 train TEST train train Fold 3 train train TEST train Fold 4 train train train TEST Every region serves as the held-out test set exactly once; no region ever appears in both train and test within a fold. Each score is a genuine transfer to an unseen region.

Production-Ready Code with a Per-Region Score Table

cross_val_score returns a bare array of numbers with no indication of which region produced which score — useless when you need to know where the model fails. Use cross_validate and pair each fold with its held-out region id to build an interpretable table.

import numpy as np
import pandas as pd
from sklearn.model_selection import LeaveOneGroupOut, cross_validate
from sklearn.ensemble import RandomForestRegressor


def leave_one_region_out_report(
    estimator,
    X: np.ndarray,
    y: np.ndarray,
    groups: np.ndarray,
    scoring: str = "r2",
) -> pd.DataFrame:
    """Run leave-one-region-out CV and return a per-region score table.

    Parameters
    ----------
    estimator : fitted-capable sklearn estimator
    X : feature matrix, shape (n_samples, n_features)
    y : target vector, shape (n_samples,)
    groups : region id per sample, shape (n_samples,)
    scoring : any sklearn scoring string

    Returns
    -------
    pandas.DataFrame indexed by region, with test size and score,
    plus a size-weighted aggregate row.
    """
    groups = np.asarray(groups)
    if groups.shape[0] != X.shape[0]:
        raise ValueError("groups length must match number of samples.")

    logo = LeaveOneGroupOut()
    # The held-out region for fold i is the unique group in its test_idx
    held_out = [
        np.unique(groups[test_idx])[0]
        for _, test_idx in logo.split(X, y, groups)
    ]
    sizes = [np.sum(groups == r) for r in held_out]

    cv = cross_validate(
        estimator, X, y,
        groups=groups, cv=logo,
        scoring=scoring, n_jobs=-1,
        return_train_score=False,
    )

    table = pd.DataFrame({
        "region": held_out,
        "n_test": sizes,
        scoring: cv["test_score"],
    }).sort_values(scoring).reset_index(drop=True)

    # Size-weighted mean is more honest than an unweighted fold average
    weighted = np.average(table[scoring], weights=table["n_test"])
    table.loc[len(table)] = ["WEIGHTED_MEAN", table["n_test"].sum(), weighted]
    return table


rf = RandomForestRegressor(n_estimators=300, random_state=42, n_jobs=-1)
report = leave_one_region_out_report(rf, X, y, region_id, scoring="r2")
print(report.to_string(index=False))

A representative table makes the deployment risk legible at a glance:

        region  n_test     r2
     ecoreg_23     118 -0.041
     ecoreg_08     342  0.211
     ecoreg_15     890  0.517
     ecoreg_11    1204  0.603
 WEIGHTED_MEAN    2554  0.488

The negative score on ecoreg_23 is the headline result: the model does worse than predicting the mean when transferred there. That region is either genuinely out-of-distribution or under-sampled, and no unweighted average would have surfaced it.

Handling Imbalanced Group Sizes

Real regions are never equal in sample count, and that imbalance distorts naive aggregation two ways. First, an unweighted mean lets a 118-sample region swing the headline number as much as a 1204-sample region. The size-weighted mean in the function above corrects this. Second, tiny regions produce high-variance, sometimes meaningless scores — an R² computed on 12 samples is noise. Before assigning groups, merge regions below a minimum sample threshold into their nearest neighbour, or drop them from the transfer evaluation and note the reduced coverage.

When you have many small regions and do not need a fold per region, switch to GroupKFold, which packs whole regions into a fixed n_splits without ever splitting a region across the boundary. It keeps the transfer semantics while capping compute and stabilizing fold sizes. For choosing the spatial granularity of groups in the first place, the reasoning in Choosing Block Size for Spatial Block Cross-Validation applies directly: too-fine regions leak, too-coarse regions leave too little to train on.

Verification

Two invariants must hold for the split to be a valid region-held-out evaluation: each region is held out exactly once, and no region appears in both train and test within any fold. Assert them explicitly rather than trusting the API.

import numpy as np
from collections import Counter
from sklearn.model_selection import LeaveOneGroupOut


def verify_leave_one_region_out(X, groups) -> None:
    """Assert LORO invariants: each region tested once, no overlap per fold."""
    groups = np.asarray(groups)
    logo = LeaveOneGroupOut()
    held_out_counts = Counter()

    for train_idx, test_idx in logo.split(X, groups=groups):
        train_regions = set(groups[train_idx])
        test_regions = set(groups[test_idx])

        # Invariant 1: no region in both train and test
        overlap = train_regions & test_regions
        assert not overlap, f"Region(s) {overlap} leaked across the split"

        # Each fold holds out exactly one region
        assert len(test_regions) == 1, "A fold held out more than one region"
        held_out_counts[test_regions.pop()] += 1

    # Invariant 2: every region held out exactly once
    assert all(c == 1 for c in held_out_counts.values()), \
        "Some region was held out more or less than once"
    assert set(held_out_counts) == set(np.unique(groups)), \
        "Not every region was used as a test fold"
    print(f"Verified: {len(held_out_counts)} regions, each held out exactly once.")


verify_leave_one_region_out(X, region_id)

Once the split is proven sound, pair the per-region scores with a residual audit. The techniques in Evaluating Model Performance with Spatial Metrics let you check whether the residuals in a poorly-transferring region are spatially structured — the difference between a model that is merely uncalibrated there and one that is missing a driver entirely.

FAQ

When should I use LeaveOneGroupOut instead of GroupKFold?

Use LeaveOneGroupOut when you want one fold per region and a per-region score table — it holds out exactly one group at a time, giving you as many folds as you have regions. Use GroupKFold when you have many small regions and want to cap the number of folds at n_splits; it packs whole regions into k folds without splitting any region across the train/test boundary. Both guarantee no region appears in both training and test, so the transfer semantics are preserved either way.

How do I assign region group labels to my samples?

Perform a spatial join between your sample points and a polygon layer of regions (watersheds, ecoregions, or administrative units) with geopandas.sjoin using predicate="within". The joined region identifier becomes the groups array passed to LeaveOneGroupOut. Ensure both layers share the same CRS before joining, and handle points that fall outside every polygon so they do not become a silent NaN group.

What if my regions have very different numbers of samples?

Imbalanced group sizes are expected with real regions. Do not average fold scores with equal weight — a region with 20 samples should not count the same as one with 2000. Report a per-region score table, weight the aggregate mean by test-set size, and inspect the worst-performing regions individually. If a region is too small to yield a stable score, merge it with an adjacent region before assigning groups.

Part of: Spatial Cross-Validation Strategies — Training Geospatial Predictive Models in Python