My Model Scores 0.95 in Testing and Fails in the Field

Problem statement

A model with an excellent cross-validation score produces a map that domain experts reject. The score was not a lie — it answered a different question from the one deployment asks.

Measured on a real problem — Sentinel-2 NDVI on Snowdon predicted from terrain by a random forest on 3,000 cells — scored three ways (Examples 2 and 4):

random 5-fold CV, inside the sampled area      +0.374
true error over the whole unsampled map        +0.348
the same model applied to the adjacent half    -0.409

Random CV was within 0.03 R² for the task it matched, and a different task — trained in the western half of the window, applied to the eastern — gave a negative R², worse than predicting the training mean.

Quick answer

Work through five checks, in order of how often they are the cause:

1. is the deployment task the same as the CV task?
2. do the covariates in the new area resemble the training range?
3. did anything leak — coordinates, a proxy for the target, duplicated rows?
4. is the target distribution the same?
5. is the model extrapolating, and how does it behave when it does?

for i, name in enumerate(feature_names):
    train, new = X_train[:, i], X_new[:, i]
    lo, hi = np.percentile(train, [1, 99])
    outside = float(((new < lo) | (new > hi)).mean())
    print(f"  {name:12} {outside:6.1%} outside the training range")

  elevation     12.0% outside the training range
  slope          3.3% outside the training range
  northness      2.3% outside the training range
  eastness       1.8% outside the training range
  tpi            5.3% outside the training range

CV task against deployment, shift, leakage, extrapolation

1. Compare the CV task with the deployment task

The most common cause, and the easiest to miss because both are called "prediction".

If you validated by holding out random rows and you deploy to a new region, you measured interpolation and deployed extrapolation. The measured gap between those two on one dataset was +0.374 against −0.409.

Note that block CV does not close it: on the same western training cells it promised +0.033 (measured here).

2. Check covariate shift

A model can only interpolate the feature space it saw. If the new area's elevations, slopes or spectral values lie outside the training range, every prediction there is extrapolation.

A feature with more than about 10% of its new values outside the training 1st-to-99th percentile band is a strong signal, and tree models handle it particularly badly. For the eastern half, 12.0% of elevations were outside the band learnt in the west; every other terrain feature was under 6% (Example 1).

3. Look for leakage

Three kinds, in decreasing order of subtlety:

• Spatial — random splits put near neighbours in both folds.
• Target proxy — a covariate derived from the target, or measured at the same time by the same process.
• Duplicated rows — repeated samples from one site, which appear in both folds under any splitter.

The tell for a proxy is a single feature dominating importance and a score that is too good. A soil model with R² = 0.98 from remote sensing alone is usually predicting the survey campaign rather than the soil.

4. Compare the target distributions

print(f"train: mean {y_train.mean():.3f}, sd {y_train.std():.3f}")
print(f"new:   mean {y_new.mean():.3f}, sd {y_new.std():.3f}")

train: mean 0.745, sd 0.145
new:   mean 0.713, sd 0.145

If the new area's target has a different mean or spread, the model was fitted on a different problem. That is not fixable by tuning.

5. Understand how your model extrapolates

A random forest cannot predict outside the range of the training target — it averages training values in a leaf. In a new area it saturates at a plausible-looking constant, producing a flat map rather than an obviously wrong one.

Measured, adding coordinates as features made the transfer less wrong (−0.081 against −0.409) precisely because every new x value was outside the training range, so the model flattened: its eastern predictions had a standard deviation of 0.087, against 0.134 without coordinates and 0.145 in the data (Example 4). A more "conservative" model, by accident.

Linear models do the opposite: they extrapolate confidently and diverge.

A deployment audit and an applicability mask

The data used on this page

A 10 km window on Snowdon, 4.16–4.00° W, 53.03–53.12° N. The target is NDVI from the least cloudy Sentinel-2 L2A scene of June–September 2025 on Earth Search, on the scene's 20 m grid; the features are terrain from the Copernicus DEM GLO-30 resampled to the same grid — elevation, slope, northness, eastness and a topographic position index over 300 m. Cells the scene classifies as water, cloud, cloud shadow, snow or no data are dropped. X, y and coords hold every remaining cell. The deployment the page audits is a model trained on 3,000 random cells of the western half (X_train, y_train) and applied to every cell of the eastern half (X_new, y_new).

import numpy as np
import rasterio
import rasterio.merge
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.warp import reproject, transform_bounds
from scipy import ndimage

AOI = (-4.16, 53.03, -4.00, 53.12)
items = Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=list(AOI), datetime="2025-06-01/2025-10-01",
    query={"eo:cloud_cover": {"lt": 60}}).items()
item = min(items, key=lambda i: i.properties["eo:cloud_cover"])

