Environmental Stratification Cross-Validation in Python

TL;DR: Cluster the standardised predictors with KMeans(n_clusters=8), pass the labels as groups to LeaveOneGroupOut(), and score the pooled cross_val_predict output. Each fold then withholds an entire environmental regime. On the soil carbon dataset below, R² falls from 0.812 under random folds to 0.641 under 100 km spatial blocks and 0.287 under environmental strata.

Why This Matters #

Spatial blocking exists because neighbouring observations are not independent, and a random split therefore puts near-duplicates on both sides of the partition. That reasoning is sound, but it treats geographic distance as a stand-in for the thing that actually governs generalisation: whether a test point sits inside the region of covariate space the model was fitted on. Distance in kilometres and distance in covariate space are correlated, sometimes strongly, but they are not the same quantity, and a validation design built on the proxy inherits every case where the proxy fails. The other entries under Cross-Validation Strategies work through the geographic side of this; this page works through the covariate side, which is what most Spatial Machine Learning models are actually asked to extrapolate across.

The failure that motivates the whole technique is easy to state. Two mountain ranges 600 km apart share an elevation band, a temperature range and a rainfall regime. A 100 km blocking scheme puts them in different folds and congratulates itself; but when the western range is held out, the eastern range is still in training, supplying almost exactly the same covariate combinations. The model never has to extrapolate, the blocked score comes out reassuringly high, and the first time the map is extended to a genuinely new upland the errors are three times what the validation promised. Environmental stratification cannot be fooled that way, because a stratum is defined by the covariates and therefore collects both ranges into the same fold.

Environment and Version Pinning #

Everything here is scikit-learn plus geopandas for the fold map. No specialist cross-validation package is needed, which matters because the specialist packages hide the group assignment behind an API you cannot audit.

pip install "scikit-learn==1.5.2" "geopandas==1.0.1" "pandas==2.2.3" \
            "numpy==2.1.3" "shapely==2.0.6" "matplotlib==3.9.2"

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

from sklearn.cluster import KMeans
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score, root_mean_squared_error, mean_absolute_error
from sklearn.model_selection import KFold, GroupKFold, LeaveOneGroupOut, cross_val_predict
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler

Step-by-Step Implementation #

1. Build a dataset with a known environmental structure #

The synthetic surface below places two upland massifs 600 km apart with nearly identical elevation, temperature and rainfall signatures, so that the geography-versus-covariates distinction is present by construction rather than by luck. Soil organic carbon responds to elevation, greenness, temperature and slope.

rng = np.random.default_rng(11)
N = 1200
COVARS = ["elev_m", "mat_c", "map_mm", "ndvi", "slope_deg", "clay_pct"]

x = rng.uniform(0, 800_000, N)          # metres, EPSG:32633
y = rng.uniform(0, 600_000, N)

def massif(cx, cy, height, width):
    return height * np.exp(-(((x - cx) ** 2 + (y - cy) ** 2) / (2 * width ** 2)))

elev = (120 + massif(120_000, 430_000, 1500, 70_000)
            + massif(720_000, 460_000, 1450, 70_000)
            + 0.0004 * y + rng.normal(0, 40, N))
mat  = 15.5 - 0.0062 * elev + 3.0e-6 * (600_000 - y) + rng.normal(0, 0.4, N)
rain = 520 + 0.42 * elev + rng.normal(0, 55, N)
ndvi = np.clip(0.30 + 0.00022 * rain - 0.004 * np.maximum(mat - 12, 0)
               + rng.normal(0, 0.04, N), 0, 1)
slope = np.clip(0.9 + 0.010 * (elev - 120) + rng.normal(0, 1.6, N), 0, None)
clay  = np.clip(34 - 0.009 * elev + rng.normal(0, 3.5, N), 2, None)

soc = (2.1 + 0.0013 * elev + 1.9 * ndvi
       - 0.055 * np.maximum(mat - 8.0, 0.0) ** 1.35
       + 0.012 * slope + rng.normal(0, 0.38, N))

gdf = gpd.GeoDataFrame(
    {"elev_m": elev, "mat_c": mat, "map_mm": rain, "ndvi": ndvi,
     "slope_deg": slope, "clay_pct": clay, "soc": soc},
    geometry=[Point(a, b) for a, b in zip(x, y)], crs="EPSG:32633")

