How to fit a species distribution model in Python

On this page

Problem statement

A presence-only distribution model has three moving parts: the presences, a background sample that defines "available", and a classifier. The classifier is the least important of the three — swapping logistic regression for gradient boosting changes the map less than moving the background from "the continent" to "the region the species could reach".

This guide fits a lynx model on WorldClim covariates with scikit-learn, evaluates it both ways, and shows how much of the score is an artefact of where the background came from.

Quick answer

Presences, background, covariates, classifier:

import glob
import numpy as np, pandas as pd, geopandas as gpd, rasterio
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
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"]     # screened for collinearity

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")
pres = lynx.dropna(subset=["decimalLongitude", "decimalLatitude"])
pres_xy = list(zip(pres.decimalLongitude, pres.decimalLatitude))

land = gpd.read_file("bio/ne_land.gpkg")
minx, miny, maxx, maxy = (-170, 35, -50, 72)                 # the species' accessible area
rng = np.random.default_rng(0)
bg = []
poly = land.clip([minx, miny, maxx, maxy]).union_all()
while len(bg) < 10_000:
    cand = gpd.GeoSeries(gpd.points_from_xy(rng.uniform(minx, maxx, 8000),
                                            rng.uniform(miny, maxy, 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))]
ok = X.notna().all(axis=1).to_numpy()
X, y = X[ok], y[ok]
print(f"{int(y.sum()):,} presences, {int((y == 0).sum()):,} background points, "
      f"{X.shape[1]} covariates")

from sklearn.preprocessing import PolynomialFeatures

models = {
    "logistic, linear terms": make_pipeline(StandardScaler(),
                                            LogisticRegression(max_iter=2000)),
    "logistic, + quadratics": make_pipeline(
        StandardScaler(), PolynomialFeatures(2, include_bias=False),
        LogisticRegression(max_iter=5000)),
    "gradient boosting": HistGradientBoostingClassifier(random_state=0),
}
for name, model in models.items():
    auc = cross_val_score(model, X, y, cv=5, scoring="roc_auc")
    print(f"  {name:<24} random 5-fold AUC {auc.mean():.3f} ± {auc.std():.3f}")

6,343 presences, 9,978 background points, 6 covariates
  logistic, linear terms   random 5-fold AUC 0.542 ± 0.050
  logistic, + quadratics   random 5-fold AUC 0.869 ± 0.041
  gradient boosting        random 5-fold AUC 0.928 ± 0.027

A linear logistic regression scores 0.54 — barely better than a coin toss. Adding quadratic terms to the same classifier lifts it to 0.87, and gradient boosting reaches 0.93. The difference is not the algorithm's cleverness: it is that a species occupies a middle range of climate, and a model with only linear terms can express "colder is always better" but not "cold, but not the coldest".

Three inputs; the background decides what the model compares, and the feature shape decides what it can say.

Accessible area, screened covariates, background on land

1. Define the accessible area first

The background sample answers "available to this species". A continental background makes a regional species look perfectly predicted; a background restricted to its reachable range makes the model work for its living.

2. Screen the covariates

Six to eight weakly correlated variables chosen on ecological grounds beat all nineteen bioclim layers.

3. Sample background on land, in the right numbers

Ten thousand background points is the usual convention. Fewer makes the class boundary noisy; many more mostly costs time.

4. Let the model express an optimum

A species occupies a middle range of temperature, not "as much as possible". A logistic regression in raw covariates cannot say that; add quadratic terms, splines, or use a learner that finds the shape itself.

5. Evaluate spatially, not randomly

Random folds put a presence and its neighbour in different folds. The score that results measures interpolation, not prediction.

6. Predict onto the grid and look at it

A map with a colour ramp from zero to one is the output most readers will judge; check that it is not extrapolating outside the covariate range it was fitted on.

The same presences, two definitions of "available", two different models.

How much the background choice matters

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 much the background choice matters

def build(minx, miny, maxx, maxy, n=10_000, seed=0):
    rng = np.random.default_rng(seed)
    poly = land.clip([minx, miny, maxx, maxy]).union_all()
    pts = []
    while len(pts) < n:
        cand = gpd.GeoSeries(gpd.points_from_xy(rng.uniform(minx, maxx, 8000),
                                                rng.uniform(miny, maxy, 8000)), crs=4326)
        pts += [(p.x, p.y) for p in cand[cand.within(poly)]]
    Xb = pd.concat([sample_layers(pres_xy), sample_layers(pts[:n])], ignore_index=True)
    yb = np.r_[np.ones(len(pres_xy)), np.zeros(n)]
    good = Xb.notna().all(axis=1).to_numpy()
    return Xb[good], yb[good]

extents = {"whole world": (-180, -60, 180, 84),
           "North America": (-170, 35, -50, 72),
           "Canada and the northern US": (-141, 42, -52, 70)}
model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000))
for label, box in extents.items():
    Xb, yb = build(*box)
    auc = cross_val_score(model, Xb, yb, cv=5, scoring="roc_auc")
    print(f"{label:<28} AUC {auc.mean():.3f}")

whole world                  AUC 0.771
North America                AUC 0.542
Canada and the northern US   AUC 0.523

The AUC falls as the background narrows, and the narrow model is the useful one. A continent-wide background asks "is this lynx habitat or is it Mexico", which any model can answer.

Example 2 — what the model learned

fitted = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)).fit(X, y)
coef = pd.Series(fitted[-1].coef_[0], index=X.columns).sort_values(key=abs, ascending=False)
print("standardised coefficients:")
for name, c in coef.items():
    print(f"  {name:>6} {c:>7.2f}  {'+' if c > 0 else '−'}")

