Feature Importance Says Coordinates Are the Best Predictor
On this page
Problem statement
You add x and y to a model, the score improves, and the importance table puts them near the top. The physical story you were building — slope drives vegetation, elevation drives temperature — quietly changes.
Measured on a real problem, predicting Sentinel-2 NDVI on Snowdon from terrain with a random forest on 3,000 cells (importances measured here):
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
Coordinates took second and fourth place and pushed TPI from third to sixth. And the model genuinely improved: the true map R² rose from 0.348 to 0.509 (Example 1).
Both things are true at once. The coordinates are informative here and meaningless anywhere else.
Quick answer
Decide from the deployment, not from the importance table:
# predicting inside the sampled area: coordinates help
X = np.column_stack([terrain, coords]) # true map R2 0.509
# predicting anywhere else: they are a lookup table
X = terrain # true map R2 0.348The same features under block cross-validation, which simulates moving:
terrain only 0.258
terrain + x,y 0.283
A much smaller gain, because a coordinate learned in one place transfers poorly to another.
What coordinates stand in for and what they cost
1. Understand what the coordinates are standing in for
A tree given x and y can split the study area into regions and fit a different local mean in each. That captures everything your measured covariates miss: soil, land-use history, microclimate, management.
It is a genuine improvement in predictive accuracy for the sampled area — measured, 0.348 to 0.509 — and it is not a physical explanation.
2. Notice that importance is relative
Adding coordinates did not make slope less important in any physical sense. It made slope's share smaller, because importance sums to one.
TPI fell from 0.127 to 0.077, which means the coordinates absorbed a good part of what TPI was explaining — the two carry overlapping information about local position.
Importance is a statement about this fitted model on this data, not about the world.
3. Check the transfer cost before keeping them
in-area R² block CV R²
terrain only 0.348 0.258
terrain + x,y 0.509 0.283
The in-area gain is 0.161; the block-CV gain is 0.025. The difference is the part that does not travel.
If the model will only ever predict inside the sampled area, keep them. If it must move, the honest comparison is the block-CV column.
4. Prefer covariates that travel
Replace raw coordinates with variables that encode why a location differs:
- distance to the coast, to a river, to a city
- elevation, slope, aspect, topographic position
- long-term climate normals
- geology or soil class
These carry the same information where it is causal, and they mean the same thing in a new area. A model using them can transfer; one using x and y cannot.
Test the replacement the same way, because a smooth covariate over a small area is a coordinate in disguise. Across this 10 km window, distance to the coast runs from 6.9 to 17.5 km in one direction, and adding it behaved much like adding x: the in-area R² rose from 0.348 to 0.420, and block CV fell slightly, from 0.258 to 0.249 (Example 3). It will only travel to places whose coast does the same thing.
5. Use permutation importance, on a held-out set
Tree impurity importance is biased towards high-cardinality continuous features — which coordinates are, maximally. Permutation importance on held-out data measures what the model actually uses to predict data it has not seen. Here, on 20,000 unsampled cells, it agreed with impurity and gave the coordinates slightly more, 26% of the total against 23% (Example 2): the model really does lean on location here, it is not an artefact of the importance measure.
Permutation importance and covariates that travel
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_grid and y_grid hold every remaining cell — the map to predict — and X, y and coords a random 3,000 of them, the training sample.
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"]
X_grid, y_grid, coords_grid = X, y, coords
sample = np.random.default_rng(0).choice(len(y_grid), 3000, replace=False)
X, y, coords = X_grid[sample], y_grid[sample], coords_grid[sample]
print(f"{item.id}: {shape[1]} x {shape[0]} cells of 20 m, {valid.sum():,} usable "
f"({valid.mean():.1%}); training on {len(y):,}")S2C_30UVD_20250712_0_L2A: 543 x 509 cells of 20 m, 270,536 usable (97.9%); training on 3,000
Example 1 — measuring what coordinates cost and buy
import numpy as np
from sklearn.cluster import KMeans
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from sklearn.model_selection import GroupKFold, KFold
def coordinate_effect(X, y, coords, X_holdout=None, y_holdout=None,
coords_holdout=None, n_blocks=10):
"""In-area accuracy against transferability, with and without x, y."""
def model():
return RandomForestRegressor(n_estimators=200, min_samples_leaf=2,
n_jobs=-1, random_state=0)
def cv(features, groups=None):
splitter = (GroupKFold(5).split(features, y, groups) if groups is not None
else KFold(5, shuffle=True, random_state=0).split(features))
predicted = np.empty(len(y))
for train, test in splitter:
predicted[test] = model().fit(features[train], y[train]) \
.predict(features[test])
return r2_score(y, predicted)
blocks = KMeans(n_blocks, random_state=0, n_init=10).fit_predict(coords)
with_xy = np.column_stack([X, coords])
rows = []
for label, features in (("terrain only", X), ("terrain + x,y", with_xy)):
random = cv(features)
spatial = cv(features, groups=blocks)
true = None
if X_holdout is not None:
holdout = (np.column_stack([X_holdout, coords_holdout])
if label.endswith("x,y") else X_holdout)
true = r2_score(y_holdout, model().fit(features, y).predict(holdout))
rows.append((label, random, spatial, true))
print(f" {label:16} random {random:+.3f} block {spatial:+.3f}"
+ (f" true {true:+.3f}" if true is not None else ""))
gain_in = rows[1][1] - rows[0][1]
gain_out = rows[1][2] - rows[0][2]
print(f" coordinates gain {gain_in:+.3f} in-area, {gain_out:+.3f} "
"under block CV — the difference is what does not travel")
return rows
rest = np.setdiff1d(np.arange(len(y_grid)), sample)
rows = coordinate_effect(X, y, coords, X_grid[rest], y_grid[rest], coords_grid[rest]) terrain only random +0.374 block +0.258 true +0.348
terrain + x,y random +0.504 block +0.283 true +0.509
coordinates gain +0.130 in-area, +0.025 under block CV — the difference is what does not travel
Example 2 — permutation importance, which is harder to fool
import numpy as np
from sklearn.inspection import permutation_importance
def honest_importance(model, X_test, y_test, feature_names, n_repeats=10):
"""Permutation importance on held-out data, with the spread shown."""
result = permutation_importance(model, X_test, y_test,
n_repeats=n_repeats, random_state=0,
n_jobs=-1)
order = np.argsort(-result.importances_mean)
print(f" permutation importance over {n_repeats} repeats:")
for i in order:
mean = result.importances_mean[i]
std = result.importances_std[i]
bar = "#" * int(max(mean, 0) * 60)
print(f" {feature_names[i]:14} {mean:+.4f} ± {std:.4f} {bar}")
coordinate_share = sum(result.importances_mean[i]
for i, name in enumerate(feature_names)
if name in ("x", "y", "lon", "lat", "easting",
"northing"))
total = max(result.importances_mean.sum(), 1e-12)
if coordinate_share / total > 0.2:
print(f" ! coordinates account for {coordinate_share / total:.0%} of "
"the importance — the model will not transfer")
return result
names = FEATURES + ["x", "y"]
with_xy = RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1,
random_state=0).fit(np.column_stack([X, coords]), y)
held = np.random.default_rng(4).choice(rest, 20_000, replace=False)
print(" impurity:", ", ".join(f"{n} {v:.3f}" for n, v in sorted(zip(names, with_xy.feature_importances_),
key=lambda t: -t[1])))
honest_importance(with_xy, np.column_stack([X_grid[held], coords_grid[held]]), y_grid[held], names) impurity: elevation 0.443, y 0.118, slope 0.113, x 0.109, northness 0.080, tpi 0.077, eastness 0.059
permutation importance over 10 repeats:
elevation +0.5317 ± 0.0040 ###############################
y +0.1572 ± 0.0035 #########
slope +0.1203 ± 0.0029 #######
x +0.1127 ± 0.0019 ######
northness +0.0473 ± 0.0024 ##
tpi +0.0376 ± 0.0018 ##
eastness +0.0189 ± 0.0011 #
! coordinates account for 26% of the importance — the model will not transfer
Permutation importance measures the drop in score when a feature is shuffled, on data the model has not seen. Impurity importance measures how often a feature was split on during training, which favours continuous high-cardinality features regardless of usefulness.
Example 3 — replacing coordinates with covariates that travel
import numpy as np
from scipy.spatial import cKDTree
def positional_covariates(coords, coast=None, rivers=None, cities=None,
elevation=None):
"""Variables that encode why a location differs, not which one it is."""
features, names = [], []
for label, targets in (("dist_coast", coast), ("dist_river", rivers),
("dist_city", cities)):
if targets is None:
continue
distances, _ = cKDTree(np.asarray(targets)).query(coords, k=1)
features.append(distances)
names.append(label)
if elevation is not None:
features.append(np.asarray(elevation))
names.append("elevation")
if not features:
raise ValueError("no positional covariates supplied")
print(f" {len(names)} transferable covariates: {names}")
print(" these mean the same thing in a new area; x and y do not")
return np.column_stack(features), names
# Natural Earth's 1:10m coastline around Wales, densified to 100 m, in the grid's CRS
import geopandas as gpd
from shapely.geometry import box
coast = gpd.read_file("https://naciscdn.org/naturalearth/10m/physical/ne_10m_coastline.zip",
bbox=(-6, 52, -2, 55)).clip(box(-6, 52, -2, 55)).to_crs(crs)
vertices = np.vstack([np.asarray(g.coords) for g in coast.geometry.segmentize(100).explode().geometry])
dist, _ = positional_covariates(coords_grid, coast=vertices)
print(f" distance to the coast over the window: {dist.min() / 1000:.1f} to {dist.max() / 1000:.1f} km")
coordinate_effect(np.column_stack([X, dist[sample]]), y, coords,
np.column_stack([X_grid[rest], dist[rest]]), y_grid[rest], coords_grid[rest]) 1 transferable covariates: ['dist_coast']
these mean the same thing in a new area; x and y do not
distance to the coast over the window: 6.9 to 17.5 km
terrain only random +0.439 block +0.249 true +0.420
terrain + x,y random +0.530 block +0.293 true +0.536
coordinates gain +0.092 in-area, +0.043 under block CV — the difference is what does not travel
Distance to the coast is 40 km in one place and 40 km in another, and in both it implies a similar maritime influence. An easting of 422,000 implies nothing outside its own UTM zone.
Why x and y win inside a study area
Why coordinates are so effective within a study area
Spatial data is autocorrelated: nearby places have similar values. A model that knows where a point is can exploit that directly, without needing any covariate to explain it.
For a tree, x and y are ideal split variables — continuous, high-cardinality, and able to isolate any region. Given enough depth, a tree with coordinates approximates a nearest-neighbour interpolator.
That is why the in-area gain is real and large: 0.348 to 0.509 in the measurement here. It is also why the model has stopped being a model of the process and become a model of the place.
Why impurity importance overstates them
Tree impurity importance counts how much each feature reduced impurity across the splits it was used for. Features with many possible split points get more opportunities.
Continuous coordinates have as many split points as there are distinct values — the maximum. Categorical features with few levels have the fewest.
So impurity importance is biased towards coordinates by construction, before any question of usefulness. Permutation importance on held-out data avoids this, at the cost of being slower. When the two agree, as they did here, the model's dependence on location is genuine.
Why the importance table is not a physical explanation
Three reasons, all of which apply here.
It is relative. Adding a feature reduces every other share, so TPI falling from 0.127 to 0.077 says nothing about TPI's physical relevance.
It is model-specific. A different model on the same data gives a different ranking.
Correlated features share arbitrarily. Coordinates and TPI both encode local position, and the split between them depends on the fitting, not on the world.
Use importance to understand the model. Use domain reasoning and controlled experiments to understand the process.
Why replacing coordinates is worth the effort
A model with coordinates is a very good interpolator for one area. A model with transferable covariates is a weaker interpolator that can be applied elsewhere, checked against new data, and criticised on physical grounds.
The measured cost of dropping coordinates is 0.161 R² in-area. Under block CV the gap is only 0.025 — and in a genuine new region, coordinates are worth nothing at all.
Which trade to make depends entirely on what the model is for, and that is the question the importance table cannot answer.
Impurity bias, correlated features and latitude
- Importance is relative and sums to one. Adding a feature reduces every other share.
- Impurity importance favours continuous high-cardinality features — coordinates maximally.
- Use permutation importance on held-out data.
- Coordinates help in-area and not elsewhere — 0.161 against 0.025 under block CV here.
- Prefer distance-to-feature covariates that mean the same thing anywhere.
- Correlated features split importance arbitrarily.
- A latitude term is not the same as
y— it has a physical meaning through insolation. - Report the feature list with the model; two models with different features are different models.
Internal links
- Spatial machine learning explained — the wider set of consequences
- Spatial leakage explained: why random cross-validation lies — the measurements behind this
- Spatial features explained: turning geometry into columns — covariates that travel
- My model scores 0.95 in testing and fails in the field — where this leads
- How to do spatial cross-validation in Python — measuring the transfer cost
- How to build spatial features for a machine learning model — constructing the alternatives
- Extrapolation in space explained: the area of applicability — what to ship instead
- Spatial autocorrelation explained — why coordinates work at all
FAQ
Why do coordinates dominate my feature importance?
Because they are continuous and high-cardinality, which impurity importance favours, and because spatial autocorrelation makes location genuinely predictive within the study area.
Does that mean my model is wrong?
No. Measured, adding coordinates raised the true in-area R² from 0.348 to 0.509. It means the model has learned this place rather than the process.
Should I remove the coordinates?
If the model must work elsewhere, yes. If it will only ever predict inside the sampled area, keeping them is a legitimate accuracy gain.
What should I use instead?
Covariates that encode why a location differs: distance to coast, river or city, elevation, aspect, climate normals. These mean the same thing in a new area — but check them under block CV too, because over a small window a smooth covariate behaves like a coordinate.
Why did another feature's importance drop when I added coordinates?
Importance sums to one, so shares redistribute. TPI fell from 0.127 to 0.077 here without becoming less physically relevant.
Is permutation importance better?
Yes, for this purpose. It measures the score drop when a feature is shuffled on held-out data, rather than counting splits during training. Here it agreed with impurity importance, which is itself informative: the dependence on location was real.
Can I trust feature importance as an explanation?
Only as an explanation of the model. It is relative, model-specific, and splits correlated features arbitrarily.