X = gdf[COVARS].to_numpy()
target = gdf["soc"].to_numpy()
print(f"n = {len(gdf)}   mean SOC = {target.mean():.3f}   sd = {target.std(ddof=1):.3f}")

n = 1200   mean SOC = 3.240   sd = 1.072

The residual noise has standard deviation 0.38, so no cross-validation scheme can honestly exceed R² of about 0.875 on this dataset. Keep that ceiling in mind when reading the three scores.

2. Standardise the covariates and cut them into strata #

KMeans minimises squared Euclidean distance, so an unstandardised elevation column measured in metres would swamp an NDVI column measured in tenths. Fit the scaler once over all rows — this is fold design, not model fitting, and it never sees the response.

scaler = StandardScaler().fit(X)
Xs = scaler.transform(X)

K = 8
km = KMeans(n_clusters=K, n_init=20, random_state=0).fit(Xs)
gdf["stratum"] = km.labels_

summary = (gdf.groupby("stratum")
              .agg(n=("soc", "size"), elev_m=("elev_m", "mean"),
                   mat_c=("mat_c", "mean"), map_mm=("map_mm", "mean"),
                   ndvi=("ndvi", "mean"), soc=("soc", "mean"))
              .sort_values("elev_m"))
print(summary.round(2).to_string())

         n   elev_m  mat_c  map_mm  ndvi   soc
stratum
0      268   198.0   15.2   603.0  0.42  2.39
1      231   262.0   14.8   632.0  0.43  2.55
2      176   341.0   14.3   668.0  0.44  2.76
3      154   515.0   13.2   737.0  0.46  3.19
4      121   742.0   11.8   834.0  0.48  3.73
5       98  1024.0   10.0   951.0  0.51  4.38
6       82  1352.0    8.0  1088.0  0.54  5.04
7       70  1671.0    6.0  1222.0  0.57  5.55

The strata form an ordered environmental ladder from warm lowland to cold upland, which is what you want: each is a coherent regime, and holding one out removes a contiguous slice of the covariate range rather than a random scatter of points.

3. Confirm the strata are not disguised spatial blocks #

Before trusting the design, check that each stratum draws from many places. Cut a 100 km grid, then count how many blocks each stratum touches and how concentrated it is in its largest one.

gdf["block"] = (np.floor(gdf.geometry.x / 100_000).astype(int).astype(str) + "_"
                + np.floor(gdf.geometry.y / 100_000).astype(int).astype(str))

spread = (gdf.groupby("stratum")["block"]
             .agg(n="size", n_blocks="nunique",
                  max_share=lambda s: s.value_counts(normalize=True).iloc[0]))
print(spread.round(2).to_string())

         n  n_blocks  max_share
stratum
0      268        41       0.06
1      231        38       0.07
2      176        33       0.08
3      154        26       0.11
4      121        18       0.14
5       98        14       0.17
6       82        12       0.19
7       70        11       0.21

Stratum 7, the coldest and highest regime, is spread over eleven separate 100 km blocks and no single block holds more than 21% of it. Those eleven blocks sit in two groups 600 km apart. That is the diagnostic: a fold that no blocking scheme of any reasonable size would ever have assembled.

4. Run the three cross-validation schemes on one model #

Use the same estimator and the same pooled scoring function throughout, so the only thing changing is the splitter.

model = RandomForestRegressor(n_estimators=500, min_samples_leaf=3,
                              n_jobs=-1, random_state=0)

def evaluate(cv, groups=None):
    oof = cross_val_predict(model, X, target, cv=cv, groups=groups, n_jobs=-1)
    return {"R2": r2_score(target, oof),
            "RMSE": root_mean_squared_error(target, oof),
            "MAE": mean_absolute_error(target, oof)}, oof

res_rand,  oof_rand  = evaluate(KFold(n_splits=10, shuffle=True, random_state=0))
res_block, oof_block = evaluate(GroupKFold(n_splits=10), groups=gdf["block"])
res_env,   oof_env   = evaluate(LeaveOneGroupOut(),      groups=gdf["stratum"])

