Forestry GIS Python · Forestry · Ecology

Computing Range Shift and Centroid Metrics from SDM Projections

Two suitability rasters — present and future — contain a report’s worth of numbers, and the ones people quote are almost always the same handful: how far the range moved, in which direction, how much it grew or shrank, and how much of it stayed put. Computing them takes a page of numpy; computing them in a way that survives scrutiny takes a projection choice and an extrapolation mask. This guide does both, inside SDM Projection and Transferability, part of Species Distribution Modeling with MaxEnt.

When to use each metric

Metric Units Answers Caveat
Centroid shift km, bearing Where is the range moving? A single point summarising a possibly bimodal distribution
Suitable area change ha or % Is the range growing or shrinking? Entirely dependent on the binary threshold
Range overlap (Jaccard/Schoener) 0–1 How much is the same ground? Insensitive to where the disagreement is
Elevation shift m Is it moving upslope? Needs a DTM on the same grid
Latitudinal shift km Poleward movement Distinct from centroid shift when the range is not compact
Stable / gain / loss areas ha Which ground matters for management? The most operationally useful, and least often reported

Report the gain/loss/stable decomposition alongside any headline shift figure. A centroid that moved 240 km north is compatible with a range that shifted wholesale and with one that stayed put while gaining a distant northern patch, and only the decomposition distinguishes them.

Decomposing a projected range shift into stable, lost and gained area Two overlapping range outlines are drawn on a map frame with a north arrow. The present range is the lower outline and the future range the upper one. Their intersection is shaded as stable area, the part of the present range outside the future one is marked as lost, and the part of the future range outside the present one is marked as gained. Two dots mark the suitability-weighted centroid of each range, and an arrow between them shows the shift distance and bearing. N lost gained stable shift 240 km, bearing 026° present range future range suitability-weighted centroid The same centroid shift is consistent with a wholesale move and with a static range plus a distant gain.

The centroid, and the projection it depends on

A suitability-weighted centroid is the mean position weighted by suitability:

Computing it in a geographic CRS averages degrees, which weights a degree of longitude at 60° north the same as one at the equator, so the centroid drifts. Every metric on this page must be computed in an equal-area projection — the discipline set out in Choosing an Equal-Area CRS for Forest Area Statistics.

Minimal reproducible example

import numpy as np
import rasterio


def weighted_centroid(path: str, threshold: float | None = None,
                      mess_path: str | None = None) -> dict:
    """Suitability-weighted centroid in the raster's own (equal-area) CRS."""
    with rasterio.open(path) as src:
        assert not src.crs.is_geographic, "compute centroids in an equal-area CRS"
        s = src.read(1, masked=True).astype("float64").filled(np.nan)
        transform, crs = src.transform, src.crs

    if mess_path:                                  # drop extrapolated pixels
        with rasterio.open(mess_path) as m:
            s = np.where(m.read(1) >= 0, s, np.nan)
    if threshold is not None:
        s = np.where(s >= threshold, s, np.nan)

    rows, cols = np.indices(s.shape)
    ok = np.isfinite(s) & (s > 0)
    if not ok.any():
        return {"status": "no_suitable_area"}

    w = s[ok]
    xs, ys = rasterio.transform.xy(transform, rows[ok], cols[ok])
    xs, ys = np.asarray(xs), np.asarray(ys)
    return {
        "x": float((w * xs).sum() / w.sum()),
        "y": float((w * ys).sum() / w.sum()),
        "crs": str(crs),
        "n_px": int(ok.sum()),
        "total_weight": float(w.sum()),
    }


def shift_metrics(present: dict, future: dict) -> dict:
    """Distance and bearing between two centroids in a projected CRS."""
    dx = future["x"] - present["x"]
    dy = future["y"] - present["y"]
    dist_km = float(np.hypot(dx, dy) / 1000.0)
    bearing = float((np.degrees(np.arctan2(dx, dy)) + 360.0) % 360.0)
    return {
        "shift_km": round(dist_km, 1),
        "bearing_deg": round(bearing, 1),
        "northward_km": round(dy / 1000.0, 1),
        "eastward_km": round(dx / 1000.0, 1),
    }