with rasterio.open(item.assets["scl"].href) as ds:                 # the 20 m grid
    win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
                                       transform=ds.transform).round_offsets().round_lengths()
    scl, transform, crs = ds.read(1, window=win), ds.window_transform(win), ds.crs
    bounds = rasterio.windows.bounds(win, ds.transform)
shape = scl.shape
band = {}
for name in ("red", "nir"):                                         # 10 m bands averaged onto it
    with rasterio.open(item.assets[name].href) as ds:
        band[name] = ds.read(1, window=rasterio.windows.from_bounds(*bounds, transform=ds.transform),
                             out_shape=shape, resampling=Resampling.average).astype("float32") / 1e4
ndvi = (band["nir"] - band["red"]) / (band["nir"] + band["red"])

tiles = [rasterio.open(f"https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{t}_00_DEM/"
                       f"Copernicus_DSM_COG_10_N53_00_{t}_00_DEM.tif") for t in ("W005", "W004")]
mosaic, mosaic_transform = rasterio.merge.merge(tiles, bounds=(-4.18, 53.02, -3.98, 53.13))
elev = np.full(shape, np.nan, dtype="float32")
reproject(mosaic[0], elev, src_transform=mosaic_transform, src_crs=tiles[0].crs,
          dst_transform=transform, dst_crs=crs, resampling=Resampling.bilinear)
dz_dy, dz_dx = np.gradient(elev, -20, 20)                           # rows run south
slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
aspect = np.arctan2(-dz_dx, -dz_dy)                                 # direction the slope faces, from north
tpi = elev - ndimage.uniform_filter(elev, size=15)                  # against the 300 m neighbourhood

valid = np.isin(scl, (2, 4, 5, 7)) & np.isfinite(ndvi) & np.isfinite(elev)
rows, cols = np.nonzero(valid)
xs, ys = rasterio.transform.xy(transform, rows, cols)
coords = np.column_stack([xs, ys])
X = np.column_stack([elev[valid], slope[valid], np.cos(aspect[valid]), np.sin(aspect[valid]), tpi[valid]])
y = ndvi[valid]
FEATURES = ["elevation", "slope", "northness", "eastness", "tpi"]
west = coords[:, 0] < np.median(coords[:, 0])
pick = np.random.default_rng(0).choice(np.nonzero(west)[0], 3000, replace=False)
X_train, y_train, X_new, y_new = X[pick], y[pick], X[~west], y[~west]
feature_names = FEATURES
print(f"{item.id}: {shape[1]} x {shape[0]} cells of 20 m, {valid.sum():,} usable "
      f"({valid.mean():.1%}); {len(X_train):,} training cells in the west, {len(X_new):,} to predict in the east")

S2C_30UVD_20250712_0_L2A: 543 x 509 cells of 20 m, 270,536 usable (97.9%); 3,000 training cells in the west, 135,292 to predict in the east

Example 1 — a deployment audit

import numpy as np
from scipy.stats import ks_2samp


def deployment_audit(X_train, X_new, feature_names, y_train=None, y_new=None):
    """Everything that differs between training and deployment."""
    print(f"  {len(X_train):,} training rows, {len(X_new):,} deployment rows")

    problems = []
    for i, name in enumerate(feature_names):
        train, new = X_train[:, i], X_new[:, i]
        lo, hi = np.percentile(train, [1, 99])
        outside = float(((new < lo) | (new > hi)).mean())
        shift = float((new.mean() - train.mean()) / (train.std() + 1e-12))
        statistic, p_value = ks_2samp(train, new)

        flag = ""
        if outside > 0.10:
            flag = "  <- extrapolating"
            problems.append(f"{name}: {outside:.0%} outside the training range")
        elif abs(shift) > 1.0:
            flag = "  <- shifted"
            problems.append(f"{name}: mean shifted {shift:+.1f} sd")

        print(f"  {name:14} {outside:6.1%} outside  "
              f"{shift:+6.2f} sd shift  KS {statistic:.3f}{flag}")

    if y_train is not None and y_new is not None:
        print(f"  target: train mean {y_train.mean():.3f} "
              f"sd {y_train.std():.3f}; "
              f"new mean {y_new.mean():.3f} sd {y_new.std():.3f}")

    for p in problems:
        print(f"  ! {p}")
    return problems


deployment_audit(X_train, X_new, feature_names, y_train, y_new)

  3,000 training rows, 135,292 deployment rows
  elevation       12.0% outside   +0.50 sd shift  KS 0.177  <- extrapolating
  slope            3.3% outside   +0.58 sd shift  KS 0.257
  northness        2.3% outside   -0.08 sd shift  KS 0.071
  eastness         1.8% outside   +0.13 sd shift  KS 0.108
  tpi              5.3% outside   -0.00 sd shift  KS 0.102
  target: train mean 0.745 sd 0.145; new mean 0.713 sd 0.145
  ! elevation: 12% outside the training range

