Forestry GIS Python · Forestry · Ecology

Projecting MaxEnt Models to Future Climate Scenarios in Python

The mechanics of a future projection are ten lines of code, and almost all the risk sits in the data plumbing around them: which scenario, which climate model, which period, and whether the future stack really contains the variables the model thinks it does. This guide covers that plumbing, inside SDM Projection and Transferability, part of Species Distribution Modeling with MaxEnt.

When to use each scenario structure

Design Runs Reports Use when
One GCM, one scenario 1 A single map Never for publication; fine for a pipeline smoke test
Multi-GCM, one scenario 5–10 Median plus spread The usual minimum for a management report
One GCM, multi-scenario 3–4 Sensitivity to emissions When the question is explicitly about policy pathways
Multi-GCM × multi-scenario × period 50+ Median, agreement, and a scenario comparison A full assessment; plan the storage first

Inter-model spread usually exceeds inter-scenario spread until around mid-century, which has a practical consequence: if you can only afford one axis of variation, vary the climate model rather than the emissions scenario for near-term work, and add scenarios for end-of-century questions.

Inter-model spread against inter-scenario difference over time Projected suitable area is plotted against time from 2030 to 2090 for two emission scenarios. Each scenario is drawn as a band whose width is the spread across climate models. Near 2030 the two bands overlap almost completely, so the choice of scenario matters far less than the choice of climate model. By 2090 the bands have separated and the scenario becomes the dominant source of difference. 2030 2050 2070 2090 period suitable area low emissions high emissions bands overlap: model choice dominates bands separate: scenario dominates Band width is the spread across climate models; band separation is the scenario effect.

Minimal reproducible example

import itertools
from pathlib import Path

import numpy as np
import rasterio

VARIABLES = ["bio1", "bio4", "bio12", "bio15", "bio5", "bio6"]


def stack_paths(root: str, gcm: str, scenario: str, period: str) -> list[str]:
    """Paths for one climate future, in the model's variable order."""
    return [f"{root}/{gcm}/{scenario}/{period}/{v}.tif" for v in VARIABLES]


def verify_order(paths: list[str]) -> None:
    """Fail loudly if a path does not carry the variable the model expects."""
    for expected, path in zip(VARIABLES, paths):
        assert Path(path).stem == expected, \
            f"variable order mismatch: expected {expected}, got {Path(path).stem}"
        assert Path(path).exists(), f"missing predictor: {path}"


def project_future(model, root: str, gcm: str, scenario: str, period: str,
                   train_min: np.ndarray, train_max: np.ndarray,
                   out_dir: str, clamp: bool = True) -> str:
    """Project a fitted model onto one GCM-scenario-period combination."""
    paths = stack_paths(root, gcm, scenario, period)
    verify_order(paths)

    srcs = [rasterio.open(p) for p in paths]
    try:
        ref = srcs[0]
        arr = np.stack([s.read(1, masked=True).filled(np.nan) for s in srcs])
        n_vars, ny, nx = arr.shape
        flat = arr.reshape(n_vars, -1).T
        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")

        profile = ref.profile | {"count": 1, "dtype": "float32", "nodata": np.nan,
                                 "compress": "deflate", "predictor": 2, "tiled": True}
        out_path = f"{out_dir}/{gcm}_{scenario}_{period}.tif"
        Path(out_dir).mkdir(parents=True, exist_ok=True)
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(out.reshape(ny, nx), 1)
            dst.update_tags(gcm=gcm, scenario=scenario, period=period,
                            clamped=str(clamp), variables=",".join(VARIABLES))
    finally:
        for s in srcs:
            s.close()
    return out_path

The update_tags call costs nothing and solves the problem that otherwise appears six months later: a directory of fifty projection rasters whose filenames encode a convention nobody wrote down. Embedding the GCM, scenario, period, clamping state and variable list in the file itself makes every output self-describing.

Running the full grid is then a loop:

def project_grid(model, root: str, gcms: list[str], scenarios: list[str],
                 periods: list[str], train_min, train_max, out_dir: str) -> list[str]:
    """Every combination, with failures recorded rather than aborting the run."""
    written, failed = [], []
    for gcm, scenario, period in itertools.product(gcms, scenarios, periods):
        try:
            written.append(project_future(model, root, gcm, scenario, period,
                                          train_min, train_max, out_dir))
        except (AssertionError, OSError) as exc:
            failed.append({"gcm": gcm, "scenario": scenario,
                           "period": period, "error": str(exc)})
    if failed:
        print(f"{len(failed)} combination(s) failed:")
        for f in failed:
            print(f"  {f['gcm']}/{f['scenario']}/{f['period']}: {f['error']}")
    return written

Collecting failures rather than raising on the first one matters at this scale. A grid of fifty projections that aborts on combination three has wasted an hour and told you about one missing file; the version above tells you about all of them at once.

Parameter reference

Parameter Type Typical Notes
VARIABLES list 4–10 bioclim Must match training order exactly; fewer, uncorrelated variables transfer better
clamp bool True Produce both; the difference is a diagnostic
train_min / train_max array from training stack Per-variable, in the training data’s units
GCM count int 5–10 Below three there is no meaningful spread
scenario str 2–4 Emissions pathways; label with the standard identifiers
period str 2–4 20-year windows are conventional
output dtype str float32 float64 doubles storage for no gain

