Spatial Machine Learning Explained: Why Location Breaks the Usual Rules

Problem statement

Every introductory machine learning course assumes rows are independent and identically distributed. Spatial data violates that assumption in a specific, measurable way: nearby rows are similar, so a dataset of 3,000 cells does not contain 3,000 independent observations.

That single fact breaks three things people rely on:

• Cross-validation scores stop meaning what they usually mean — in either direction, depending on what you will predict.
• Feature importance starts rewarding whatever encodes location.
• Transfer to a new area becomes unpredictable from any internal statistic.

Measured on a real problem — predicting vegetation index from terrain over a 10 km window in Snowdonia:

trained and validated inside the sampled area   R²  +0.35
the same model applied to the adjacent half     R²  -0.41

The same model, the same features, ground 5 km away. Worse than predicting the mean.

Quick answer

Three habits handle most of it:

# 1. know which task you are doing, and validate that task
#    interpolation inside a sampled area  -> random k-fold is fine
#    extrapolation to new ground          -> hold out a real region

# 2. be deliberate about coordinates as features
X = terrain_features                      # transfers, weakly
X = np.column_stack([terrain, coords])    # fits better, transfers worse

# 3. report where the model is applicable, not just how good it is
applicable = feature_space_distance(X_new, X_train) < threshold

The task, the real sample size, coordinates, applicability

1. Establish what the model is for

There are three spatial prediction tasks, and they are not variations on one problem:

• Interpolation — fill gaps between samples in an area you have surveyed. This is the easy case, and ordinary random cross-validation estimates it well.
• Extrapolation in space — apply a model fitted here to ground over there. This is the hard case, and no internal validation estimates it.
• Extrapolation in time — apply a model fitted on past data to the future. Different again, and it needs a chronological split.

Almost every disappointment in applied spatial ML comes from validating the first and deploying the second.

2. Expect the sample size to be smaller than the row count

Neighbouring cells are highly correlated, so a dense sample contains far fewer independent observations than rows.

The measured consequence: going from 100 to 10,000 training cells improved the true map R² from +0.159 to +0.401, and the curve flattens fast — tripling from 3,000 to 10,000 cells bought 0.053 (Example 4). By one rough estimate, 3,000 cells in this 10 km window hold about eight independent looks at the landscape (Example 2).

That is the signature of correlated data. Extra rows from the same neighbourhoods add very little, and any confidence interval computed from the row count is far too narrow.

3. Treat coordinates as a deliberate choice

Adding raw x and y to a random forest improved the in-area map from R² 0.348 to 0.509 — a genuine improvement for interpolation, because coordinates proxy for unmeasured local structure.

The same features under block cross-validation improved much less, 0.258 to 0.283, because a coordinate learned here does not mean anything over there.

Neither "always include coordinates" nor "never include coordinates" is right. Include them for interpolation; drop them, or replace them with transferable covariates, for anything that must move.

4. Expect feature importance to shift

terrain only    elevation 0.490, slope 0.163, tpi 0.127,
                northness 0.113, eastness 0.108
terrain + x,y   elevation 0.443, y 0.118, slope 0.113, x 0.109,
                northness 0.080, tpi 0.077, eastness 0.059

Adding coordinates pushed y into second place and x into fourth, and pushed TPI from third down to sixth (Example 4). The physical story changed because two location proxies absorbed variance that the terrain variables had been explaining.

Importance is a statement about this fitted model on this data, not about the world.

5. Report applicability, not just accuracy

A model's accuracy is a single number about the data it saw. What a user of the map needs is where the prediction is supported — which is a raster, not a number.

Splits, effective sample size and covariate transfer on Snowdon

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 spatial train/test split that means something

import numpy as np
from sklearn.cluster import KMeans


def spatial_split(coords, test_fraction=0.3, n_blocks=10, seed=0):
    """Hold out whole blocks, so the test set is genuinely unseen ground."""
    blocks = KMeans(n_clusters=n_blocks, random_state=seed,
                    n_init=10).fit_predict(coords)

    rng = np.random.default_rng(seed)
    order = rng.permutation(n_blocks)
    sizes = np.array([(blocks == b).sum() for b in order])
    cumulative = np.cumsum(sizes) / sizes.sum()
    held_out = set(order[cumulative <= test_fraction].tolist())
    if not held_out:
        held_out = {int(order[0])}

    test = np.isin(blocks, list(held_out))
    print(f"  {len(held_out)} of {n_blocks} blocks held out, "
          f"{test.mean():.1%} of rows")
    return ~test, test


sample = np.random.default_rng(0).choice(len(y), 3000, replace=False)
train, test = spatial_split(coords[sample])

  2 of 10 blocks held out, 20.2% of rows

