Forestry GIS Python · Forestry · Ecology

Choosing a Binary Threshold for SDM Suitability Maps in Python

A continuous suitability surface answers “how suitable?”; every decision made from it answers “suitable or not?”. The threshold that converts one into the other is usually chosen in a line of code and rarely justified, yet it determines the mapped area, the reported range shift, and how much land a conservation designation covers. This guide makes that choice explicit, inside Model Validation & AUC Metrics, part of Species Distribution Modeling with MaxEnt.

When to use each rule

Rule Optimises Tends to give Use when
Maximum TSS Sensitivity + specificity − 1 A balanced, mid-range cut The general default for presence-background data
Maximum Kappa Agreement beyond chance Similar to TSS, less stable at low prevalence Comparability with older literature
Equal sensitivity and specificity Equal error rates Slightly lower than TSS The two error types cost the same
10th percentile training presence Omits the 10% most marginal presences A conservative, wider map Records include some misidentification or positional error
Minimum training presence Zero training omission A very wide map The species is genuinely everywhere it was recorded
Cost-weighted A stated cost ratio Whatever the ratio implies The consequences of the two errors differ, which is usual

Maximum TSS is the reasonable default. Reach for the 10th percentile rule when the occurrence data is imperfect enough that the most marginal records are probably errors, and for a cost-weighted rule whenever the two mistakes have genuinely different consequences — which, in conservation planning, they almost always do.

Where the common threshold rules fall on the same model Sensitivity falls and specificity rises as the suitability threshold increases from zero to one. Their crossing point marks the equal sensitivity-specificity threshold. A third curve, the true skill statistic, peaks slightly to the left of that crossing and marks the maximum TSS threshold. A fourth vertical line further left marks the tenth percentile training presence threshold, which is more permissive than either. 0 0.5 1.0 0 0.25 0.50 0.75 1.0 suitability threshold sensitivity specificity TSS 10th pct training presence max TSS equal sens/spec

Minimal reproducible example

import numpy as np


def threshold_candidates(y_true: np.ndarray, scores: np.ndarray,
                         presence_scores: np.ndarray | None = None,
                         n_steps: int = 500) -> dict:
    """Compute the standard threshold rules from evaluation data.

    y_true: 1 for presence, 0 for background; scores: model output for those points.
    presence_scores: suitability at the TRAINING presences, for percentile rules.
    """
    cuts = np.linspace(scores.min(), scores.max(), n_steps)
    pos, neg = y_true == 1, y_true == 0
    n_pos, n_neg = int(pos.sum()), int(neg.sum())
    assert n_pos > 0 and n_neg > 0, "need both classes to choose a threshold"

    sens = np.array([(scores[pos] >= c).sum() / n_pos for c in cuts])
    spec = np.array([(scores[neg] < c).sum() / n_neg for c in cuts])
    tss = sens + spec - 1.0

    out = {
        "max_tss": float(cuts[int(np.argmax(tss))]),
        "max_tss_value": float(tss.max()),
        "equal_sens_spec": float(cuts[int(np.argmin(np.abs(sens - spec)))]),
        "sensitivity_at_max_tss": float(sens[int(np.argmax(tss))]),
        "specificity_at_max_tss": float(spec[int(np.argmax(tss))]),
    }
    if presence_scores is not None and presence_scores.size:
        out["min_training_presence"] = float(np.nanmin(presence_scores))
        out["pct10_training_presence"] = float(np.nanpercentile(presence_scores, 10))
    return out


def cost_weighted_threshold(y_true, scores, cost_fn: float = 5.0,
                            cost_fp: float = 1.0, n_steps: int = 500) -> dict:
    """Threshold minimising a stated cost of false negatives against false positives."""
    cuts = np.linspace(scores.min(), scores.max(), n_steps)
    pos, neg = y_true == 1, y_true == 0
    costs = np.array([cost_fn * (scores[pos] < c).sum()
                      + cost_fp * (scores[neg] >= c).sum() for c in cuts])
    i = int(np.argmin(costs))
    return {"threshold": float(cuts[i]), "cost": float(costs[i]),
            "cost_ratio_fn_fp": cost_fn / cost_fp}

The cost-weighted rule is the honest one for applied work, because it forces the cost ratio into the open. Declaring that missing an occupied stand is five times worse than flagging an unoccupied one is a defensible management statement; silently accepting maximum TSS asserts a ratio of one without saying so.

What the choice costs in hectares

The threshold does not shift the map slightly — it can change mapped area by a factor of two or more, because suitability distributions are dense in the mid-range where the candidate thresholds sit.

import numpy as np
import rasterio


def area_by_threshold(raster_path: str, thresholds: dict) -> list[dict]:
    """Mapped suitable area under each candidate threshold."""
    with rasterio.open(raster_path) as src:
        arr = src.read(1, masked=True).astype("float64").filled(np.nan)
        cell_ha = abs(src.transform.a) * abs(src.transform.e) / 10_000.0
    valid = np.isfinite(arr)
    total = float(valid.sum() * cell_ha)
    rows = []
    for name, t in sorted(thresholds.items(), key=lambda kv: kv[1]):
        ha = float((valid & (arr >= t)).sum() * cell_ha)
        rows.append({"rule": name, "threshold": round(t, 4),
                     "area_ha": round(ha, 1),
                     "pct_of_domain": round(100.0 * ha / max(total, 1), 1)})
    return rows
