Forestry GIS Python · Forestry · Ecology

Detecting Extrapolation with MESS Maps in Python

A projected suitability map looks identical whether the model is interpolating within familiar conditions or extrapolating into climates it has never seen. The multivariate environmental similarity surface — MESS — is the layer that tells the two apart, pixel by pixel, and it costs one pass over the predictor stack. This guide computes it and uses it, inside SDM Projection and Transferability, part of Species Distribution Modeling with MaxEnt.

When to use each novelty measure

Measure Detects Cost Limitation
MESS Values outside any single variable’s training range One pass Blind to novel combinations of familiar values
MoD (most dissimilar variable) Which variable is responsible Free with MESS Descriptive only
Mahalanobis distance Novel combinations, assuming ellipsoidal training data Covariance inverse Assumes multivariate normality
ExDet Both univariate and combinational novelty Two passes Less widely implemented
Convex hull membership Exact novelty in low dimensions Exponential in dimension Impractical above ~6 variables

Compute MESS always — it is cheap and it is what reviewers expect. Add Mahalanobis or ExDet when the projection domain is a different region rather than a different time, because that is where novel combinations of individually familiar values are common: the same temperature and the same rainfall, in a pairing that does not occur in the training area.

What MESS detects and what it misses A scatter plot of two predictor variables. The training data forms a diagonal cloud, indicating the two variables are correlated. Dashed lines mark the univariate range of each variable. Point A sits inside the cloud and is correctly flagged as familiar. Point B lies beyond the upper range of the horizontal variable and is correctly flagged as novel. Point C lies inside both univariate ranges but far off the diagonal, in a combination that never occurs in training, and MESS reports it as familiar. mean annual temperature annual precipitation A — familiar, MESS > 0 B outside range, MESS < 0 C — novel combination MESS says familiar MESS is univariate per variable, so point C passes every individual range test it applies.

Minimal reproducible example

import numpy as np