The Kolmogorov-Smirnov statistic catches distribution changes that a mean shift misses — a feature with the same mean and a different shape.

Example 2 — checking for leakage

import numpy as np
from sklearn.metrics import r2_score


def leakage_checks(X, y, coords, feature_names, model_factory):
    """Three tests that catch the common leaks."""
    # 1. a single dominant feature is suspicious
    model = model_factory().fit(X, y)
    importance = getattr(model, "feature_importances_", None)
    if importance is not None:
        order = np.argsort(-importance)
        top = importance[order[0]]
        print(f"  top feature: {feature_names[order[0]]} ({top:.1%})")
        if top > 0.6:
            print("  ! one feature dominates — check it is not derived from "
                  "the target")

    # 2. duplicate or near-duplicate coordinates
    from scipy.spatial import cKDTree
    distances, _ = cKDTree(coords).query(coords, k=2)
    duplicates = int((distances[:, 1] < 1e-6).sum())
    if duplicates:
        print(f"  ! {duplicates:,} rows share a location with another row — "
              "these will land in different folds")

    # 3. how much better is random CV than a spatial split?
    from sklearn.cluster import KMeans
    from sklearn.model_selection import GroupKFold, KFold

    def score(splitter):
        predicted = np.empty(len(y))
        for train, test in splitter:
            predicted[test] = model_factory().fit(X[train], y[train]) \
                                             .predict(X[test])
        return r2_score(y, predicted)

    random = score(KFold(5, shuffle=True, random_state=0).split(X))
    blocks = KMeans(10, random_state=0, n_init=10).fit_predict(coords)
    spatial = score(GroupKFold(5).split(X, y, blocks))

    print(f"  random CV {random:+.3f}, block CV {spatial:+.3f}, "
          f"gap {random - spatial:+.3f}")
    if random - spatial > 0.3:
        print("  ! large gap — the model relies heavily on spatial proximity")
    return {"random": random, "block": spatial}


from sklearn.ensemble import RandomForestRegressor

forest = lambda: RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1, random_state=0)
sample = np.random.default_rng(0).choice(len(y), 3000, replace=False)
leakage_checks(X[sample], y[sample], coords[sample], feature_names, forest)

  top feature: elevation (49.0%)
  random CV +0.374, block CV +0.258, gap +0.116

A large random-block gap is not proof of a problem, and a small one is not an all-clear: here it was 0.116, on a model whose random CV (+0.374) was close to its true in-area score (+0.348) and which still failed in the other half of the window. It is a signal that the model leans on location, which matters if you will deploy elsewhere and does not if you will not.

Example 3 — an applicability mask to ship with the map

import numpy as np
from scipy.spatial import cKDTree
from sklearn.preprocessing import StandardScaler


def applicability_mask(X_train, X_new, quantile=0.95, k=5):
    """Flag predictions whose covariates are unlike anything in training."""
    scaler = StandardScaler().fit(X_train)
    train_scaled = scaler.transform(X_train)
    new_scaled = scaler.transform(X_new)

    tree = cKDTree(train_scaled)
    train_distances, _ = tree.query(train_scaled, k=k + 1)
    threshold = float(np.quantile(train_distances[:, 1:].mean(axis=1), quantile))

    new_distances, _ = tree.query(new_scaled, k=k)
    new_mean = new_distances.mean(axis=1)
    inside = new_mean <= threshold

    print(f"  threshold: {quantile:.0%} of training points are within "
          f"{threshold:.3f} (scaled feature-space distance)")
    print(f"  {inside.mean():.1%} of deployment locations are inside the "
          "area of applicability")
    return inside, new_mean


inside, distance = applicability_mask(X_train, X_new)
east_model = forest().fit(X_train, y_train).predict(X_new)
for label, part in (("inside", inside), ("outside", ~inside)):
    print(f"  {label} the mask: R2 {r2_score(y_new[part], east_model[part]):+.3f}, "
          f"median error {np.median(np.abs(y_new[part] - east_model[part])):.3f} over {part.sum():,} cells")

  threshold: 95% of training points are within 0.946 (scaled feature-space distance)
  74.5% of deployment locations are inside the area of applicability
  inside the mask: R2 -0.891, median error 0.067 over 100,796 cells
  outside the mask: R2 -0.342, median error 0.117 over 34,496 cells

Shipping this mask with the prediction is the practical answer to "we cannot validate transfer", and it has a limit worth seeing. Here 74.5% of the east was inside it, and the errors there were smaller — a median of 0.067 against 0.117 outside. But R² was negative on both sides (−0.891 inside, −0.342 outside, each against its own variance): the model was worse than the mean everywhere in the east. The mask marks where the covariates are unfamiliar; it cannot mark where the relationship between them and the target has changed.