Mapped area under four threshold rules, same model Four horizontal bars showing suitable area for one suitability raster under different threshold rules. The minimum training presence rule gives the largest area at about five hundred and twenty thousand hectares, the tenth percentile rule about four hundred thousand, maximum TSS about two hundred and forty thousand, and equal sensitivity and specificity about one hundred and sixty thousand. The span between the extremes is more than three times. min training presence 10th percentile max TSS equal sens/spec 520 000 ha 400 000 ha 240 000 ha 160 000 ha mapped suitable area One model, one raster, four defensible rules — and a 3.3× span in the number that gets quoted.

Parameter reference

Parameter Type Default Notes
evaluation data — held-out, spatially blocked Thresholds fitted on training data are optimistic
n_steps int 500 Resolution of the threshold search; more is cheap
cost_fn / cost_fp float stated ratio The declaration that makes a cost-weighted rule defensible
percentile float 10 5–20; higher is more conservative about marginal records
presence_scores array training presences Percentile rules need the training set, not the held-out one
reporting — rule + value + area All three, always

Expected output and verification

Compute every candidate, report the one you chose and why, and publish the area under each so a reader can substitute their own preference.

def threshold_report(y_true, scores, presence_scores, raster_path,
                     cost_fn: float = 5.0) -> dict:
    """Every candidate threshold with its mapped area, plus the chosen one."""
    cands = threshold_candidates(y_true, scores, presence_scores)
    cw = cost_weighted_threshold(y_true, scores, cost_fn=cost_fn)
    thresholds = {k: v for k, v in cands.items()
                  if k in ("max_tss", "equal_sens_spec",
                           "pct10_training_presence", "min_training_presence")}
    thresholds["cost_weighted"] = cw["threshold"]
    areas = area_by_threshold(raster_path, thresholds)
    span = max(a["area_ha"] for a in areas) / max(min(a["area_ha"] for a in areas), 1)
    return {"candidates": cands, "cost_weighted": cw,
            "areas": areas, "area_span_ratio": round(span, 2)}


rep = threshold_report(y_test, scores_test, presence_scores, "suitability.tif")
assert rep["candidates"]["max_tss_value"] > 0.4, \
    f"TSS too low for any threshold to be meaningful: {rep['candidates']}"
print(f"mapped area varies {rep['area_span_ratio']}x across candidate rules")

The TSS assertion is the gate that matters. Below about 0.4 the model barely separates presences from background, and no threshold produces a useful binary map — the honest output in that case is the continuous surface with a note that it does not support a binary interpretation.

The standing test is one the threshold machinery must fail:

import numpy as np

rng = np.random.default_rng(20260809)
# A model with no skill: scores drawn from the same distribution for both classes.
y = np.r_[np.ones(200), np.zeros(2000)].astype(int)
s = rng.beta(2, 5, size=y.size)

null = threshold_candidates(y, s)
assert null["max_tss_value"] < 0.15, \
    f"a skill-free model produced a usable TSS ({null['max_tss_value']:.2f})"

Feeding it scores that carry no information must produce a TSS near zero. An implementation that reports 0.4 on random data has a class-labelling or indexing bug, and every real threshold it produces is meaningless in the same way.

The two errors, their consequences, and the threshold they imply A two by two matrix of predicted against actual presence. The false negative cell notes that an occupied stand is not surveyed and may be logged, while the false positive cell notes that survey effort is spent on an unoccupied stand. Beneath, three cost ratios are listed with the threshold each implies: an equal ratio gives the maximum TSS threshold, a five to one ratio moves it lower, and a twenty to one ratio lower still. predicted suitable predicted unsuitable actually present actually absent correct surveyed and found false negative stand logged unsurveyed false positive survey effort wasted correct no effort spent cost ratio 1 : 1 cost ratio 5 : 1 cost ratio 20 : 1 threshold 0.42 — maximum TSS threshold 0.28 threshold 0.14 — nearly everything flagged Choosing maximum TSS silently asserts the first row.

Common pitfalls

  • Choosing the threshold on training data. Thresholds are as prone to optimism as any other fitted quantity; derive them from held-out, spatially blocked evaluation data.
  • Reporting area without the rule. An area figure with no stated threshold cannot be compared with anything, including a later version of itself.
  • Percentile rules from the wrong set. The 10th percentile rule is defined on training presence suitabilities, not on held-out ones.
  • Maximum kappa at low prevalence. Kappa is unstable when presences are rare relative to background, and its optimum can jump between runs.
  • Thresholding then averaging an ensemble. Threshold each member and average the binary maps to get agreement; thresholding an averaged surface discards the disagreement.

Frequently Asked Questions

Should I threshold at all?

Often not. A continuous suitability surface carries strictly more information, and many downstream uses — ranking sites, weighting a connectivity model, prioritising survey effort — work directly on it. Threshold when the decision is genuinely binary, and keep the continuous layer alongside.

Which rule do reviewers expect?

Maximum TSS and the 10th percentile training presence rule are both common and both accepted. What is expected without exception is that you name the rule and report the threshold value. A paper reporting “suitable area” with neither is asking to be queried.

How do I threshold when the model is an ensemble?

Threshold each member with its own rule, then combine the binary maps into an agreement layer. The commonly used shortcut — threshold the ensemble median — produces a single map that conceals whether the members agreed, which is the main thing an ensemble was built to reveal.

Up: Model Validation & AUC Metrics · Species Distribution Modeling with MaxEnt