def mess_and_mod(train: np.ndarray, proj: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """MESS surface and the index of the most dissimilar variable per pixel.

    train: (n_train, n_vars) values at the training background;
    proj:  (n_pixels, n_vars) values in the projection domain.
    """
    n_train, n_vars = train.shape
    sims = np.empty((proj.shape[0], n_vars), dtype="float64")

    for j in range(n_vars):
        ref = np.sort(train[:, j])
        lo, hi = float(ref[0]), float(ref[-1])
        rng = max(hi - lo, 1e-12)
        v = proj[:, j]
        # Percentile of v within the training distribution.
        f = np.searchsorted(ref, v, side="right") / n_train * 100.0
        sims[:, j] = np.where(
            f == 0.0, (v - lo) / rng * 100.0,
            np.where(f <= 50.0, 2.0 * f,
                     np.where(f < 100.0, 2.0 * (100.0 - f),
                              (hi - v) / rng * 100.0)))

    mod = np.argmin(sims, axis=1).astype("int16")
    return sims.min(axis=1).astype("float32"), mod

The four branches encode the whole definition. Below the training minimum the score goes negative in proportion to how far below; in the lower half it rises to 100 at the median; in the upper half it falls symmetrically; above the maximum it goes negative again. A value of 100 means the pixel sits exactly at the training median for that variable; zero means it sits exactly on the boundary; negative means outside.

The most-dissimilar-variable layer is free and is the more actionable of the two. A MESS map says “novel”; the MoD map says “because precipitation is outside anything the model has seen”, which is a sentence you can act on — by adding data, dropping the variable, or restricting the projection domain.

import rasterio


def write_mess(mess: np.ndarray, mod: np.ndarray, shape, profile,
               variables: list[str], out_path: str) -> str:
    """Write MESS and the most-dissimilar-variable index as a two-band raster."""
    profile = profile | {"count": 2, "dtype": "float32", "nodata": np.nan,
                         "compress": "deflate", "tiled": True}
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(mess.reshape(shape).astype("float32"), 1)
        dst.write(mod.reshape(shape).astype("float32"), 2)
        dst.set_band_description(1, "mess")
        dst.set_band_description(2, "most_dissimilar_variable_index")
        dst.update_tags(variables=",".join(variables))
    return out_path

Parameter reference

Parameter Type Typical Notes
train sample array 10,000 background points The reference set defines the envelope; use the same background the model trained on
MESS threshold float 0 Negative is outside the range; stricter masks use −10 or −25
variables list model’s predictors Must match the projection stack exactly, in order
chunk size int 10⁶ pixels The percentile search is memory-hungry on a full continental stack
training source — background, not presences Presence-only training envelopes are far too narrow

The last row matters and is often wrong. MESS should be computed against the background the model saw — the environmental envelope it was fitted across — not against the presence records. Using presences alone produces a much narrower envelope and flags most of the projection domain as novel.

Expected output and verification

A correct MESS surface has properties you can assert. Computing it against the training domain itself is the strongest single test.

import numpy as np


def verify_mess(train: np.ndarray) -> dict:
    """MESS computed on the training data itself must be non-negative everywhere."""
    mess_self, _ = mess_and_mod(train, train)

    # A deliberately out-of-range probe must come back negative.
    far = train.copy()[:200]
    far[:, 0] = train[:, 0].max() + 5.0 * train[:, 0].std()
    mess_far, mod_far = mess_and_mod(train, far)

    return {
        "self_min": round(float(mess_self.min()), 3),
        "self_median": round(float(np.median(mess_self)), 1),
        "probe_max": round(float(mess_far.max()), 3),
        "probe_mod_is_var0": bool((mod_far == 0).all()),
    }


v = verify_mess(train_background)
assert v["self_min"] >= -1e-6, f"MESS negative inside the training data: {v}"
assert v["probe_max"] < 0, f"an out-of-range probe was not flagged: {v}"
assert v["probe_mod_is_var0"], f"the wrong variable was blamed: {v}"

Those three assertions are the standing test this machinery needs. The second and third deliberately feed it input it must reject — a probe pushed five standard deviations beyond the training maximum of variable zero — and check both that it refuses and that it attributes the refusal to the right variable. An implementation that scores everything positive will pass the first assertion and fail the other two.

Once the surface exists, the reporting use is to quantify how much of the projected suitable area depends on extrapolation:

import numpy as np


def novelty_report(mess: np.ndarray, suitability: np.ndarray,
                   threshold: float, mod: np.ndarray,
                   variables: list[str]) -> dict:
    """How much of the suitable area is novel, and which variable drives it?"""
    ok = np.isfinite(mess) & np.isfinite(suitability)
    suitable = ok & (suitability >= threshold)
    novel = ok & (mess < 0)
    counts = np.bincount(mod[novel & suitable], minlength=len(variables))
    return {
        "suitable_px": int(suitable.sum()),
        "novel_px": int(novel.sum()),
        "suitable_and_novel_frac": round(
            float((suitable & novel).sum() / max(suitable.sum(), 1)), 3),
        "driver": variables[int(counts.argmax())] if counts.sum() else None,
        "driver_share": round(float(counts.max() / max(counts.sum(), 1)), 3)
        if counts.sum() else None,
    }
MESS distributions for a near-term and a far-future projection Two overlaid distributions of MESS values across the same projection domain. The near-term projection is concentrated above zero with only a small tail below it. The far-future projection is shifted substantially to the left, with a large portion of its mass below zero, marking pixels whose conditions fall outside the training range. A vertical line at zero separates familiar from novel. −100 −50 0 50 100 MESS value pixels novel familiar 2030s 2080s The leftward shift is the projection moving off the edge of what the model was fitted on. Which variable makes each region novel A schematic projection domain is divided into three regions by which predictor is most dissimilar from the training range. The southern lowland is driven by mean annual temperature exceeding the training maximum. A central band is driven by precipitation seasonality. A northern upland region is driven by diurnal temperature range. A note records that the attribution tells you what to do: extend the training extent, drop the variable, or accept the limit. diurnal range precipitation seasonality mean annual temperature what to do extend training extent north check the variable earns its place accept the limit and mask A novelty map says where; the most-dissimilar-variable map says what to change.

Common pitfalls

  • Computing MESS against presence records. The envelope is then the species’ realised niche rather than the model’s training range, and almost everything reads as novel.
  • Masking rather than flagging. Setting novel pixels to NoData hides that the model made a prediction there. Keep the suitability and publish MESS beside it.
  • Treating MESS as sufficient. It is univariate per variable and cannot see a novel combination of familiar values, which is the dominant failure mode for spatial transfer.
  • Recomputing MESS per model run. It depends only on the predictor stacks, not on the fitted model, so one computation serves an entire ensemble sharing a climate future.
  • Comparing MESS values across different variable sets. Adding a variable can only lower the minimum, so a six-variable and a twelve-variable MESS are not comparable.

Frequently Asked Questions

What MESS threshold should mark a pixel as unreliable?

Zero is the definitional boundary and the usual choice. Some workflows use −10 or −25 to allow marginal extrapolation, which is defensible for variables whose response is known to be monotone and indefensible for quadratics. Whatever you pick, publish the continuous surface too, so a reader can apply their own.

Does clamping remove the need for MESS?

No — it changes what happens at novel pixels but not the fact that they are novel. Clamping holds the response flat at the boundary, which is a reasonable default, but a flat response is still an assumption rather than an observation. MESS is what records where that assumption was invoked.

Can I use MESS to choose which variables to keep?

Indirectly, and usefully. If the most-dissimilar-variable layer shows that one predictor drives almost all the novelty, ask whether that variable is earning its place. Dropping a weakly contributing variable that is responsible for most of the extrapolation often improves transferability more than any amount of regularisation tuning.

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