Ten blocks and a 30% hold-out is a reasonable default; because whole blocks are held out, the test share lands below the target — 20.2% here. Five blocks makes the score unstable — a measured standard deviation of 0.106 across samples (here) — while blocks smaller than the autocorrelation range, about 2.5 km on this data, start sharing information again.

Example 2 — how many independent observations do you have?

import numpy as np
from scipy.spatial import cKDTree


def effective_sample_size(coords, values, range_m):
    """Rows within one autocorrelation range of each other are not independent."""
    tree = cKDTree(coords)
    neighbours = tree.query_ball_point(coords, r=range_m)
    counts = np.array([len(n) for n in neighbours])

    effective = float((1.0 / counts).sum())
    print(f"  {len(coords):,} rows")
    print(f"  median neighbours within {range_m:.0f} m: {np.median(counts):.0f}")
    print(f"  effective sample size: {effective:,.0f} "
          f"({effective / len(coords):.1%} of the rows)")
    return effective


effective_sample_size(coords[sample], y[sample], range_m=2500)   # the NDVI variogram range

  3,000 rows
  median neighbours within 2500 m: 452
  effective sample size: 8 (0.3% of the rows)

This is a rough estimate, not a theorem: here each of the 3,000 cells had a median of 452 others within the 2.5 km range, which counts the sample as about eight independent looks. It is enough to stop you quoting a confidence interval based on 270,000 pixels when the landscape offers a handful of independent neighbourhoods.

Example 3 — checking the covariates transfer before you rely on them

import numpy as np


def covariate_shift(X_train, X_new, names):
    """Are the new area's covariates inside the range the model was fitted on?"""
    rows = []
    for i, name in enumerate(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))
        rows.append({"feature": name, "outside_train_range": round(outside, 4),
                     "mean_shift_sd": round(shift, 3)})
        flag = "  <-- extrapolating" if outside > 0.1 else ""
        print(f"  {name:12} {outside:6.1%} outside train range, "
              f"mean shifted {shift:+.2f} sd{flag}")
    return rows


covariate_shift(X_train, X_new, FEATURES)

  elevation     12.0% outside train range, mean shifted +0.50 sd  <-- extrapolating
  slope          3.3% outside train range, mean shifted +0.58 sd
  northness      2.3% outside train range, mean shifted -0.08 sd
  eastness       1.8% outside train range, mean shifted +0.13 sd
  tpi            5.3% outside train range, mean shifted -0.00 sd

Run this before applying a model anywhere new. A covariate with more than about 10% of its values in the new area outside the training range means the model is extrapolating — elevation, at 12.0%, for the eastern half here — and tree models handle that particularly badly: they cannot predict outside the range of the training target at all.

Example 4 — rows, coordinates and importance, measured

from sklearn.cluster import KMeans
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from sklearn.model_selection import GroupKFold

forest = lambda: RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1, random_state=0)

print("training cells -> true R2 over every other cell")
for n in (100, 300, 1000, 3000, 10000):
    pick = np.random.default_rng(0).choice(len(y), n, replace=False)
    others = np.setdiff1d(np.arange(len(y)), pick)
    print(f"  {n:6,} {r2_score(y[others], forest().fit(X[pick], y[pick]).predict(X[others])):+.3f}")

rest = np.setdiff1d(np.arange(len(y)), sample)
blocks = KMeans(10, random_state=0, n_init=10).fit_predict(coords[sample])
for label, features, names in (("terrain only", X, FEATURES),
                               ("terrain + x, y", np.column_stack([X, coords]), FEATURES + ["x", "y"])):
    Xs, ys = features[sample], y[sample]
    p = np.empty(len(ys))
    for tr, te in GroupKFold(5).split(Xs, ys, blocks):
        p[te] = forest().fit(Xs[tr], ys[tr]).predict(Xs[te])
    fitted = forest().fit(Xs, ys)
    order = np.argsort(-fitted.feature_importances_)
    print(f"{label:15} in-area R2 {r2_score(y[rest], fitted.predict(features[rest])):+.3f}, "
          f"10-block CV {r2_score(ys, p):+.3f}")
    print("    " + ", ".join(f"{names[i]} {fitted.feature_importances_[i]:.3f}" for i in order))

north = coords[:, 1] >= 5_880_000
print(f"NDVI north of 5,880,000: mean {y[north].mean():.3f} over {north.sum():,} cells; "
      f"south: {y[~north].mean():.3f} over {(~north).sum():,}")

training cells -> true R2 over every other cell
     100 +0.159
     300 +0.219
   1,000 +0.313
   3,000 +0.348
  10,000 +0.401
terrain only    in-area R2 +0.348, 10-block CV +0.258
    elevation 0.490, slope 0.163, tpi 0.127, northness 0.113, eastness 0.108
terrain + x, y  in-area R2 +0.509, 10-block CV +0.283
    elevation 0.443, y 0.118, slope 0.113, x 0.109, northness 0.080, tpi 0.077, eastness 0.059