Example 4 — the three scores, and what coordinates do

from sklearn.model_selection import KFold

rest = np.setdiff1d(np.arange(len(y)), sample)
with_xy = np.column_stack([X, coords])
for label, features in (("terrain only", X), ("terrain + x, y", with_xy)):
    Xs, ys = features[sample], y[sample]
    p = np.empty(len(ys))
    for train, test in KFold(5, shuffle=True, random_state=0).split(Xs):
        p[test] = forest().fit(Xs[train], ys[train]).predict(Xs[test])
    in_area = r2_score(y[rest], forest().fit(Xs, ys).predict(features[rest]))
    moved = forest().fit(features[pick], y[pick]).predict(features[~west])
    print(f"{label:15} random CV {r2_score(ys, p):+.3f}   true in-area {in_area:+.3f}   "
          f"west -> east {r2_score(y_new, moved):+.3f}   east prediction sd {moved.std():.3f} "
          f"(the east itself {y_new.std():.3f})")

terrain only    random CV +0.374   true in-area +0.348   west -> east -0.409   east prediction sd 0.134 (the east itself 0.145)
terrain + x, y  random CV +0.504   true in-area +0.509   west -> east -0.081   east prediction sd 0.087 (the east itself 0.145)

How a score can be right and useless

Why the score can be exactly right and still useless

Cross-validation estimates generalisation to data drawn from the same distribution as the training data. That is a precise statement and a narrow one.

Random CV on a dense random sample estimates interpolation within the sampled extent, and it does so well — +0.374 against a true +0.348 in the measurement here, and +0.339 against +0.344 averaged over ten samples.

Deployment to a new region asks about a different distribution. No resampling of the training data can estimate it, because the training data contains no information about the difference.

Why tree models fail quietly

A random forest predicts the mean of the training targets in a leaf. It therefore cannot produce a value outside the training target range, and for a feature value beyond the training range it reuses the most extreme split it learned.

Predictions in a new area consequently converge on a constant. The map looks plausible — smooth, in range, no obvious artefacts — and carries no information.

That is worse than an obviously wrong map, which at least gets noticed. It is why an applicability mask matters more for tree models than for linear ones.

Why coordinates make transfer worse and the score better

Adding raw x and y let the model carve the study area into regions with different local means. That is genuinely useful for interpolation — the measured in-area R² rose from 0.348 to 0.509.

It is a lookup table for one place. Applied elsewhere, the splits are meaningless.

The measured transfer results are instructive: without coordinates, −0.409; with them, −0.081. The version with coordinates scored better because it had been flattened towards a constant, not because it transferred.

Why domain experts spot it first

A model's failure in a new area is usually a pattern failure: the map is smooth where it should be varied, or puts a class where everyone knows it does not occur.

Metrics do not capture pattern. An expert looking at the map sees it immediately, and the fastest route to catching this class of failure is to show a domain expert the map before the metrics.

Range overlap, saturation and target proxies

• A good score for the wrong task is still a good score. State the task.
• Check covariate range overlap before deploying anywhere new.
• Tree models saturate outside their range; linear models diverge.
• Coordinates as features destroy transferability while improving in-area accuracy.
• Duplicated sites leak under any splitter.
• A dominant feature may be a target proxy.
• A negative R² means worse than the mean — report it.
• Show the map to a domain expert before trusting the metric.

Internal links

• Spatial leakage explained: why random cross-validation lies — the measurements behind this
• How to do spatial cross-validation in Python — choosing a scheme
• Extrapolation in space explained: the area of applicability — the mask to ship
• Spatial machine learning explained — why location breaks the usual rules
• Feature importance says coordinates are the best predictor — the coordinate trap
• How to evaluate a spatial model honestly — what to report
• My prediction raster is striped, blocky or full of NoData — the visual symptoms
• Spatial autocorrelation explained — the structure underneath

FAQ

Why does my model score well and fail in practice?

Almost always because the validation task and the deployment task differ. Random CV measures interpolation; deploying to a new region is extrapolation, and the measured gap was +0.374 against −0.409.

Does spatial cross-validation fix it?

Not for a new region. Block CV promised +0.033 where a real transfer test gave −0.409.

How do I know if I am extrapolating?

Compare each covariate's deployment values against the training 1st-to-99th percentile band. More than about 10% outside means extrapolation.

Why does my prediction map look flat in the new area?

A tree model cannot predict outside its training target range, so beyond the training feature space it saturates at a constant.

Should I include coordinates as features?

For interpolation inside the sampled area, they help — 0.348 to 0.509 here. For anything that must transfer, no.

What is a target proxy?

A covariate derived from, or measured alongside, the target. It produces implausibly high scores and does not exist at prediction time.

What should I ship with a prediction map?

An area-of-applicability mask showing where the covariates resemble the training data, and the validation scheme with the score.