Do range edges sit where climate changes abruptly?

Each step() call filters the table and writes the before/after count to cleaning_audit.json, which is where the funnel table above comes from. The issue-flag check is a substring match against the five flags listed in config.COORD_KILL_FLAGS.

def clean(df, audit):
    """Base cleaning, applied to every fork.  Records counts at each step."""
    def step(name, mask):
        before = len(step.df)
        step.df = step.df[mask].copy()
        audit[name] = {"removed": int(before - len(step.df)), "remaining": int(len(step.df))}
        print(f"  {name:28s} -removed {before-len(step.df):>9,}  -> {len(step.df):>10,}")
    step.df = df

    step("has_species", step.df["species"].notna() & (step.df["species"].str.len() > 0))
    step("has_coordinates",
         step.df["decimalLatitude"].notna() & step.df["decimalLongitude"].notna())

    # coordinate-invalidating issue flags only (see PREREGISTRATION.md 3)
    iss = step.df["issue"].fillna("")
    bad = np.zeros(len(step.df), dtype=bool)
    for f in C.COORD_KILL_FLAGS:
        bad |= iss.str.contains(f, regex=False).to_numpy()
    step("coord_issue_flags", ~bad)

    step("basis_of_record", ~step.df["basisOfRecord"].isin(C.DROP_BASIS_OF_RECORD))

    lat, lon = step.df["decimalLatitude"], step.df["decimalLongitude"]
    step("bbox_north_america",
         (lon >= C.BBOX["lon_min"]) & (lon <= C.BBOX["lon_max"]) &
         (lat >= C.BBOX["lat_min"]) & (lat <= C.BBOX["lat_max"]))

    # exact (0,0)-style and integer-degree "centroid" artefacts are covered by ZERO_COORDINATE and
    # by the bbox; we additionally drop records whose lat AND lon are both exact integers, a
    # classic sign of a country/state centroid rather than a locality.
    both_int = (lat == lat.round()) & (lon == lon.round())
    step("integer_degree_centroids", ~both_int)
    return step.df

Morphology is done with distance transforms: grow = cells within r of the range, shrink = cells more than r inside it. components() labels connected pieces (8-neighbour) and keeps the biggest until 95% of cells are covered. depth_of() pads the mask so a range touching the array edge isn't credited with infinite depth. In build_occ, note the order: land mask, then closing, then land mask again (closing can reach into the ocean), then fragments.

def dilate(m, r):
    if r <= 0: return m
    return ndi.distance_transform_edt(~m) <= r


def erode(m, r):
    if r <= 0: return m
    return ndi.distance_transform_edt(m) > r


def closing(m, r):
    if r <= 0: return m
    return erode(dilate(m, r), r)


def components(m, mode):
    if mode == "all" or not m.any():
        return m
    lab, n = ndi.label(m, structure=np.ones((3, 3)))
    if n <= 1:
        return m
    sizes = np.bincount(lab.ravel())[1:]
    if mode == "largest":
        return lab == (int(np.argmax(sizes)) + 1)
    order = np.argsort(sizes)[::-1]          # main95
    keep, cum, tot = [], 0, sizes.sum()
    for oi in order:
        keep.append(int(oi) + 1); cum += sizes[oi]
        if cum >= 0.95 * tot:
            break
    return np.isin(lab, keep)


def depth_of(occ):
    """Distance from each occupied cell to the nearest unoccupied cell.  Padded with False so a
    range abutting the array edge is not credited with infinite depth."""
    p = np.pad(occ, 1, constant_values=False)
    return ndi.distance_transform_edt(p)[1:-1, 1:-1]