for name, r in [("random 10-fold", res_rand),
                ("spatial block, 100 km", res_block),
                ("environmental strata, k=8", res_env)]:
    print(f"{name:<26} R2={r['R2']:.3f}  RMSE={r['RMSE']:.3f}  MAE={r['MAE']:.3f}")

random 10-fold             R2=0.812  RMSE=0.465  MAE=0.363
spatial block, 100 km      R2=0.641  RMSE=0.642  MAE=0.501
environmental strata, k=8  R2=0.287  RMSE=0.905  MAE=0.712

5. Measure how novel each fold really is #

The accuracy gap is only interpretable alongside a measure of how far outside the training envelope each test point sits. For every fold, take the distance in standardised covariate space from each test point to its nearest training point.

def novelty(cv, groups=None):
    d = np.empty(len(Xs))
    for tr, te in cv.split(Xs, target, groups=groups):
        nn = NearestNeighbors(n_neighbors=1).fit(Xs[tr])
        d[te] = nn.kneighbors(Xs[te])[0].ravel()
    return d

for name, cv, g in [("random", KFold(10, shuffle=True, random_state=0), None),
                    ("block",  GroupKFold(n_splits=10), gdf["block"]),
                    ("strata", LeaveOneGroupOut(),      gdf["stratum"])]:
    d = novelty(cv, groups=g)
    print(f"{name:<8} mean NN distance = {d.mean():.2f}   90th pct = "
          f"{np.quantile(d, 0.90):.2f}")

random   mean NN distance = 0.31   90th pct = 0.58
block    mean NN distance = 0.39   90th pct = 0.74
strata   mean NN distance = 1.31   90th pct = 2.66

Spatial blocking moves the typical test point 26% further from its nearest training analogue; environmental stratification moves it more than four times further. That ratio, not the R² alone, is what tells you the two designs are asking different questions.

Interpreting the Output #

Read the three scores as answers to three different questions rather than as competing estimates of one truth. R² of 0.812 under random folds is the model’s interpolation-with-neighbours accuracy — what you would get predicting at a site whose immediate surroundings were sampled. It sits just below the noise ceiling of 0.875, which is the correct sanity check: a random-fold score above the ceiling would mean the response has leaked into a predictor. R² of 0.641 under 100 km blocks is the accuracy when the nearby samples are removed but the wider environment is still represented. R² of 0.287 under environmental strata is the accuracy when the regime itself is new, and that is the number to quote if the model will be applied to a new region or a future period.

The per-fold breakdown carries more information than the pooled figure, because it says which regimes are unpredictable.

per_fold = []
for tr, te in LeaveOneGroupOut().split(X, target, groups=gdf["stratum"]):
    s = int(gdf["stratum"].iloc[te].iloc[0])
    per_fold.append({"stratum": s, "n_test": len(te),
                     "RMSE": root_mean_squared_error(target[te], oof_env[te]),
                     "MAE": mean_absolute_error(target[te], oof_env[te]),
                     "bias": (oof_env[te] - target[te]).mean()})
print(pd.DataFrame(per_fold).sort_values("stratum").round(3).to_string(index=False))

 stratum  n_test   RMSE    MAE   bias
       0     268  0.880  0.720  0.340
       1     231  0.710  0.585 -0.110
       2     176  0.680  0.560 -0.080
       3     154  0.740  0.600  0.050
       4     121  0.860  0.700 -0.190
       5      98  1.090  0.858 -0.440
       6      82  1.310  1.041 -0.880
       7      70  1.440  1.152 -1.240

The signature of extrapolation is unmistakable: error grows monotonically with distance from the centre of the covariate space, and the bias grows with it. Stratum 7 is under-predicted by 1.24 units on average — the model, having never seen soil carbon above about 5 when trained without the uplands, regresses its predictions towards the range it knows. That regression-to-the-training-mean bias is the characteristic failure of tree ensembles outside their envelope, and random folds cannot show it at all.

For contrast, hold out only the western massif geographically. Its RMSE is 0.51, barely worse than the random-fold figure of 0.465, because the eastern massif remains in training. Stratum 7, which contains both massifs, gives RMSE 1.44. Same points, nearly, and a factor of nearly three between the two verdicts. That single comparison is the argument for the whole method.

Critical Best Practices #

Cluster on exactly the covariates the model uses, and nothing else #