pres_env, bg_env = X[y == 1], X[y == 0]
print(f"\n{'variable':>8}{'presences':>12}{'background':>12}{'difference':>12}")
for c in X.columns:
    print(f"{c:>8}{pres_env[c].mean():>12.1f}{bg_env[c].mean():>12.1f}"
          f"{pres_env[c].mean() - bg_env[c].mean():>12.1f}")

standardised coefficients:
    bio6   -1.77  −
    bio1    1.17  +
    bio4   -0.72  −
    bio5   -0.08  −
   bio15    0.06  +
   bio12   -0.03  −

variable   presences  background  difference
    bio1        -1.2        -0.9        -0.3
    bio4      1223.5      1187.4        36.1
    bio5        21.3        21.1         0.2
    bio6       -23.0       -21.8        -1.2
   bio12       551.1       574.0       -22.9
   bio15        47.7        46.9         0.8

The differences between presences and background are tiny — a third of a degree in annual mean temperature, a millimetre or two of rain — which is why the linear model can barely separate them. The coefficients are still readable: the strongest is a negative weight on the minimum temperature of the coldest month, which for a linear term means "colder winters, more lynx". That is true across the range this background covers and would stop being true if the background included the high Arctic, and reading it as ecology rather than as a description of this particular comparison is the standard way to over-interpret a distribution model.

Example 3 — predicting onto a grid

with rasterio.open(layers[0]) as src:
    window = src.window(-141, 42, -52, 70)
    transform = src.window_transform(window)
    shape = (int(window.height), int(window.width))

stack = {}
for path in layers:
    name = "bio" + path.split("_")[-1].split(".")[0]
    if name not in USE:
        continue
    with rasterio.open(path) as src:
        band = src.read(1, window=src.window(-141, 42, -52, 70)).astype("float32")
        band[band < -1e30] = np.nan
    stack[name] = band.ravel()
grid = pd.DataFrame(stack)[USE]
valid = grid.notna().all(axis=1).to_numpy()

pred = np.full(len(grid), np.nan, dtype="float32")
pred[valid] = fitted.predict_proba(grid[valid])[:, 1]
surface = pred.reshape(shape)

print(f"prediction grid {shape}, {valid.sum():,} land cells")
print(f"suitability: median {np.nanmedian(surface):.3f}, "
      f"90th percentile {np.nanpercentile(surface, 90):.3f}")
for cut in (0.5, 0.7, 0.9):
    print(f"  cells above {cut}: {np.nansum(surface > cut):>7,} "
          f"({np.nansum(surface > cut) / valid.sum():.1%} of land)")

prediction grid (168, 534), 64,971 land cells
suitability: median 0.397, 90th percentile 0.501
  cells above 0.5:   6,722 (10.3% of land)
  cells above 0.7:       0 (0.0% of land)
  cells above 0.9:       0 (0.0% of land)

The background is the model

Why the background is the model

A presence-background model estimates the ratio between the environment where the species was recorded and the environment that was available. "Available" is not in the data — it is asserted by the background sample, and the assertion carries all the ecology: which region the species could reach, over what time, and through what barriers. Two analysts with the same presences and different accessible areas will produce different maps and both be doing it correctly.

Why AUC on random folds is inflated

Occurrence records are spatially clustered, so a random split puts records from the same cluster in both the training and the test fold. The model can then get the test records right by remembering their neighbours, and the score measures how well it interpolates within clusters. Spatially blocked folds are the standard remedy and always score lower.

Why the response shape matters more than the algorithm

The gap between 0.54 and 0.93 above is not a gap between algorithms, it is a gap between a model that can express an environmental optimum and one that cannot. Lynx occupy a band of winter temperature, not its extreme; a linear term can only say "colder is better" or "warmer is better", and neither is true. Adding quadratic terms to the same logistic regression recovers most of the difference — which is why MaxEnt, whose feature classes include hinges and quadratics by default, has been the field standard for twenty years.

The same presences scored against three definitions of available: the widest one looks best.

MaxEnt, prevalence and extrapolation

MaxEnt is a penalised log-linear presence-background model; its outputs are comparable to these, and elapid brings it to Python.
Prevalence is unidentifiable from presence-only data, so predicted values are relative suitability, not probability of occurrence.
Target-group backgrounds use the records of similar species instead of random points, which carries the sampling bias into the background.
Extrapolation: check whether prediction cells fall outside the covariate range of the training data; MESS maps do this.
Thresholding a continuous surface into presence/absence needs a rule (maximum sensitivity+specificity, 10th percentile of presences), and the rule changes the area.
Future projections need the same variables from a climate model, and the extrapolation question becomes central.

Internal links

Species distribution models explained — the assumptions behind this code
How to extract climate covariates at occurrence points — the inputs
A species distribution model scores well and predicts badly — the evaluation
How to thin occurrence records spatially — the presences this expects
Sampling bias in occurrence data explained — what a target-group background corrects
How to run a weighted site suitability analysis in Python — the same map from expert weights instead

FAQ

What is a presence-background model?

A classifier fitted to distinguish the environment at occurrence records from the environment at a sample of available locations.

How many background points do I need?

Ten thousand is the usual convention; the class boundary gets noisy below a few thousand.

Which algorithm should I use?

It matters less than the background. Logistic regression is interpretable, gradient boosting is flexible, MaxEnt is the field standard.

Is the output a probability of occurrence?

No — it is relative suitability. Prevalence cannot be identified from presence-only data.

Why does my AUC drop when I restrict the background?

Because the model is being asked a harder, more useful question. The lower score is the more honest one.

Do I need to thin the presences first?

Usually yes, to keep clustered records from dominating the fit.