SDM Projection and Transferability in Python

A species distribution model fitted on today’s climate and evaluated by cross-validation tells you how well it interpolates. Projecting it to 2070, or to a mountain range it was never trained on, asks something entirely different: whether the fitted relationships hold outside the conditions that produced them. That question is about transferability, and it is the part of species distribution modelling where confident maps most often outrun their evidence. This page covers doing it honestly — projecting to new conditions, measuring how far outside the training envelope each pixel sits, and reporting the uncertainty alongside the map. It belongs to Species Distribution Modeling with MaxEnt and consumes the model produced by MaxEnt Model Training & Tuning.

The concrete problem: you have a tuned MaxEnt model for a forest tree or an indicator species, a stack of future climate scenarios, and a request for a map of where the species will be suitable in fifty years. Producing that map is easy. Producing one whose caveats are legible in the map itself is the work.

Prerequisites

• A tuned, validated model with recorded feature classes and regularisation, from Tuning MaxEnt Regularization with betamultiplier Grid Search.
• A training predictor stack and one or more projection stacks with identical variable names, order, units and grid.
• Python ≥ 3.10 with rasterio, numpy, geopandas, scikit-learn, and either elapid or a MaxEnt binary wrapper.
• At least three climate models per scenario if the projection is to future climate — a single GCM is one opinion, not a projection.
• A spatially blocked validation already run, from Spatial Block Cross-Validation for SDM in Python — random cross-validation over-states transferability badly.
• A decision on clamping, made and recorded before the first projection is run.

Identical variable order is worth checking mechanically. Most projection errors are not subtle statistical failures; they are a stack whose bands were written in a different sequence, so the model reads mean diurnal range where it expects annual precipitation and produces a map that is entirely wrong and entirely plausible.

Concept: Novelty, Clamping and What They Cost

A model fitted between 4 °C and 18 °C mean annual temperature has learned nothing about 22 °C. Two things can happen when it meets that value.

Unclamped, the fitted response continues along whatever trajectory the last hinge or quadratic implies. For a quadratic that means the response turns over and can fall to zero; for a hinge it means unbounded extrapolation. Neither is evidence.

Clamped, values beyond the training range are held at the boundary value. The response is flat rather than fabricated, which is more conservative but hides where the clamping happened.