Note the argument order in arctan2(dx, dy) — reversed from the usual mathematical convention, because a compass bearing is measured clockwise from north rather than counter-clockwise from east. Getting this wrong produces a bearing reflected about the diagonal, which for a northward shift reads as eastward and is easy to miss.

Area change and the decomposition

import numpy as np
import rasterio


def area_decomposition(present_path: str, future_path: str, threshold: float,
                       mess_path: str | None = None) -> dict:
    """Stable, gained and lost area between two binarised projections."""
    with rasterio.open(present_path) as p, rasterio.open(future_path) as f:
        assert p.transform == f.transform, "rasters must share a grid"
        pres = p.read(1, masked=True).filled(np.nan)
        fut = f.read(1, masked=True).filled(np.nan)
        cell_ha = abs(p.transform.a) * abs(p.transform.e) / 10_000.0

    valid = np.isfinite(pres) & np.isfinite(fut)
    if mess_path:
        with rasterio.open(mess_path) as m:
            valid &= (m.read(1) >= 0)

    a = valid & (pres >= threshold)
    b = valid & (fut >= threshold)
    stable, gained, lost = (a & b), (~a & b), (a & ~b)
    union = (a | b).sum()
    return {
        "present_ha": round(float(a.sum() * cell_ha), 1),
        "future_ha": round(float(b.sum() * cell_ha), 1),
        "stable_ha": round(float(stable.sum() * cell_ha), 1),
        "gained_ha": round(float(gained.sum() * cell_ha), 1),
        "lost_ha": round(float(lost.sum() * cell_ha), 1),
        "net_change_pct": round(float(100.0 * (b.sum() - a.sum()) / max(a.sum(), 1)), 1),
        "jaccard": round(float(stable.sum() / max(union, 1)), 3),
    }

Every one of those numbers depends on threshold, which is why it belongs in the output and why the threshold choice deserves its own treatment in Choosing a Binary Threshold for SDM Suitability Maps. Reporting “the range shrinks by 34%” without the threshold is reporting an opinion.

The elevation shift is the same weighted mean applied to a terrain raster, and in mountainous regions it is often the more ecologically meaningful number:

def elevation_shift(present_path: str, future_path: str, dtm_path: str,
                    threshold: float) -> dict:
    """Suitability-weighted mean elevation of the range, present and future."""
    with rasterio.open(dtm_path) as d:
        dtm = d.read(1, masked=True).filled(np.nan)

    def weighted_mean_elev(path):
        with rasterio.open(path) as src:
            s = src.read(1, masked=True).filled(np.nan)
        ok = np.isfinite(s) & np.isfinite(dtm) & (s >= threshold)
        if not ok.any():
            return float("nan")
        return float((s[ok] * dtm[ok]).sum() / s[ok].sum())

    e0, e1 = weighted_mean_elev(present_path), weighted_mean_elev(future_path)
    return {"present_m": round(e0, 1), "future_m": round(e1, 1),
            "upslope_shift_m": round(e1 - e0, 1)}

Parameter reference

Parameter Type Typical Notes
CRS — equal-area Mandatory; degrees give a drifting centroid and meaningless areas
threshold float from a stated rule Every area figure depends on it; report it
mess_path str supplied Excluding extrapolated pixels usually shrinks the shift
weighting — continuous suitability Binary weighting discards magnitude information
DTM grid — matched Elevation shift needs the terrain on the same grid
ensemble handling — per member Compute metrics per member, then report their spread

That last row is the one most often skipped. Computing metrics on the ensemble median raster gives one number with no uncertainty; computing them per member and reporting the median and range gives the same central figure plus an honest interval.

Expected output and verification

The metrics are simple arithmetic, so verification is about catching the projection and orientation errors that make them wrong in ways that still look plausible.

import numpy as np