NDVI north of 5,880,000: mean 0.705 over 165,467 cells; south: 0.769 over 105,069

Correlated rows, coordinates and tree extrapolation

Why correlated rows inflate confidence rather than accuracy

Two cells twenty metres apart have almost the same elevation, slope and vegetation. As training data the second one adds almost nothing; as a count it doubles your apparent sample.

Every statistic computed from the row count is therefore too confident: standard errors too small, confidence intervals too narrow, p-values too significant. Accuracy itself is unaffected — the model is as good as it is — but every statement about how certain you are is wrong.

That is why the measured curve flattens: 100 to 10,000 training cells moved the true R² from 0.159 to 0.401, and the last 7,000 of those rows added 0.053. They were mostly repeats of information already present.

Why coordinates help and hurt at the same time

A random forest given x and y can carve the study area into regions and fit a different local mean in each. That is a legitimate and powerful way to capture spatial structure the measured covariates miss — soil, land use history, microclimate.

It is also, precisely, a lookup table for this area. Applied elsewhere, the splits are meaningless: the model can learn "south of northing 5,880,000 is greener" — true here, a mean NDVI of 0.769 against 0.705 to the north (Example 4) — which is a fact about this window rather than about the world.

The measurement shows both effects. In-area map R² rose from 0.348 to 0.509; block CV, which simulates moving, rose much less, 0.258 to 0.283.

There is a middle path: replace raw coordinates with covariates that travel — distance to coast, elevation, aspect, rainfall — which encode the reason a location differs rather than its identity.

Why tree models fail strangely under extrapolation

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

The consequence is peculiar rather than catastrophic-looking: predictions in a new area saturate at a plausible-looking constant. Nothing raises an error, and the map looks flat rather than wrong.

That is visible in the transfer test. With coordinates included, the new-region R² was −0.081 — the model had flattened towards a constant over the new region, because every eastern x value was outside the western training range. Without coordinates, it confidently applied a wrong relationship and scored −0.409. The more "conservative" model was the one that had been flattened.

Why applicability matters more than accuracy

A single accuracy number describes the average over the data you happened to have. A map is used pixel by pixel, and the pixels vary enormously in how well they are supported.

The useful deliverable is therefore two rasters: the prediction, and a mask of where the covariates resemble the training data closely enough for the prediction to mean anything. See Extrapolation in space explained.

Row counts, block sizes, scale and spatial imbalance

• Row count is not sample size. Estimate the effective sample size from the autocorrelation range.
• Coordinates are a deliberate choice, not a default. Include for interpolation, exclude for transfer.
• Tree models cannot extrapolate the target. Predictions saturate at the training range.
• Linear models extrapolate confidently and wrongly. Different failure, same cause.
• Scale matters. A model fitted on 20 m cells does not apply to 1 km cells; the covariates mean different things.
• Class imbalance is usually spatial. Rare classes cluster, so a random split can put an entire class in one fold.
• Blocks must exceed the autocorrelation range, or spatial CV still leaks.
• Report the CV scheme with the score. A bare R² is not interpretable.

Internal links

• Spatial leakage explained: why random cross-validation lies — the measurements behind this page
• How to do spatial cross-validation in Python — the implementation
• Spatial features explained: turning geometry into columns — covariates that travel
• Extrapolation in space explained: the area of applicability — the mask to ship with the map
• Feature importance says coordinates are the best predictor — the coordinate trap in detail
• My model scores 0.95 in testing and fails in the field — the symptom
• Spatial autocorrelation explained — the structure underneath all of this
• How to evaluate a spatial model honestly — what to report instead of one number

FAQ

What makes machine learning on spatial data different?

Nearby observations are correlated, so rows are not independent. That breaks cross-validation, inflates confidence, and makes transfer to new areas unpredictable from internal statistics.

Should I include coordinates as features?

For predicting inside the sampled area, yes — they raised the true map R² from 0.348 to 0.509 here. For a model that must work elsewhere, no: coordinates are a lookup table for one place.

Why does more training data barely help?

Because correlated rows repeat information. Going from 3,000 to 10,000 training cells raised the true map R² only from 0.348 to 0.401.

Can I use a random forest for spatial prediction?

Yes, and it is a sensible default. Be aware it cannot predict outside the range of the training target, so in a new area its predictions saturate rather than error.

How do I know if my model will work in a new region?

Test it on a real held-out region. No cross-validation scheme computed from one region can tell you — block CV promised +0.033 where reality was −0.409.

What is the effective sample size?

Roughly, the number of independent neighbourhoods your samples cover, rather than the number of rows. Estimate it from the autocorrelation range before quoting any confidence interval.

Is spatial machine learning just interpolation with more steps?

For the interpolation task, largely yes, and a good interpolator is often competitive. Machine learning earns its place when covariates carry information that geometry alone cannot.