$$\mathrm{MESS}(x)=\underset{j}{\min }100\times \{\begin{array}{ll}\frac{x_j-{\min }_j}{{\max }_j-{\min }_j}& x_j<{\min }_j\\ \text{percentile-based interior score}& {\min }_j\le x_j\le {\max }_j\\ \frac{{\max }_j-x_j}{{\max }_j-{\min }_j}& x_j>{\max }_j\end{array}$$

The multivariate environmental similarity surface reduces that to one number per pixel: negative means at least one predictor is outside its training range, and the magnitude says how far. Reporting it alongside the projection is the difference between “suitable” and “suitable, on an extrapolation the data cannot support”.

Step-by-Step Python Pipeline

Step 1 — Verify the projection stack matches the training stack

import numpy as np
import rasterio


def assert_stacks_compatible(train_paths: list[str], proj_paths: list[str]) -> dict:
    """Refuse to project unless variable order, grid and CRS all agree."""
    assert len(train_paths) == len(proj_paths), "different number of predictors"
    info = []
    for tp, pp in zip(train_paths, proj_paths):
        with rasterio.open(tp) as t, rasterio.open(pp) as p:
            assert t.crs == p.crs, f"CRS differs: {tp} vs {pp}"
            assert abs(abs(t.transform.a) - abs(p.transform.a)) < 1e-6, \
                f"resolution differs: {tp} vs {pp}"
            info.append({
                "variable": tp.split("/")[-1],
                "train_range": (float(t.read(1, masked=True).min()),
                                float(t.read(1, masked=True).max())),
                "proj_range": (float(p.read(1, masked=True).min()),
                               float(p.read(1, masked=True).max())),
            })
    return {"n_vars": len(info), "variables": info}

Print that table before every projection. It is the cheapest possible check and it catches the reordered-band error, the units error (precipitation in mm against cm), and the scaling error (temperature ×10, as several climate products ship it) all at once.

Step 2 — Project, with and without clamping

import numpy as np


def project(model, stack: np.ndarray, train_min: np.ndarray, train_max: np.ndarray,
            clamp: bool = True) -> np.ndarray:
    """Apply a fitted model to a predictor stack of shape (n_vars, y, x)."""
    n_vars, ny, nx = stack.shape
    flat = stack.reshape(n_vars, -1).T                       # (n_pixels, n_vars)
    valid = np.isfinite(flat).all(axis=1)

    x = flat[valid]
    if clamp:
        x = np.clip(x, train_min, train_max)

    out = np.full(flat.shape[0], np.nan, dtype="float32")
    out[valid] = model.predict(x).astype("float32")
    return out.reshape(ny, nx)

Produce both. The difference between the clamped and unclamped surfaces is itself a diagnostic: where they agree, extrapolation is not driving the answer; where they diverge by more than a few hundredths of suitability, the map is being written by the model’s behaviour outside its data.

Step 3 — Compute the extrapolation surface

The MESS calculation is short, and writing it out is worth it because the interior score is where implementations differ.

import numpy as np


def mess(train: np.ndarray, proj: np.ndarray) -> np.ndarray:
    """Multivariate environmental similarity surface.

    train: (n_train, n_vars) reference values; proj: (n_pixels, n_vars).
    Negative output = at least one variable outside the training range.
    """
    n_train = train.shape[0]
    sims = np.empty_like(proj, dtype="float64")
    for j in range(train.shape[1]):
        ref = np.sort(train[:, j])
        lo, hi = ref[0], ref[-1]
        rng = max(hi - lo, 1e-12)
        v = proj[:, j]
        # Fraction of reference values below v, as a percentage.
        f = np.searchsorted(ref, v, side="right") / n_train * 100.0
        s = np.where(f == 0, (v - lo) / rng * 100.0,
            np.where(f <= 50, 2.0 * f,
            np.where(f < 100, 2.0 * (100.0 - f),
                     (hi - v) / rng * 100.0)))
        sims[:, j] = s
    return sims.min(axis=1)

The minimum across variables is deliberate: a pixel is novel if any predictor is outside its range, because the model has no information about that combination regardless of how ordinary the other variables look.

Step 4 — Ensemble across climate models

A single general circulation model produces a confident, precise, and largely arbitrary map. Running several and reporting the spread is the minimum defensible practice.

import numpy as np


def ensemble(projections: list[np.ndarray]) -> dict:
    """Median suitability and agreement across an ensemble of projections."""
    arr = np.stack(projections)                      # (n_models, y, x)
    median = np.nanmedian(arr, axis=0)
    iqr = (np.nanpercentile(arr, 75, axis=0) - np.nanpercentile(arr, 25, axis=0))
    # Agreement: fraction of models placing the pixel above the binary threshold.
    return {"median": median.astype("float32"),
            "iqr": iqr.astype("float32"),
            "n_models": len(projections)}


def agreement(projections: list[np.ndarray], threshold: float) -> np.ndarray:
    """Fraction of ensemble members calling a pixel suitable."""
    arr = np.stack(projections)
    return (arr >= threshold).mean(axis=0).astype("float32")

Publish the agreement layer with the median. A pixel where eight of nine models agree is a different claim from one where five do, and a median map alone erases that distinction entirely.

Validation & Verification

Transferability cannot be validated by the data that trained the model, which is the whole difficulty. Three approaches give real evidence, in descending order of strength.

Independent data from the projection domain. If any occurrence records exist in the target region or period, evaluate against them directly. This is the only true test, and even a small sample is worth more than any internal statistic.

Spatially blocked cross-validation. Holding out geographic blocks approximates transfer to a new region. AUC from blocked folds is typically 0.05–0.15 lower than from random folds, and that gap is roughly the amount by which random validation over-states transferability.

Temporal or environmental hold-out. Fit on the cooler half of the environmental range and test on the warmer half. It directly measures whether the response extrapolates, which is the question a climate projection asks.

import numpy as np


def transferability_report(auc_random: float, auc_blocked: float,
                           mess_map: np.ndarray, suitable: np.ndarray) -> dict:
    """Summarise how much of the projected suitable area rests on extrapolation."""
    novel = mess_map < 0
    return {
        "auc_random": round(auc_random, 3),
        "auc_blocked": round(auc_blocked, 3),
        "optimism": round(auc_random - auc_blocked, 3),
        "suitable_area_px": int(np.nansum(suitable)),
        "suitable_on_novel_frac": round(
            float(np.nansum(suitable & novel) / max(np.nansum(suitable), 1)), 3),
        "median_mess_in_suitable": round(float(np.nanmedian(mess_map[suitable])), 1),
    }


rep = transferability_report(auc_random, auc_blocked, mess_map, suitable_mask)
assert rep["optimism"] < 0.15, f"random validation is badly optimistic: {rep}"
assert rep["suitable_on_novel_frac"] < 0.5, \
    f"over half the projected suitable area is extrapolation: {rep}"

The second assertion is the one that should stop a publication. If most of the future suitable area sits on environmental combinations that do not exist in the training data, the map is a statement about the model, not about the species.

Failure Modes & Gotchas

• Predictor bands in a different order between stacks. Silent, catastrophic, and caught in one line by a compatibility check.
• Different units or scaling. Temperature ×10 and precipitation in centimetres are both common in climate products and neither raises an error.
• Projecting with hinge features unclamped. Hinges extrapolate without bound, so a novel pixel can receive an arbitrarily high suitability.
• A single GCM presented as “the projection”. Inter-model spread is usually larger than the difference between emission scenarios at mid-century.
• Random cross-validation used to justify transfer. It measures interpolation among spatially autocorrelated points and routinely over-states transferability by 0.1 AUC or more.
• Reporting range shift without the extrapolation mask. A centroid that moved 300 km north may have moved entirely into pixels the model knows nothing about.

Performance & Scale Notes

Projection is cheap; the ensemble multiplies it. A continental stack at 1 km is roughly 10⁷ pixels per variable, so a 19-variable stack across nine GCMs and three periods is 27 projections over 190 million values. Two habits keep it tractable: project in windows rather than whole rasters, and compute MESS once per unique climate stack rather than once per model run, since it depends only on the predictors.

import numpy as np
import rasterio
from rasterio.windows import Window


def project_windowed(model, stack_paths: list[str], out_path: str,
                     train_min, train_max, block: int = 1024):
    """Project a model across a large stack window by window."""
    srcs = [rasterio.open(p) for p in stack_paths]
    try:
        profile = srcs[0].profile | {"count": 1, "dtype": "float32", "nodata": np.nan,
                                     "compress": "deflate", "predictor": 2, "tiled": True}
        with rasterio.open(out_path, "w", **profile) as dst:
            for j in range(0, srcs[0].height, block):
                for i in range(0, srcs[0].width, block):
                    win = Window(i, j, min(block, srcs[0].width - i),
                                 min(block, srcs[0].height - j))
                    arr = np.stack([s.read(1, window=win, masked=True).filled(np.nan)
                                    for s in srcs])
                    dst.write(project(model, arr, train_min, train_max), 1, window=win)
    finally:
        for s in srcs:
            s.close()
    return out_path

Frequently Asked Questions

Related

• MaxEnt Model Training & Tuning — produces the model this projects
• Projecting MaxEnt Models to Future Climate Scenarios in Python — the projection mechanics in detail
• Detecting Extrapolation with MESS Maps in Python — the novelty surface in detail
• Computing Range Shift and Centroid Metrics from SDM Projections — turning projections into reportable numbers
• Spatial Block Cross-Validation for SDM in Python — the validation that estimates transferability honestly

Up: Species Distribution Modeling with MaxEnt