def build_occ(cnt_sub, land_sub, thin, min_rec, comp, clo_cells, agg):
    """Occupied-cell mask for one fork combination, inside the species' padded bounding box."""
    if thin == "grid10km" and agg > 1:
        h, w = cnt_sub.shape
        ph, pw = -(-h // agg) * agg, -(-w // agg) * agg
        pad = np.zeros((ph, pw), np.int32); pad[:h, :w] = cnt_sub
        blk = pad.reshape(ph // agg, agg, pw // agg, agg).sum((1, 3))
        occ = np.repeat(np.repeat(blk >= min_rec, agg, 0), agg, 1)[:h, :w]
    else:
        occ = cnt_sub >= min_rec
    occ &= land_sub                      # climate is undefined off the land mask
    if not occ.any():
        return occ
    occ = closing(occ, clo_cells)
    occ &= land_sub
    occ = components(occ, comp)
    return occ

For each direction it averages three central differences with weights 1, 2, 1. A difference only counts if both of its cells are land, and the weights are renormalised over whatever counts. With a full window this is the textbook Horn (1981) operator.

def horn_gradient(z, mask, cell_km, min_valid=0.6):
    """Horn (1981) 3x3 weighted central difference, generalised to partial windows.

    dz/dx is a weighted mean of the three left-right differences (row weights 1,2,1); a row
    contributes only if BOTH of its flanking cells are land, and the weights are renormalised over
    the contributing rows.  With a complete window this is exactly the standard Horn operator.

    Returns gradient magnitude in <units of z> per km.
    """
    zf = np.where(mask, z, 0.0).astype(np.float64)
    v = mask.astype(np.float64)
    wts = [1.0, 2.0, 1.0]

    num = np.zeros_like(zf); den = np.zeros_like(zf)
    for r, wr in zip((-1, 0, 1), wts):
        zr = _shift(zf, r, 1, 0.0); zl = _shift(zf, r, -1, 0.0)
        vr = _shift(v, r, 1, 0.0);  vl = _shift(v, r, -1, 0.0)
        ok = vr * vl
        num += wr * ok * (zr - zl); den += wr * ok
    dzdx = np.divide(num, den * 2.0 * cell_km, out=np.full_like(num, np.nan), where=den > 0)

    num = np.zeros_like(zf); den = np.zeros_like(zf)
    for csh, wc in zip((-1, 0, 1), wts):
        zd = _shift(zf, 1, csh, 0.0); zu = _shift(zf, -1, csh, 0.0)
        vd = _shift(v, 1, csh, 0.0);  vu = _shift(v, -1, csh, 0.0)
        ok = vd * vu
        num += wc * ok * (zd - zu); den += wc * ok
    dzdy = np.divide(num, den * 2.0 * cell_km, out=np.full_like(num, np.nan), where=den > 0)

    vf = correlate(v, np.ones((3, 3)), mode="constant", cval=0.0) / 9.0
    g = np.hypot(dzdx, dzdy)
    g[(vf < min_valid) | ~mask] = np.nan
    return g, vf

Rather than looping over ~900,000 windows, it builds the sums a plane fit needs (Σ1, Σx, Σxz, …) for every cell with one convolution each, then solves all the 3 × 3 normal equations in batches. The residual sum of squares is Σz² − β·r, so roughness never needs the individual residuals. Ocean cells get weight zero, which is how partial coastal windows work.

def planar_fit(z, mask, k, cell_km, min_valid=0.6, components=False):
    """Weighted least-squares fit of a PLANE  z ~ a + b*x + c*y  over each k x k window,
    using only the land cells in that window.

    Returns
      slope_km : sqrt(b^2+c^2) converted to <units of z> per km   -- the 'planar' gradient fork
      rms      : root-mean-square residual from the fitted plane, in units of z
                 -- 'how badly does a flat plane describe this neighbourhood', i.e. the
                    departure-from-linearity that the hypothesis is actually about
      vf       : valid fraction of the window

    Note: the residual RMS is invariant to any invertible linear rescaling of the (x,y)
    coordinates, so it is identical whether x,y are measured in cells or in km.  Only the SLOPE
    needs the metric conversion.  This is what lets the same routine serve the native-lat/lon fork.
    """
    c0 = k // 2
    q = np.arange(k) - c0
    KX = np.tile(q, (k, 1)).astype(np.float64)          # column offset (x)
    KY = KX.T.copy()                                    # row offset (y)
    ones = np.ones((k, k))

    v = mask.astype(np.float64)
    zv = np.where(mask, z, 0.0).astype(np.float64)

    def cor(a, K):
        return correlate(a, K, mode="constant", cval=0.0)

    S0 = cor(v, ones); Sx = cor(v, KX); Sy = cor(v, KY)
    Sxx = cor(v, KX * KX); Syy = cor(v, KY * KY); Sxy = cor(v, KX * KY)
    Sz = cor(zv, ones); Szx = cor(zv, KX); Szy = cor(zv, KY); Szz = cor(zv * zv, ones)

    vf = S0 / (k * k)
    good = (vf >= min_valid) & mask & (S0 >= 3)
    idx = np.flatnonzero(good.ravel())

    slope = np.full(z.shape, np.nan, dtype=np.float64)
    rms = np.full(z.shape, np.nan, dtype=np.float64)
    bx = np.full(z.shape, np.nan, dtype=np.float64)
    by = np.full(z.shape, np.nan, dtype=np.float64)

    flat = [a.ravel() for a in (S0, Sx, Sy, Sxx, Syy, Sxy, Sz, Szx, Szy, Szz)]
    CH = 400_000
    for s in range(0, idx.size, CH):
        ii = idx[s:s + CH]
        s0, sx, sy, sxx, syy, sxy, sz, szx, szy, szz = (a[ii] for a in flat)
        M = np.empty((ii.size, 3, 3))
        M[:, 0, 0] = s0;  M[:, 0, 1] = sx;  M[:, 0, 2] = sy
        M[:, 1, 0] = sx;  M[:, 1, 1] = sxx; M[:, 1, 2] = sxy
        M[:, 2, 0] = sy;  M[:, 2, 1] = sxy; M[:, 2, 2] = syy
        r = np.stack([sz, szx, szy], axis=1)
        det = np.linalg.det(M)
        okd = np.abs(det) > 1e-8 * (np.abs(s0) ** 3 + 1.0)
        beta = np.full((ii.size, 3), np.nan)
        if okd.any():
            # numpy>=2 treats a trailing (N,3) as a matrix, not a batch of vectors
            beta[okd] = np.linalg.solve(M[okd], r[okd][:, :, None])[:, :, 0]
        rss = szz - np.einsum("ij,ij->i", beta, r)
        slope.ravel()[ii] = np.hypot(beta[:, 1], beta[:, 2]) / cell_km
        rms.ravel()[ii] = np.sqrt(np.maximum(rss, 0.0) / s0)
        bx.ravel()[ii] = beta[:, 1]     # d z / d(column), units of z per cell
        by.ravel()[ii] = beta[:, 2]     # d z / d(row),    units of z per cell
    if components:
        return slope, rms, vf, bx, by
    return slope, rms, vf

def classify(score, q, tail):
    """Return boolean masks (margin, core) for the outer/inner q of a marginality score.

    Ties matter here and are not hidden.  `depth` is a Euclidean distance transform, so a large
    block of cells sits at exactly 1.0 (every cell touching the range edge).  Taking
    `score <= quantile(score, q)` therefore returns MORE than a fraction q of cells whenever the
    q-th quantile falls on a tie -- for a geometrically thin range it can be a third of them.
    We do not break the tie: a cell one step from the edge is not more marginal than its
    neighbour just because of its array index.  The realised fractions are recorded per species
    and reported, and the equal-n rank alternative is run as a fork.
    """
    lo, hi = np.quantile(score, q), np.quantile(score, 1 - q)
    if tail == "low":
        return score <= lo, score >= hi
    return score >= hi, score <= lo

Species with fewer than 20 valid cells in either class get no value. The residual measure is already on a log scale and can be negative, so it is compared by difference rather than by log ratio; B's first version used the ratio for it too, which silently dropped half the species, and B found and fixed this before reading results.

def effect(vm, vc, stat, already_log=False):
    """Per-species margin-vs-core effect.

    'median'/'mean'  -> log ratio, so that +0.10 means the margin value is ~10.5% higher and the
                        statistic is symmetric under swapping the two groups.
    'cliffs_delta'   -> rank-based, bounded [-1,1], no distributional assumption, and insensitive
                        to the heavy right tail that both gradient and roughness have.

    `already_log=True` is required for metrics that are THEMSELVES residuals on a log scale
    (normrough5_resid).  Such a metric is centred on zero and is negative for roughly half of all
    cells, so log(median_margin / median_core) is undefined whenever either median is negative.
    Silently dropping those species is not a missing-data problem, it is a selection effect: the
    survivors are exactly the species whose margins AND cores are both rougher than the landscape
    regression predicts.  For these metrics the effect is the DIFFERENCE of medians, which is
    already a log ratio because the underlying quantity is a log.
    """
    vm = vm[np.isfinite(vm)]; vc = vc[np.isfinite(vc)]
    if vm.size < 20 or vc.size < 20:
        return np.nan, np.nan, np.nan
    if stat == "cliffs_delta":
        return cliffs_delta(vm, vc), np.median(vm), np.median(vc)
    f = np.median if stat == "median" else np.mean
    a, b = float(f(vm)), float(f(vc))
    if not (np.isfinite(a) and np.isfinite(b)):
        return np.nan, a, b
    if already_log:
        return a - b, a, b
    if a <= 0 or b <= 0:
        return np.nan, a, b
    return float(np.log(a / b)), a, b

def across_species(lrs, n_boot=10000, seed=1):
    """Inference with the SPECIES as the unit of analysis.

    Grid cells inside one range are strongly spatially autocorrelated; pooling ~10^5 of them gives a
    p-value that measures sample size, not nature.  So each species contributes exactly one number
    and we ask whether that distribution is centred away from zero.
    """
    x = np.asarray([v for v in lrs if np.isfinite(v)], float)
    out = {"n_species": int(x.size)}
    if x.size < 3:
        return {**out, "median": np.nan, "mean": np.nan, "ci_lo": np.nan, "ci_hi": np.nan,
                "wilcoxon_p": np.nan, "frac_positive": np.nan, "sign_p": np.nan}
    rng = np.random.default_rng(seed)
    boot = np.median(x[rng.integers(0, x.size, (n_boot, x.size))], axis=1)
    try:
        wp = float(stats.wilcoxon(x, alternative="two-sided").pvalue)
    except Exception:
        wp = np.nan
    npos = int((x > 0).sum())
    return {**out,
            "median": float(np.median(x)), "mean": float(np.mean(x)),
            "ci_lo": float(np.percentile(boot, 2.5)), "ci_hi": float(np.percentile(boot, 97.5)),
            "wilcoxon_p": wp,
            "frac_positive": npos / x.size,
            "sign_p": float(stats.binomtest(npos, x.size, 0.5).pvalue)}

Depth depends only on the shape, not on where it sits, so margin and core are computed once per rotation (36 angles) and reused for every placement. That trick is what makes 999 draws × 100 species × 12 measures run in about 25 minutes.

def rotations(occ_sub, q, n_angles, rng):
    """Pre-rotate a species mask and pre-classify margin/core for each angle."""
    out = []
    for ang in np.linspace(0, 360, n_angles, endpoint=False):
        rot = occ_sub if ang == 0 else ndi.rotate(occ_sub, ang, order=0, reshape=True,
                                                  prefilter=False, mode="constant", cval=False)
        rot = rot.astype(bool)
        if rot.sum() < C.MIN_OCCUPIED_CELLS:
            continue
        rr, cc = np.nonzero(rot)
        dep = depth_of(rot)[rr, cc].astype(np.float64)
        mm, co = classify(dep, q, "low")
        out.append({"shape": rot.shape, "r": rr, "c": cc, "m": mm, "k": co})
    return out

Up to 60 random positions are tried until 90% of the shape lands on land; if none works the draw is recorded as missing (10 of ~100,000 draws failed). The p-value compares the observed across-species median with the 999 null medians.

# ---- N1 displacement --------------------------------------------------------------
r0, c0 = row.min(), col.min()
sub = np.zeros((row.max() - r0 + 1, col.max() - c0 + 1), bool)
sub[row - r0, col - c0] = True
rots = rotations(sub, q, N_ANGLES, rng)
drawn = 0
while drawn < C.N_NULL_DRAWS:
    R = rots[rng.integers(len(rots))]
    h, w = R["shape"]
    if h >= H or w >= W:
        placement_fail += 1; drawn += 1; continue
    ok = False
    for _ in range(60):
        dr = rng.integers(0, H - h); dc = rng.integers(0, W - w)
        rr = R["r"] + dr; cc = R["c"] + dc
        if land[rr, cc].mean() >= C.NULL_MIN_LAND_OVERLAP:
            ok = True; break
    if not ok:
        placement_fail += 1
        for mk in M: n1[mk][drawn, si] = np.nan
        drawn += 1; continue
    for mk, arr in M.items():
        vals = arr[rr, cc].astype(np.float64)
        lr, _, _ = effect(vals[R["m"]], vals[R["k"]], "median",
                          already_log=mk[1] in LOG_METRICS)
        n1[mk][drawn, si] = lr
    drawn += 1

Both loops list three measures; normrough5_resid, which B's own null script computed, is not among them.

items, nrows = [], []
for v in ["bio1", "bio12", "pc12"]:
    for m in ["grad_horn3", "rough5", "normrough5_ratio"]:

for v in ["bio1", "bio12", "pc12"]:
    for m in ["grad_horn3", "rough5", "normrough5_ratio"]:
        k = f"{v}|{m}"; nv = NM[k]; cb = CB[k]; ce = CE[k]

The ratio divides roughness by steepness × cell size. The residual fits log(roughness) = a + b·log(steepness) across all land cells and keeps what is left, so the adjustment uses the slope the data actually have (0.75 for precipitation) rather than assuming 1.

def measure(self, var, metric):
    """Return the 2-D metric raster, deriving the normalised-roughness variants on demand."""
    if metric in ("grad_horn3", "grad_planar5", "rough5", "rough7", "rough9", "lap"):
        return self.layer(f"{var}__{metric}")
    g = self.layer(f"{var}__grad_horn3"); r = self.layer(f"{var}__rough5")
    ck = self.cell_km if np.isfinite(self.cell_km) else 1.0
    if metric == "normrough5_ratio":
        den = g * ck
        out = np.full(g.shape, np.nan, np.float32)
        ok = np.isfinite(g) & np.isfinite(r) & (den > 0)
        out[ok] = r[ok] / den[ok]
        return out
    # normrough5_resid: regress log(rough) on log(grad) over ALL land cells, keep residuals
    key = f"__resid__{var}"
    if key in self._cache:
        return self._cache[key]
    ok = np.isfinite(g) & np.isfinite(r) & (g > 0) & (r > 0)
    X = np.log(g[ok]); Y = np.log(r[ok])
    A = np.vstack([np.ones_like(X), X]).T
    coef, *_ = np.linalg.lstsq(A, Y, rcond=None)
    out = np.full(g.shape, np.nan, np.float32)
    out[ok] = Y - (coef[0] + coef[1] * X)
    self._cache[key] = out
    self._cache[key + "__coef"] = coef
    return out

Start from the primary, then add one entry per alternative value of each choice. The same loop in main() runs every specification against every variable and measure: 21 × 3 × 8 = 504 summaries.

def specs():
    """The primary specification plus every single-fork deviation from it."""
    base = dict(grid="ea", res=C.RES_PRIMARY, resamp="bilinear", occ=C.OCC_PRIMARY_KEY,
                cm="depth", q=C.D2_PRIMARY, stat="median", cls="quantile")
    out = [("PRIMARY", dict(base))]
    for r in C.RESOLUTIONS:
        if r != C.RES_PRIMARY: out.append((f"res={r}", {**base, "res": r}))
    out.append(("resamp=nearest", {**base, "resamp": "nearest"}))
    out.append(("grid=native_latlon", {**base, "grid": "nat"}))
    for v in C.occ_variants():
        if v != C.OCC_PRIMARY_KEY: out.append((f"occ:{v}", {**base, "occ": v}))
    for cm in ["alpha_depth", "cdist", "maha", "dens"]:
        out.append((f"marginality={cm}", {**base, "cm": cm}))
    for q in C.FORK_D2_QUANTILES:
        if q != C.D2_PRIMARY: out.append((f"q={q}", {**base, "q": q}))
    for st in ["mean", "cliffs_delta"]:
        out.append((f"stat={st}", {**base, "stat": st}))
    out.append(("classify=equal_n", {**base, "cls": "equal_n"}))
    return out

cd code
python3 01_clean_occurrences.py   # ~4 min; needs the 1 GB GBIF csv (its output is shipped)
python3 02_build_climate.py       # needs WorldClim 2.1 bioclim at 2.5m, 5m, 10m in data/climate/
python3 03_derivatives.py
python3 04_species_grids.py
python3 05_validate_ranges.py     # needs network once (taxonomic checklists)
python3 06_main_analysis.py       # the main result
python3 07_nulls.py --res 2.5m    # ~25 min
python3 08_verification.py
python3 09_figures.py && python3 10_report.py

# laptop profile: 12 species, 10', 99 null draws, ~20 s
export RANGELIMITS_PROFILE=laptop
for s in 02_build_climate 03_derivatives 04_species_grids 05_validate_ranges \
         06_main_analysis 07_nulls 08_verification 09_figures 10_report
do python3 ${s}.py || break; done

Everything on this page is read from B's output files; nothing in the analysis is recomputed. The script collects per-species effects, the fork summaries, B's saved null draws (taking the across-species median of each draw, as B's own p-values do), the verification results, a random sample of land cells for the scatter, and renders the range maps from B's saved masks using B's own classification rule. The code excerpts are pulled verbatim from B's scripts by function name at build time.

"""Pull arm B's own outputs into one JSON blob for the rewritten report.
Reads only files B's pipeline wrote; recomputes nothing except (a) medians of B's saved
null draws, (b) a random sample of land cells for the steepness-roughness scatter, and
(c) range maps rendered from B's saved species masks."""
import json, io, base64
import numpy as np, pandas as pd
from PIL import Image

B = "arm_B/data"
OUT = "B_data.json"
VARS = ["bio1", "bio12", "pc12"]
METRICS = ["grad_horn3", "rough5", "normrough5_ratio", "normrough5_resid"]
PRIMARY_OCC = "thin=cell|minrec=1|comp=main95|clokm=80"

out = {}

# ---- per-species effects, primary spec -------------------------------------------
per = pd.read_csv(f"{B}/results/per_species_LR.csv.gz")
pp = per[(per.spec == "PRIMARY") & per.metric.isin(METRICS)]
species = sorted(pp.species.unique())
out["species"] = species
out["lr"] = {f"{v}|{m}": [None if not np.isfinite(x) else round(float(x), 4)
                          for x in pp[(pp["var"] == v) & (pp.metric == m)]
                          .set_index("species").reindex(species).LR]
             for v in VARS for m in METRICS}

# ---- spec summaries (21 forks) ----------------------------------------------------
ss = pd.read_csv(f"{B}/results/spec_summary.csv")
ss = ss[ss.metric.isin(METRICS)]
out["specs"] = list(dict.fromkeys(ss.spec))
out["summ"] = {f"{r.spec}#{r['var']}|{r.metric}":
               [round(r["median"], 4), round(r.ci_lo, 4), round(r.ci_hi, 4),
                float(f"{r.wilcoxon_p:.4g}"), round(r.frac_positive, 3), round(r["mean"], 4)]
               for _, r in ss.iterrows()}

# ---- nulls: across-species median per draw, from B's saved draws -------------------
nj = json.load(open(f"{B}/results/nulls_2.5m.json"))
nd = np.load(f"{B}/results/null_draws_2.5m.npz")
out["nulls"] = {}
for v in VARS:
    for m in METRICS:
        k = f"{v}|{m}"
        d = np.nanmedian(nd[f"N1|{k}"], axis=1)
        d = d[np.isfinite(d)]
        s = nj["measures"][k]
        out["nulls"][k] = {"draws": [round(float(x), 4) for x in d],
                           "obs": round(s["observed_median_LR"], 4),
                           "centre": round(s["N1_null_median"], 4),
                           "lo": round(s["N1_null_p2.5"], 4), "hi": round(s["N1_null_p97.5"], 4),
                           "p": round(s["N1_p_onesided"], 3)}

# ---- verification: cluster bootstrap, coastal, effort ------------------------------
ver = json.load(open(f"{B}/results/verification.json"))
out["cluster"] = {k: [round(x, 4) for x in v["cluster_boot_ci"]]
                  for k, v in ver["A_cluster_bootstrap"].items()}
out["coastal"] = {k: round(v["median_LR_inland_only"], 4) for k, v in ver["E_coastal"].items()}
out["nonindep"] = ver["A_nonindependence"]

# ---- range validation + records per species ----------------------------------------
rv = json.load(open(f"{B}/results/range_validation.json"))["2.5m"]["per_species"]
aud = json.load(open(f"{B}/interim/cleaning_audit.json"))["per_species_after_cleaning"]
out["sp_meta"] = {s: {"records": aud.get(s),
                      "recall": rv.get(s, {}).get(PRIMARY_OCC, {}).get("recall"),
                      "jaccard": rv.get(s, {}).get(PRIMARY_OCC, {}).get("jaccard")}
                  for s in species}

# ---- steepness vs roughness scatter sample (bio12) --------------------------------
cl = np.load(f"{B}/interim/climate/ea_2.5m_bilinear.npz")
der = np.load(f"{B}/interim/derived/ea_2.5m_bilinear.npz")
land = cl["mask"]
g = der["bio12__grad_horn3"].astype(float); r = der["bio12__rough5"].astype(float)
ok = land & np.isfinite(g) & np.isfinite(r) & (g > 0) & (r > 0)
X, Y = np.log10(g[ok]), np.log10(r[ok])
b1, b0 = np.polyfit(X, Y, 1)
idx = np.random.default_rng(1).choice(X.size, 2500, replace=False)
out["scatter"] = {"x": [round(float(x), 3) for x in X[idx]], "y": [round(float(y), 3) for y in Y[idx]],
                  "b0": round(float(b0), 4), "b1": round(float(b1), 4), "n": int(X.size),
                  "r": round(float(np.corrcoef(X, Y)[0, 1]), 3)}

# ---- range maps: B's primary masks, margin/core by B's own classify rule ------------
sd = np.load(f"{B}/interim/species/species_2.5m.npz")
H, W = land.shape
PAL = [246, 247, 245,   # 0 ocean / outside (page ground)
       221, 226, 222,   # 1 land
       168, 184, 176,   # 2 occupied
       196, 86, 27,     # 3 margin
       47, 93, 138]     # 4 core
maps = {}
for s in species:
    p = f"{s}###{PRIMARY_OCC}###"
    row, col, dep = sd[p + "row"], sd[p + "col"], sd[p + "depth"].astype(float)
    lo, hi = np.quantile(dep, 0.10), np.quantile(dep, 0.90)
    img = np.where(land, 1, 0).astype(np.uint8)
    img[row, col] = 2
    img[row[dep >= hi], col[dep >= hi]] = 4
    img[row[dep <= lo], col[dep <= lo]] = 3
    pad = 40
    r0, r1 = max(0, row.min() - pad), min(H, row.max() + pad + 1)
    c0, c1 = max(0, col.min() - pad), min(W, col.max() + pad + 1)
    crop = img[r0:r1, c0:c1]
    f = int(np.ceil(max(crop.shape) / 560))          # pool large crops, class priority kept
    if f > 1:
        h2, w2 = -(-crop.shape[0] // f) * f, -(-crop.shape[1] // f) * f
        pc = np.zeros((h2, w2), np.uint8); pc[:crop.shape[0], :crop.shape[1]] = crop
        blk = pc.reshape(h2 // f, f, w2 // f, f)
        # priority: margin(3) > core(4) > occupied(2) > land(1) > ocean(0)
        prio = np.array([0, 1, 2, 4, 3])
        crop = np.argmax(np.stack([(blk == k).any((1, 3)) * prio[k] for k in range(5)]), 0)
        crop = crop.astype(np.uint8)
    im = Image.fromarray(crop, "P"); im.putpalette(PAL + [0] * (768 - len(PAL)))
    buf = io.BytesIO(); im.save(buf, "PNG", optimize=True)
    maps[s] = {"png": base64.b64encode(buf.getvalue()).decode(), "px_km": 5 * f,
               "cells": int(row.size), "margin": int((dep <= lo).sum()), "core": int((dep >= hi).sum())}
out["maps"] = maps

json.dump(out, open(OUT, "w"), separators=(",", ":"))
print("species", len(species), "bytes", len(json.dumps(out)),
      "maps KB", sum(len(m["png"]) for m in maps.values()) // 1024,
      "slope", b1, "r", out["scatter"]["r"])