A species distribution model scores well and predicts badly
On this page
Problem statement
The model reports an AUC of 0.93 on five-fold cross-validation. Projected onto a neighbouring region it puts the species on the wrong side of a mountain range; projected forward in time it produces a map nobody believes.
The score was never a measure of that. Random cross-validation splits rows, and occurrence records are spatially clustered, so nearly every held-out record has a training record within a few kilometres carrying almost identical covariates. The model is being tested on data it has effectively seen.
Quick answer
Split space, not rows, and compare:
import glob
import numpy as np, pandas as pd, geopandas as gpd, rasterio
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import cross_val_score, GroupKFold
from gbif import search
layers = sorted(glob.glob("bio/wc/wc2.1_10m_bio_*.tif"),
key=lambda p: int(p.split("_")[-1].split(".")[0]))
USE = ["bio1", "bio4", "bio5", "bio6", "bio12", "bio15"]
def sample_layers(coords):
out = {}
for path in layers:
name = "bio" + path.split("_")[-1].split(".")[0]
if name not in USE:
continue
with rasterio.open(path) as src:
v = np.array([x[0] for x in src.sample(coords)], dtype="float32")
v[v < -1e30] = np.nan
out[name] = v
return pd.DataFrame(out)[USE]
lynx = search(taxonKey=2435263, hasCoordinate="true").dropna(
subset=["decimalLongitude", "decimalLatitude"])
pres_xy = list(zip(lynx.decimalLongitude, lynx.decimalLatitude))
land = gpd.read_file("bio/ne_land.gpkg")
rng = np.random.default_rng(0)
poly = land.clip([-170, 35, -50, 72]).union_all()
bg = []
while len(bg) < 10_000:
cand = gpd.GeoSeries(gpd.points_from_xy(rng.uniform(-170, -50, 8000),
rng.uniform(35, 72, 8000)), crs=4326)
bg += [(p.x, p.y) for p in cand[cand.within(poly)]]
bg_xy = bg[:10_000]
X = pd.concat([sample_layers(pres_xy), sample_layers(bg_xy)], ignore_index=True)
y = np.r_[np.ones(len(pres_xy)), np.zeros(len(bg_xy))]
xy = np.array(pres_xy + bg_xy)
ok = X.notna().all(axis=1).to_numpy()
X, y, xy = X[ok], y[ok], xy[ok]
block = 5.0 # degrees; the spatial fold size
groups = ((xy[:, 0] // block).astype(int) * 1000 + (xy[:, 1] // block).astype(int))
model = HistGradientBoostingClassifier(random_state=0)
random_auc = cross_val_score(model, X, y, cv=5, scoring="roc_auc")
spatial_auc = cross_val_score(model, X, y, cv=GroupKFold(5), groups=groups,
scoring="roc_auc")
print(f"random 5-fold AUC : {random_auc.mean():.3f} ± {random_auc.std():.3f}")
print(f"spatial 5-fold AUC: {spatial_auc.mean():.3f} ± {spatial_auc.std():.3f}")
print(f"drop: {random_auc.mean() - spatial_auc.mean():.3f}")
print(f"{len(np.unique(groups))} spatial blocks of {block}°")random 5-fold AUC : 0.928 ± 0.027
spatial 5-fold AUC: 0.856 ± 0.038
drop: 0.072
147 spatial blocks of 5.0°
The gap between those two numbers is the part of the score that was memorisation.
Check clustering, block by space, read the spread
1. Check the clustering before trusting any score
Median nearest-neighbour distance among presences, against the resolution of the covariates. If records are closer together than a cell, random folds are meaningless.
2. Block by space
Group records into blocks larger than the autocorrelation range and hold out whole blocks. GroupKFold with a grid-derived group label is enough.
3. Choose the block size deliberately
Too small and neighbours leak between folds; too large and each fold extrapolates to a different climate. Report the size.
4. Look at the fold-to-fold spread, not just the mean
A model that scores 0.85, 0.84, 0.83 and 0.45 across blocks has a region it cannot handle, and the mean hides it.
5. Test the thing you will do
If the model will be projected forward in time or to a new region, evaluate it that way — leave out a region, or fit on one period and test on another.
6. Report both scores
The random score is what most papers print; publishing both makes the difference visible and is the point.
Scores at several block sizes
The data used on this page
Records come from the GBIF API via gbif.search, the paging function from the download guide plus eventDate and countryCode. Climate is WorldClim 2.1 bioclim at 10′ in bio/wc/; pa.protected_areas() is British Columbia's protected-area layer (EPSG:3005). GBIF changes daily, so the setup answers each request with the reply it got on 9 October 2026.
import hashlib
import io
import os
import zipfile
import geopandas as gpd
import requests
# GBIF adds and revises records daily. So that the numbers below reproduce, every GBIF request on
# this page is answered from the reply it got on 9 October 2026 (one file per request, named by a
# SHA-1 of its URL); any other request, or all of them once these lines are deleted, goes live.
SAVED = zipfile.ZipFile(io.BytesIO(requests.get(
"https://data.spatialworkflow.io/gbif-api-2026-10-09.zip", timeout=300).content))
live_get = requests.get
def saved_get(url, params=None, **kwargs):
prepared = requests.Request("GET", url, params=params).prepare()
try:
body = SAVED.read(hashlib.sha1(prepared.url.encode()).hexdigest()[:16] + ".json")
except KeyError:
return live_get(url, params=params, **kwargs)
reply = requests.models.Response()
reply.status_code, reply._content, reply.url = 200, body, prepared.url
reply.headers["Content-Type"] = "application/json"
return reply
requests.get = saved_get
os.makedirs("bio/wc", exist_ok=True)
with open("gbif.py", "w") as f:
f.write('''
import pandas as pd, requests
API = "https://api.gbif.org/v1"
MAX_OFFSET = 100_000
FIELDS = ["key", "species", "speciesKey", "decimalLongitude", "decimalLatitude",
"coordinateUncertaintyInMeters", "year", "month", "basisOfRecord",
"datasetKey", "institutionCode", "issues", "eventDate", "countryCode"]
def count(**params):
return requests.get(f"{API}/occurrence/search", params={**params, "limit": 0}, timeout=120).json()["count"]
def search(limit=None, page=300, **params):
total = count(**params)
if limit is None and total > MAX_OFFSET:
raise ValueError(f"{total:,} records exceeds the {MAX_OFFSET:,} offset limit")
rows, offset = [], 0
target = total if limit is None else min(limit, total)
while offset < target:
body = requests.get(f"{API}/occurrence/search",
params={**params, "limit": min(page, target - offset), "offset": offset},
timeout=180).json()
rows += body["results"]
offset += len(body["results"])
if body["endOfRecords"] or not body["results"]:
break
return pd.DataFrame([{f: r.get(f) for f in FIELDS} for r in rows])
''')
with open("pa.py", "w") as f:
f.write('''
import os
import geopandas as gpd
URL = ("https://services6.arcgis.com/ubm4tcTYICKBpist/arcgis/rest/services/"
"British_Columbia_Parks_Ecological_Reserves_and_Protected_Areas/FeatureServer/0/query"
"?where=1%3D1&outFields=PROTECTED_LANDS_NAME&outSR=3005&f=geojson")
def protected_areas():
if not os.path.exists("bio/bc_protected.gpkg"):
gpd.read_file(URL).to_file("bio/bc_protected.gpkg")
return gpd.read_file("bio/bc_protected.gpkg")
''')
if not os.path.exists("bio/wc/wc2.1_10m_bio_1.tif"):
blob = requests.get("https://geodata.ucdavis.edu/climate/worldclim/2_1/base/wc2.1_10m_bio.zip", timeout=600).content
zipfile.ZipFile(io.BytesIO(blob)).extractall("bio/wc")
if not os.path.exists("bio/ne_land.gpkg"): # Natural Earth 1:110m land, the land mask
blob = requests.get("https://naciscdn.org/naturalearth/110m/physical/ne_110m_land.zip",
timeout=120).content
zipfile.ZipFile(io.BytesIO(blob)).extractall("bio/ne_land")
gpd.read_file("bio/ne_land/ne_110m_land.shp").to_file("bio/ne_land.gpkg", driver="GPKG")
print(len([n for n in os.listdir("bio/wc") if n.endswith(".tif")]), "WorldClim layers")19 WorldClim layers
Example 1 — how clustered the records are
from shapely import STRtree
pres = gpd.GeoSeries(gpd.points_from_xy([p[0] for p in pres_xy],
[p[1] for p in pres_xy]), crs=4326).to_crs(6933)
tree = STRtree(pres.values)
src, nbr = tree.query_nearest(pres.values, exclusive=True, all_matches=False)
nn = gpd.GeoSeries(pres.values[src], crs=6933).distance(
gpd.GeoSeries(pres.values[nbr], crs=6933), align=False) / 1000
print(f"{len(pres):,} presences")
print(f"nearest neighbour: median {np.median(nn):.2f} km, "
f"25th {np.percentile(nn, 25):.2f} km, 75th {np.percentile(nn, 75):.2f} km")
print(f"records closer together than one 10-minute cell (18.5 km): "
f"{(nn < 18.5).mean():.0%}")6,353 presences
nearest neighbour: median 5.89 km, 25th 1.87 km, 75th 16.66 km
records closer together than one 10-minute cell (18.5 km): 79%
Example 2 — the score at several block sizes
for deg in (1.0, 2.5, 5.0, 10.0):
g = ((xy[:, 0] // deg).astype(int) * 10_000 + (xy[:, 1] // deg).astype(int))
n_blocks = len(np.unique(g))
if n_blocks < 5:
continue
auc = cross_val_score(model, X, y, cv=GroupKFold(5), groups=g, scoring="roc_auc")
print(f"{deg:>5}° blocks ({n_blocks:>4} blocks): AUC {auc.mean():.3f} "
f"(folds {', '.join(f'{a:.2f}' for a in auc)})") 1.0° blocks (2560 blocks): AUC 0.895 (folds 0.92, 0.91, 0.82, 0.92, 0.90)
2.5° blocks ( 483 blocks): AUC 0.891 (folds 0.93, 0.91, 0.88, 0.88, 0.86)
5.0° blocks ( 147 blocks): AUC 0.856 (folds 0.92, 0.80, 0.86, 0.84, 0.85)
10.0° blocks ( 52 blocks): AUC 0.822 (folds 0.76, 0.86, 0.80, 0.91, 0.78)
Example 3 — leave one region out
regions = {"west of 115°W": xy[:, 0] < -115,
"115°W to 90°W": (xy[:, 0] >= -115) & (xy[:, 0] < -90),
"east of 90°W": xy[:, 0] >= -90}
from sklearn.metrics import roc_auc_score
for name, mask in regions.items():
if mask.sum() < 200 or (~mask).sum() < 200 or y[mask].sum() < 20:
continue
fit = HistGradientBoostingClassifier(random_state=0).fit(X[~mask], y[~mask])
score = roc_auc_score(y[mask], fit.predict_proba(X[mask])[:, 1])
print(f"trained without {name:<16} tested on it: AUC {score:.3f} "
f"({int(y[mask].sum()):,} presences held out)")trained without west of 115°W tested on it: AUC 0.652 (4,483 presences held out)
trained without 115°W to 90°W tested on it: AUC 0.846 (516 presences held out)
trained without east of 90°W tested on it: AUC 0.843 (1,344 presences held out)
Leaving out a whole region is the closest cheap approximation to what a projection actually does, and it scores lower than any of the fold-based numbers. That is not the model failing; it is the evaluation finally asking the right question.
Why autocorrelation inflates scores
Why spatial autocorrelation inflates the score
Two records 5 km apart share their covariate values almost exactly, because the covariates are a smooth 18 km grid. Put one in training and one in testing and the model can score the test record correctly by having memorised its twin. With thousands of clustered records, most of the test set has a twin, so the score approaches the model's ability to interpolate rather than to predict.
Why the block size is a judgement
The right block is larger than the distance over which records and covariates are autocorrelated. That distance is not a single number — it varies by variable and by region — so the practical approach is to compute the score at several block sizes and show the curve. It falls as blocks grow and flattens once the blocks exceed the autocorrelation range.
Why a lower score is a better model
Nothing about the model changes between the two evaluations. Only the question changes, from "can you fill in a gap among records you have seen" to "can you say something about somewhere you have not". The second question is the one every application of a distribution model asks, so the second number is the one that belongs in the abstract.
Buffered leave-one-out and environmental blocks
- Buffered leave-one-out removes a disc around each test record instead of blocking; it is stricter and slower.
- Environmental blocking splits by covariate space rather than geography, which tests extrapolation directly.
- blockCV (R) and
spacv(Python) implement several of these schemes with autocorrelation-based block sizing. - Temporal splits matter for projections: fit on older records, test on newer.
- AUC itself is background-dependent in presence-only work; compare models on the same background or use another metric.
- A very high blocked score usually means the background is too wide, not that the model is excellent.
Internal links
- How to fit a species distribution model in Python — the model being evaluated
- Species distribution models explained — why the background matters more
- How to thin occurrence records spatially — reducing clustering before fitting
- Sampling bias in occurrence data explained — the source of the clustering
- How to cross-validate an interpolated surface — the same argument for interpolation
- Spatial autocorrelation explained — why neighbours are not independent
FAQ
Why is my SDM's AUC so high?
Because random folds put clustered records in both training and testing. Blocked folds score lower.
What block size should I use?
Larger than the autocorrelation range of the records and covariates; compute the score at several sizes and report which you used.
Is a lower score a worse model?
No. It is a more honest measure of the same model on a harder, more relevant question.
Does thinning fix this?
It reduces the clustering, which narrows the gap — but does not remove the need for spatial evaluation.
What if one fold scores far worse than the others?
That fold is a region the model cannot handle. Report it rather than the mean.
Should I still report the random-fold score?
Report both. The difference between them is the finding.