Expected output and verification

Each projection needs a sanity check, and the ensemble needs a consistency check. Both are cheap.

import numpy as np
import rasterio


def check_projection(path: str, present_path: str) -> dict:
    """Compare a future projection against the present-day one."""
    with rasterio.open(path) as f, rasterio.open(present_path) as p:
        fut = f.read(1, masked=True)
        pres = p.read(1, masked=True)
        tags = f.tags()
    valid = (~fut.mask) & (~pres.mask)
    return {
        "gcm": tags.get("gcm"),
        "scenario": tags.get("scenario"),
        "period": tags.get("period"),
        "mean_future": round(float(fut[valid].mean()), 4),
        "mean_present": round(float(pres[valid].mean()), 4),
        "corr_with_present": round(
            float(np.corrcoef(fut[valid].ravel(), pres[valid].ravel())[0, 1]), 3),
        "range": (round(float(fut[valid].min()), 4), round(float(fut[valid].max()), 4)),
    }


rep = check_projection("out/gcm1_ssp370_2070.tif", "out/present.tif")
assert 0.0 <= rep["range"][0] and rep["range"][1] <= 1.0, \
    f"suitability outside 0-1: {rep}"
assert rep["corr_with_present"] > 0.3, \
    f"future barely resembles the present projection — check band order: {rep}"

The correlation assertion is the band-order detector. A future stack whose variables are correctly ordered produces a map that is recognisably a shifted version of the present one; a shuffled stack produces something uncorrelated with it, which is otherwise very hard to notice on a map you have never seen before.

For the ensemble, the useful output is agreement rather than a single median:

import glob

import numpy as np
import rasterio


def ensemble_agreement(pattern: str, threshold: float, out_path: str) -> dict:
    """Fraction of ensemble members calling each pixel suitable."""
    paths = sorted(glob.glob(pattern))
    assert len(paths) >= 3, "an ensemble needs at least three members"
    arrs, profile = [], None
    for p in paths:
        with rasterio.open(p) as src:
            arrs.append(src.read(1, masked=True).filled(np.nan))
            profile = profile or src.profile
    stack = np.stack(arrs)
    agree = np.nanmean(stack >= threshold, axis=0).astype("float32")
    median = np.nanmedian(stack, axis=0).astype("float32")

    profile.update(count=2, dtype="float32", nodata=np.nan)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(median, 1)
        dst.write(agree, 2)
        dst.set_band_description(1, "median_suitability")
        dst.set_band_description(2, "ensemble_agreement")
    return {"n_members": len(paths),
            "median_agreement": round(float(np.nanmedian(agree)), 3)}
Median suitability against ensemble agreement Two schematic map panels of the same region. The left panel shows median suitability with a large contiguous area marked suitable, giving an impression of confidence. The right panel shows ensemble agreement over the same region: a core area where most members agree, surrounded by a wide fringe where fewer than half the members call the pixel suitable. The fringe is annotated as the part of the projection that a median map alone conceals. Median suitability Ensemble agreement suitable one confident-looking polygon ≥ 7 of 9 < 5 of 9 core plus a wide low-agreement fringe Publishing only the left panel presents minority-supported area as settled. Organising fifty projection rasters so they stay readable A directory tree groups projection outputs first by climate model, then by scenario, then by period, with one GeoTIFF per combination. Beside one file the embedded metadata tags are listed: the climate model name, the scenario, the period, whether clamping was applied, and the ordered variable list. A note records that the tags make each file self-describing even if it is moved out of the tree. projections/ gcm_a/ ssp245/ 2050.tif 2070.tif ssp370/ 2050.tif gcm_b/ … embedded tags gcm = gcm_a scenario = ssp245 period = 2050 clamped = True variables = bio1,bio4,bio12,bio15,bio5,bio6 The tree is convenience; the tags are what survive a file being copied somewhere else.

Common pitfalls

  • Variable order assumed rather than verified. Directory listings sort bio1, bio10, bio11, bio12, bio2… lexically, which is not the order anyone intended.
  • Mixed climate baselines. A model trained on one climatology and projected onto anomalies from another carries the baseline difference as a spurious trend.
  • Suitability compared across differently clamped runs. Clamped and unclamped surfaces are not the same quantity; label and separate them.
  • Storing outputs with no metadata. Fifty rasters named by convention alone become unusable the moment the convention is forgotten.
  • Averaging binary maps instead of continuous ones. Thresholding first and then averaging throws away the information that makes an agreement layer meaningful.

Frequently Asked Questions

Should I threshold before or after the ensemble?

Compute the continuous median from continuous members, and compute agreement by thresholding each member separately and averaging the binary results. Thresholding the median instead loses the distinction between a pixel where every model is marginally above the cut and one where half are far above and half far below.

How do I handle a GCM that is missing one variable?

Drop that GCM from the ensemble rather than substituting from another. Mixing sources within a single stack introduces a physically inconsistent climate — one variable’s rainfall paired with another’s temperature — which is worse than a smaller ensemble.

Do I need to re-tune the model for each projection?

No. The model is fitted once on present-day data; projection applies it unchanged. Re-tuning per scenario would mean fitting to the projection, which is circular. What does change per projection is the extrapolation surface, and that should be recomputed for every distinct climate stack.

Up: SDM Projection & Transferability · Species Distribution Modeling with MaxEnt