def verify_shift(present, future, decomposition) -> dict:
    """Internal consistency checks across the reported metrics."""
    m = shift_metrics(present, future)
    checks = {
        "areas_sum": abs(
            decomposition["stable_ha"] + decomposition["lost_ha"]
            - decomposition["present_ha"]) < 1.0,
        "future_sum": abs(
            decomposition["stable_ha"] + decomposition["gained_ha"]
            - decomposition["future_ha"]) < 1.0,
        "bearing_in_range": 0.0 <= m["bearing_deg"] < 360.0,
        "components_match": abs(
            np.hypot(m["eastward_km"], m["northward_km"]) - m["shift_km"]) < 0.2,
    }
    return {**m, "checks": checks}


res = verify_shift(cent_present, cent_future, decomp)
assert all(res["checks"].values()), f"metric inconsistency: {res['checks']}"

Then feed it a case it must get right by construction — a synthetic range translated by a known distance:

import numpy as np

# A synthetic suitability field, and the same field shifted 100 km north.
ny, nx, cell = 400, 400, 1000.0            # 1 km cells
yy, xx = np.mgrid[0:ny, 0:nx]
base = np.exp(-(((xx - 200) ** 2 + (yy - 260) ** 2) / (2 * 40.0 ** 2)))
moved = np.exp(-(((xx - 200) ** 2 + (yy - 160) ** 2) / (2 * 40.0 ** 2)))

cy_base = (base * yy).sum() / base.sum()
cy_moved = (moved * yy).sum() / moved.sum()
shift_km = (cy_base - cy_moved) * cell / 1000.0
assert abs(shift_km - 100.0) < 1.0, f"known 100 km shift measured as {shift_km:.1f} km"

A test that recovers a shift you constructed is worth more than any amount of inspection, because it exercises the sign conventions and the cell-size arithmetic together — the two places where a range-shift figure goes wrong by a factor rather than by a little.

Per-member range shift against the value from a median raster A horizontal axis of centroid shift distance from zero to four hundred kilometres. Nine points, one per ensemble member, are spread between about ninety and three hundred and ten kilometres. Below them a single marker shows the value computed from the ensemble median raster, which lands near two hundred kilometres and is presented without any interval, concealing the spread visible above it. 0 100 200 300 400 centroid shift (km) per member 90–310 km across the ensemble median raster 203 km, reported without an interval The median raster gives one number; the members give the number and its uncertainty. When the centroid lands where the species cannot live A projected range splits into two disjoint blocks separated by a wide unsuitable valley. The suitability-weighted centroid, computed over both blocks, falls squarely in the middle of that valley where suitability is near zero. A note records that the centroid is still a valid summary of the distribution's mass but is not a place, and that the block decomposition must be reported alongside it. block 1 block 2 weighted centroid suitability ≈ 0.02 here The centroid summarises the mass correctly and describes no habitat at all. Report patch count and the gain, loss and stable decomposition beside every centroid figure.

Common pitfalls

  • Computing in a geographic CRS. Averaging degrees weights high latitudes wrongly and every area is meaningless.
  • arctan2 arguments in mathematical order. Bearings run clockwise from north; the reversed convention silently reflects the direction.
  • Metrics from the ensemble median raster only. It yields a central estimate with no spread, which is exactly the information a projection most needs to carry.
  • No extrapolation mask. A range that “expands” into novel climate is expanding into pixels the model knows nothing about, and the number is not evidence.
  • A centroid summarising a bimodal range. When suitability splits into two disjoint blocks, the centroid can sit in unsuitable ground between them. Report the decomposition too.

Frequently Asked Questions

Should the centroid be weighted by suitability or computed on the binary range?

Weighted, in general — it uses all the information and is far less sensitive to the threshold. The binary centroid has one advantage worth knowing: it is directly comparable with published figures, most of which use binary ranges. Compute both when comparing to literature.

How do I report a range that fragments rather than moves?

With the decomposition and with a fragmentation statistic — patch count and mean patch size from a connected-components labelling. A fragmenting range can show a small centroid shift and a modest area change while becoming ecologically far less viable, and neither headline number captures that.

What if the projection domain does not cover the whole future range?

Then every metric is truncated at the domain edge and the shift is under-estimated. Check whether suitable pixels touch the boundary; if they do, extend the domain and recompute rather than reporting a figure that the study extent has capped.

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