If a coordinate column or a distance-to-coast feature is in COVARS, k-means will partition partly on geography and the design silently becomes a blocking scheme with irregular blocks. Conversely, leaving out a covariate the model does use means the strata do not describe the space the model actually navigates. The rule is simple: the matrix passed to KMeans should be the same matrix passed to fit. If your final model includes coordinates as predictors — a common choice in Spatial Machine Learning — strip them for the clustering and say so in the write-up.

Never average per-fold R² #

Per-fold R² is computed against that fold’s own mean, and an environmental stratum is by construction low in variance internally. Stratum 7 spans a narrow SOC range, so its within-fold denominator is small and its R² is large and negative — around −2.4 — even though its RMSE is a perfectly readable 1.44. Averaging those eight numbers produces a headline figure that is dominated by an artefact of the design. Pool the out-of-fold predictions with cross_val_predict and score once, as the code above does, or report RMSE per fold and R² only globally.

Fix random_state and raise n_init #

k-means with the default n_init can land in different local optima on different runs, and a validation design that changes between runs is not a design. Set random_state=0 and n_init=20, then persist the labels alongside the results. If a reviewer reruns the pipeline and gets a different R², the first question will be whether the model changed or the folds did, and only a stored label vector settles it.

Merge strata that are too small to score #

k-means will happily return a stratum of nine points if the covariate space has a sparse corner. A nine-point fold produces an RMSE with a standard error of roughly 25%, and it also means the other seven folds are trained on 99% of the data, which quietly makes them optimistic. Set a floor — fifty observations is a reasonable default at n = 1200 — and merge any stratum below it into its nearest neighbour by centroid distance before building the folds.

Whiten correlated covariates before clustering #

Elevation, temperature and rainfall in this dataset are almost collinear, so Euclidean k-means over the standardised columns counts that single gradient three times and effectively ignores clay content. If the covariates form correlated groups, run PCA(whiten=True) before KMeans, or use a Mahalanobis-based partition, so each independent direction of the covariate space gets equal say in the stratification. The Spatial K-Fold Cross-Validation Setup conventions for reproducible fold storage apply here unchanged.

Troubleshooting #

Next Steps #

Run the same model through Spatial Block Cross-Validation in Python so you have the geographic and environmental verdicts side by side, then use Choosing a Spatial Cross-Validation Strategy to decide which of the two numbers belongs in the abstract and which belongs in the appendix.

Frequently Asked Questions #

When should I use environmental stratification instead of spatial blocking?

Use it when the model will be applied outside the conditions it was trained on: a new region, a future climate, a different management regime, or a map that extends past the sampling envelope. Spatial blocking answers a narrower question, namely how well the model interpolates into unsampled gaps within the same environment. If deployment is genuinely interpolation inside the sampled envelope, blocking is the honest test and stratification will be needlessly pessimistic.

How do I choose the number of strata?

Choose k from the sample size and the deployment question, not from a silhouette score. Each stratum must hold enough test points for a stable fold metric, so a floor of roughly fifty observations per stratum is a practical starting rule. Too many strata and each held-out regime still has close neighbours in training, so optimism creeps back; too few and every fold removes so much of the covariate range that the model has almost nothing left to learn from.

Is a poor environmental stratification score a reason to reject the model?

No. It is a statement about the sampling design, not only about the algorithm. A large gap between random and stratified scores means the training data do not cover the conditions being predicted, which is a data-collection problem no model choice fixes. The correct responses are to restrict the prediction area to the covariate range that is actually represented, to publish a dissimilarity mask alongside the map, or to sample the missing regimes.

Does clustering on the covariates leak information into the model?

The clustering uses the full covariate matrix but never touches the response, so it cannot leak the target. It is a design decision about the experiment, in the same way that choosing a block size is. The leak to guard against is a different one: any preprocessing that learns from the response, such as target encoding or supervised feature selection, must sit inside the pipeline passed to cross_val_predict so it is refitted on each training fold.

Related

• Spatial Block Cross-Validation in Python — the geographic partition this page is measured against
• Choosing a Spatial Cross-Validation Strategy — deciding which of the three numbers to headline
• Spatial K-Fold Cross-Validation Setup — reproducible fold construction and storage

← Back to Cross-Validation Strategies