worldgen: initial public history

32 files changed, 4849 insertions, 0 deletions

diff --git a/mapgen/__init__.py b/mapgen/__init__.py
new file mode 100644
index 0000000..63fd969
--- /dev/null
+++ b/mapgen/__init__.py
@@ -0,0 +1,7 @@
+"""World generator."""
+import os
+
+# OpenBLAS threads spin ~0.1 s after each call before sleeping; with solves on several threads (and other jobs on the
+# machine) that spinning steals the cores doing the work. A short timeout changes nothing in how work is split, so
+# results stay the same bits (the thread *count* would not: see pipeline/graph notes). Only before numpy loads.
+os.environ.setdefault("OPENBLAS_THREAD_TIMEOUT", "4")
diff --git a/mapgen/check.py b/mapgen/check.py
new file mode 100644
index 0000000..f3f4639
--- /dev/null
+++ b/mapgen/check.py
@@ -0,0 +1,272 @@
+"""`mapgen.py check`: physical sanity checks (fail) + environment summary (info)."""
+from __future__ import annotations
+
+import json
+from pathlib import Path
+
+import numpy as np
+
+from . import plateaus as PL, zones as ZN
+from .climate import DEFAULTS as CLIMATE_DEFAULTS
+from .config import params
+from .environment import GROUND_NAMES, LANDFORM_NAMES, REGIONS
+from .fields import gravity_mod
+from .graph import gradient
+from .ice import ICE_NAMES
+from .sphere import gc_dist_km, latlon_to_xyz
+
+
+def check_land_fraction(z, area, target, tol=0.02):
+    f = area[z > 0].sum() / area.sum()
+    return None if abs(f - target) <= tol else f"land fraction {f:.3f} not within {tol} of target {target}"
+
+
+def check_hypsometry(z, area):
+    edges = np.arange(-11000, 12500, 500)
+    h, _ = np.histogram(z, bins=edges, weights=area)
+    c = (edges[:-1] + edges[1:]) / 2
+    ocean = (c >= -7000) & (c <= -1000)
+    landb = (c >= -500) & (c <= 3000)
+    io, il = np.flatnonzero(ocean)[np.argmax(h[ocean])], np.flatnonzero(landb)[np.argmax(h[landb])]
+    valley = h[io:il + 1].min()
+    if h[io] > 1.5 * valley and h[il] > 1.5 * valley:
+        return None
+    return "hypsometry not bimodal (expected separate ocean-floor and continent peaks)"
+
+
+def check_max_elevation(z, gmod, cap=12000.0):
+    lim = cap * np.clip(1.0 / gmod, 1.0, 3.0)
+    bad = z > lim + 1.0
+    return None if not bad.any() else f"{int(bad.sum())} cells above {cap:.0f} m (outside low-g zones)"
+
+
+def check_drainage(recv, z, endorheic, terminal_ok=None):
+    roots = recv == np.arange(len(recv))
+    bad = roots & (z > 0) & ~endorheic
+    if bad.any():
+        return f"{int(bad.sum())} land cells are sinks outside endorheic basins"
+    if terminal_ok is not None:
+        dry = roots & endorheic & ~terminal_ok
+        if dry.any():
+            return f"{int(dry.sum())} endorheic terminals are neither lake nor salt flat"
+    return None
+
+
+def check_no_uphill(zf, recv, z, endorheic, lake=None):
+    ar = np.arange(len(recv))
+    sel = (z > 0) & ~endorheic & (recv != ar)
+    if lake is not None:
+        sel &= ~lake                                  # flat lake surfaces: routing direction is a convention
+    bad = sel & (zf[recv] >= zf)
+    return None if not bad.any() else f"{int(bad.sum())} cells route uphill on the filled surface"
+
+
+def _band_means(g, f, edges):
+    a = np.abs(g.lat)
+    return [np.sum(f[(a >= lo) & (a < hi)] * g.area_km2[(a >= lo) & (a < hi)]) /
+            max(g.area_km2[(a >= lo) & (a < hi)].sum(), 1e-9) for lo, hi in zip(edges, edges[1:])]
+
+
+def check_temperature(g, t_mean):
+    m = _band_means(g, t_mean, list(range(0, 91, 10)))
+    ok = all(b <= a + 0.5 for a, b in zip(m, m[1:])) and m[0] > m[-1] + 10
+    return None if ok else f"zonal-mean temperature does not fall poleward: {np.round(m, 1).tolist()}"
+
+
+def check_rain_bands(g, p_ann, hadley_edge):
+    """Dry belt ≈ Hadley edge −2..+8°, storm track ≈ edge +12..+25° (scaled from Earth's 30° edge)."""
+    h = hadley_edge
+    eq = _band_means(g, p_ann, [0, 10])[0]
+    sub = _band_means(g, p_ann, [h - 2, h + 8])[0]
+    mid = _band_means(g, p_ann, [h + 12, h + 25])[0]
+    if eq > sub and mid > sub:
+        return None
+    return (f"no subtropical dry belt: P(0-10)={eq:.0f}, P({h - 2:.0f}-{h + 8:.0f})={sub:.0f}, "
+            f"P({h + 12:.0f}-{h + 25:.0f})={mid:.0f} mm/yr")
+
+
+def check_rain_shadow(g, data):
+    z = np.maximum(np.asarray(data["elevation_eroded_m"], dtype=np.float64), 0.0) / 1000.0
+    grad = gradient(g, z)
+    land = z > 0
+    ww, lw = [], []
+    for s in ("jun", "dec"):
+        up = np.sum(np.asarray(data[f"wind_{s}"], dtype=np.float64) * grad, axis=1)
+        p = np.asarray(data[f"P_{s}"])
+        ww.append(p[land & (up > 0.02)])
+        lw.append(p[land & (up < -0.02)])
+    w, l = np.concatenate(ww), np.concatenate(lw)
+    if len(w) == 0 or len(l) == 0:
+        return None
+    return None if w.mean() > l.mean() else f"windward slopes ({w.mean():.0f}) not wetter than lee ({l.mean():.0f})"
+
+
+def _shares(values, area, names, mask):
+    tot = area[mask].sum()
+    out = []
+    for i, nm in enumerate(names):
+        s = area[mask & (values == i)].sum() / max(tot, 1e-9)
+        if s >= 0.005:
+            out.append(f"{nm} {100 * s:.1f}%")
+    return ", ".join(out)
+
+
+
+def check_plateaus(g, data, plateaus, seed) -> list:
+    """Spec §8: every plateau's top within top_m (± relief), hidden ones ≤ −500 m, island ones with land."""
+    if not plateaus:
+        return []
+    z = np.asarray(data["elevation_eroded_m"], dtype=np.float64)
+    ocean, ids = np.asarray(data["ocean"]), np.asarray(data["plateau_id"])
+    out = []
+    for k, p in enumerate(plateaus):
+        idx = np.flatnonzero(ids == k)
+        if len(idx) == 0:
+            out.append(f"plateau {p['name']}: no cells at this resolution")
+            continue
+        _, b = PL.semi_axes(p)
+        core = idx[(1.0 - PL.rho(g.xyz[idx], p, seed, g.radius_km)) * b > PL.MARGIN_KM]
+        lo, hi = p["top_m"]
+        if len(core):
+            med = float(np.median(-z[core]))
+            if not lo - PL.RELIEF_M <= med <= hi + PL.RELIEF_M:
+                out.append(f"plateau {p['name']}: median top depth {med:.0f} m outside {lo:.0f}–{hi:.0f} m")
+        if p.get("islands", False):
+            isl = [g.cell_index(c["lat"], c["lon"]) for c in PL.features(p, seed, g.radius_km)["cones"] if c["island"]]
+            if ocean[idx].all() and all(ocean[i] for i in isl):
+                out.append(f"plateau {p['name']}: no island breaks the surface")
+        elif z[idx].max() > -500.0:
+            out.append(f"plateau {p['name']}: a hidden plateau reaches {z[idx].max():.0f} m (limit −500 m)")
+    return out
+
+
+def check_zones(g, data, max_o2_points=0.23, max_g=0.014, per_km=30.0) -> list:
+    """Spec §5, the gradual rule: steepest change per day's walk of the O₂ fraction (points) and of gravity."""
+    out = []
+    for key, scale, lim, unit in (("o2_fraction", 100.0, max_o2_points, "O₂ points"), ("gravity_g", 1.0, max_g, "g")):
+        if key not in data:
+            continue
+        f = np.asarray(data[key], dtype=np.float64) * scale
+        steep = float(np.max(np.abs(f[g.dst] - f[g.src]) / g.edge_km)) * per_km
+        if steep > lim + 1e-9:
+            out.append(f"zones: {key} changes {steep:.3f} {unit} per {per_km:.0f} km (limit {lim})")
+    return out
+
+
+def edited_areas(xyz, tect, radius_km, plateau_km=300.0, patch_km=1500.0):
+    """Cells the revision changes by design: plateaus + 300 km, land patches + 1,500 km, low-gravity zones
+    (their relief)."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    near = lambda c, r: gc_dist_km(xyz, latlon_to_xyz(*c), radius_km) < r
+    m = np.zeros(len(xyz), bool)
+    for p in tect.get("plateau", []):
+        m |= near(p["center"], PL.semi_axes(p)[0] * (1 + PL.EDGE_WARP) / (1 - PL.EDGE_WARP) + plateau_km)
+    for p in tect.get("land_patch", []):
+        m |= near(p["center"], p["radius_km"] * 1.3 + patch_km)
+    for z in tect.get("zone", []):
+        if z["field"] == "gravity":
+            m |= near(z["center"], z["radius_km"] / (1 - ZN.WARP))
+    return m
+
+
+O2_FIELDS = ("m_o2_zones", "o2_fraction", "po2_bar")        # what the base's O₂ zones change (no relief)
+
+
+def o2_areas(xyz, tect, radius_km):
+    """Cells the base's O₂ zones reach: compare skips only the O₂ fields there."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    m = np.zeros(len(xyz), bool)
+    for z in tect.get("zone", []):
+        if z["field"] == "o2":
+            m |= gc_dist_km(xyz, latlon_to_xyz(*z["center"]), radius_km) < z["radius_km"] / (1 - ZN.WARP)
+    return m
+
+
+def compare(old: dict, new: dict, exclude, extra=None) -> list:
+    """Per field, over the cells outside `exclude` (and outside extra[field] for that field): share changed, largest
+    change, 99th percentile."""
+    keep0 = ~np.asarray(exclude, bool)
+    lines = [f"cells compared: {int(keep0.sum())} of {len(keep0)} (outside the edited areas)"]
+    for k in sorted(set(old) & set(new)):
+        a, b = np.asarray(old[k]), np.asarray(new[k])
+        if a.shape != b.shape or a.shape[:1] != keep0.shape or a.dtype.kind not in "biuf" or b.dtype.kind not in "biuf":
+            continue
+        keep = keep0 & ~np.asarray(extra[k], bool) if extra and k in extra else keep0
+        a, b = a.astype(np.float64)[keep], b.astype(np.float64)[keep]
+        with np.errstate(invalid="ignore"):
+            diff = np.where((a == b) | (np.isnan(a) & np.isnan(b)), 0.0, np.abs(b - a))
+        if diff.ndim > 1:
+            diff = diff.reshape(len(diff), -1).max(axis=1)
+        if not len(diff):
+            lines.append(f"{k}: no cells to compare")
+            continue
+        lines.append(f"{k}: {np.mean(diff > 0):.1%} changed, max {diff.max():.4g}, p99 {np.percentile(diff, 99):.4g}")
+    return lines
+
+
+def compare_report(root: Path, res: int, old_path: Path) -> int:
+    from . import config as C
+    _, tect = C.load(root)
+    out = root / "out" / f"r{res}"
+    with np.load(old_path) as z:
+        old = {k: z[k] for k in z.files}
+    with np.load(out / "cells.npz") as z:
+        new = {k: z[k] for k in z.files}
+    for name, arrays in ((old_path, old), (out / "cells.npz", new)):
+        if "g_ids" not in arrays:
+            print(f"error: {name} has no g_ids: not the cells.npz of a world build")
+            return 2
+    if not np.array_equal(old["g_ids"], new["g_ids"]):
+        print("error: the two worlds have different cells (another resolution?)")
+        return 2
+    radius = float(json.loads((out / "cells_meta.json").read_text())["radius_km"])
+    o2 = o2_areas(new["g_xyz"], tect, radius)
+    for line in compare(old, new, edited_areas(new["g_xyz"], tect, radius), {k: o2 for k in O2_FIELDS}):
+        print(line)
+    return 0
+
+def info_lines(ctx):
+    g, d = ctx.grid, ctx.data
+    z = np.asarray(d["elevation_eroded_m"])
+    land = ~np.asarray(d["ocean"]) if "ocean" in d else z > 0
+    a = g.area_km2
+    sk = np.asarray(d["sk_land"]) > 0.5
+    iou = (land & sk).sum() / max((land | sk).sum(), 1)
+    sites = PL.site_report(g, d, ctx.tect["plateau"]) if ctx.tect.get("plateau") and "plateau_id" in d else []
+    return [f"cells {g.n} (H3 res {ctx.res}); land {100 * a[land].sum() / a.sum():.1f}% "
+            f"({a[land].sum() / 1e6:.0f}M km², Earth: 29%, 149M km²)",
+            f"sketch land overlap (IoU): {iou:.2f}",
+            "regions: " + _shares(np.asarray(d["hold_region"]), a, REGIONS, land),
+            "landforms: " + _shares(np.asarray(d["landform"]), a, LANDFORM_NAMES, land),
+            "ground: " + _shares(np.asarray(d["ground"]), a, GROUND_NAMES, land),
+            "ice: " + _shares(np.asarray(d["ice"]), a, ICE_NAMES, np.ones(g.n, bool))] + [f"site: {s}" for s in sites]
+
+
+def run_checks(ctx):
+    g, d, cfg = ctx.grid, ctx.data, ctx.cfg
+    z = np.asarray(d["elevation_eroded_m"], dtype=np.float64)
+    if "ocean" in d:   # interior lows below sea level are land: sign-correct z for the land/sea checks
+        z = np.where(np.asarray(d["ocean"]), np.minimum(z, -0.001), np.maximum(z, 0.001))
+    endo = np.asarray(d["endorheic"])
+    results = [
+        check_land_fraction(z, g.area_km2, cfg["build"]["land_fraction"]),
+        check_hypsometry(z, g.area_km2),
+        check_max_elevation(z, gravity_mod(cfg, d["m_gravity_zones"])),
+        check_drainage(np.asarray(d["recv"]), z, endo, np.asarray(d["lake"]) | np.asarray(d["salt_flat"])),
+        check_no_uphill(np.asarray(d["z_filled_m"]), np.asarray(d["recv"]), z, endo, np.asarray(d["lake"])),
+        check_temperature(g, np.asarray(d["T_mean"])),
+        check_rain_bands(g, np.asarray(d["P_ann"]), params(cfg, "climate", CLIMATE_DEFAULTS)["hadley_edge_deg"]),
+        check_rain_shadow(g, d),
+    ]
+    results += check_plateaus(g, d, ctx.tect.get("plateau", []), ctx.seed) + check_zones(g, d)
+    return [r for r in results if r], info_lines(ctx)
+
+
+def report(ctx) -> int:
+    fails, infos = run_checks(ctx)
+    for line in infos:
+        print("info:", line)
+    for f in fails:
+        print("FAIL:", f)
+    print("check:", "OK" if not fails else f"{len(fails)} failure(s)")
+    return 1 if fails else 0
diff --git a/mapgen/climate.py b/mapgen/climate.py
new file mode 100644
index 0000000..491675c
--- /dev/null
+++ b/mapgen/climate.py
@@ -0,0 +1,249 @@
+"""Stage `climate`: seasonal insolation → temperature, 3-cell winds + monsoons, moisture transport → rain."""
+from __future__ import annotations
+
+import numpy as np
+from scipy import sparse
+
+from . import ocean as OC
+from .config import params
+from .graph import bicgstab_jacobi, distance_to, gradient, nearest_source, pmap, smooth_km
+from .grid import rowdot_at
+from .sphere import east_north, latlon_to_xyz, rotate_about, tangent_dir
+
+S0 = 1361.0
+SEA_FREEZE_C = -1.8          # the exported SST never goes below sea water's freezing point (ice-covered sea)
+SEASONS = ("jun", "dec", "eq")
+DEFAULTS = {
+    "t_a": -55.2, "t_b": 0.3206, "t_c": -2.974e-4, "land_seasonal": 0.45, "land_summer": 0.9, "continentality_summer_max": 1.5, "ocean_seasonal": 0.15,
+    "continentality_km": 1500.0, "continentality_max": 1.8, "lapse_c_per_km": 6.5,
+    "current_c": 4.0, "current_reach_km": 600.0, "current_leak_km": 150.0, "heat_transport_km": 300.0,
+    "hadley_edge_deg": 20.0, "ferrel_edge_deg": 55.0, "itcz_shift_deg": 8.0,
+    "trade_u": -6.0, "trade_v": 2.0, "westerly_u": 8.0, "westerly_v": 1.0, "polar_u": -4.0, "polar_v": 1.0,
+    "monsoon_k": 1.0, "monsoon_length_km": 1500.0, "monsoon_speed_scale": 6000.0,
+    "coriolis_min_deg": 20.0, "coriolis_span_deg": 50.0,
+    "base_rate": 0.25, "conv_rate": 2.0, "front_rate": 0.8, "front_lat_deg": 40.0, "front_width_deg": 10.0, "itcz_width_deg": 8.0, "oro_rate": 20.0, "subsidence": 0.2,
+    "min_rate": 0.05, "cc_per_c": 0.07, "recycle": 0.83, "eddy_k_m2s": 2.2e6,
+    "eddy_wind_ms": 8.0,
+    "global_mean_mm": 1000.0,
+    "lock": False, "lock_at": [0.0, 0.0], "lock_day_c": 120.0, "lock_night_c": -200.0, "lock_wind_ms": 10.0,
+    "lock_melt_c": 0.0, "lock_melt_width_c": 25.0,
+}
+
+
+def t_of_q(q, P):
+    """Radiative-equilibrium-like surface temperature (°C) from insolation (W/m²); quadratic fit to Earth-like
+    zonal means for a 20° tilt: equator 27, 45° ≈ 15, 60° ≈ 3, pole ≈ −22 (concave: damps polar-day summers)."""
+    return P["t_a"] + P["t_b"] * q + P["t_c"] * q * q
+
+
+def declinations(tilt):
+    return {"jun": tilt, "dec": -tilt, "eq": 0.0}
+
+
+def insolation(lat, decl):
+    """Daily-mean top-of-atmosphere insolation (W/m²)."""
+    phi = np.radians(np.clip(lat, -89.9999, 89.9999))
+    d = np.radians(decl)
+    h0 = np.arccos(np.clip(-np.tan(phi) * np.tan(d), -1.0, 1.0))
+    return S0 / np.pi * (h0 * np.sin(phi) * np.sin(d) + np.cos(phi) * np.cos(d) * np.sin(h0))
+
+
+def zonal_mean(g, f, bin_deg=2.0):
+    b = np.floor((g.lat + 90.0) / bin_deg).astype(np.int64)
+    s = np.bincount(b, weights=f * g.area_km2)
+    w = np.bincount(b, weights=g.area_km2)
+    return (s / np.maximum(w, 1e-12))[b]
+
+
+def current_anomaly(g, land, P):
+    """Subtropical gyres: cold water off west coasts, warm off east coasts; leaks onto coastal land."""
+    if not land.any():
+        return np.zeros(g.n)
+    d, src = nearest_source(g, np.flatnonzero(land))
+    e, _ = east_north(g.xyz)
+    east_comp = np.sum(tangent_dir(g.xyz, g.xyz[np.maximum(src, 0)]) * e, axis=1)
+    band = np.sin(np.radians((np.clip(np.abs(g.lat), 10.0, 50.0) - 10.0) * 4.5))
+    a = -P["current_c"] * np.sign(east_comp) * band * np.exp(-d / P["current_reach_km"])
+    a[land] = 0.0
+    return smooth_km(g, a, P["current_leak_km"])
+
+
+def temperatures(g, z, land, P, tilt, cur=None):
+    decl = declinations(tilt)
+    Q = {s: insolation(g.lat, d) for s, d in decl.items()}
+    q_ann = (Q["jun"] + Q["dec"] + 2 * Q["eq"]) / 4
+    t_ann = t_of_q(q_ann, P)
+    dist_ocean = distance_to(g, ~land) if (~land).any() else np.full(g.n, 1e4)
+    cont = np.clip(1.0 + dist_ocean / P["continentality_km"], 1.0, P["continentality_max"])
+    cur = current_anomaly(g, land, P) if cur is None else cur
+    lapse = -P["lapse_c_per_km"] * np.maximum(z, 0.0) / 1000.0
+    def season(s):
+        raw = t_of_q(Q[s], P)
+        land_resp = np.where(raw > t_ann, P["land_summer"] * np.minimum(cont, P["continentality_summer_max"]),
+                             P["land_seasonal"] * cont)       # land heats faster in summer (dry, low heat capacity)
+        resp = np.where(land, land_resp, P["ocean_seasonal"])
+        return smooth_km(g, t_ann + resp * (raw - t_ann) + cur, P["heat_transport_km"]) + lapse
+    return dict(zip(SEASONS, pmap(season, SEASONS))), dist_ocean
+
+
+def biotemperature(tmean, trange, n=12):
+    ph = np.linspace(0.0, 2 * np.pi, n, endpoint=False)
+    t = np.asarray(tmean)[:, None] + (np.asarray(trange)[:, None] / 2) * np.sin(ph)[None, :]
+    return np.clip(t, 0.0, 30.0).mean(axis=1)
+
+
+def band_winds(g, itcz_lat, P):
+    """3-cell surface winds (m/s) relative to the thermal equator."""
+    phi = g.lat - itcz_lat
+    a = np.abs(phi)
+    sgn = np.where(phi >= 0, 1.0, -1.0)
+    s1 = 0.5 * (1 + np.tanh((a - P["hadley_edge_deg"]) / 3.0))
+    s2 = 0.5 * (1 + np.tanh((a - P["ferrel_edge_deg"]) / 4.0))
+    u = (1 - s1) * P["trade_u"] + (s1 - s2) * P["westerly_u"] + s2 * P["polar_u"]
+    v = sgn * (-(1 - s1) * P["trade_v"] + (s1 - s2) * P["westerly_v"] - s2 * P["polar_v"])
+    e, n = east_north(g.xyz)
+    return u[:, None] * e + v[:, None] * n
+
+
+def monsoon_winds(g, T_s, land, P):
+    """Thermal lows over hot land / highs over cold land, flow deflected by Coriolis."""
+    anom = np.where(land, T_s - zonal_mean(g, T_s), 0.0)
+    press = -P["monsoon_k"] * smooth_km(g, anom, P["monsoon_length_km"])
+    flow = -gradient(g, press)
+    theta = np.radians(P["coriolis_min_deg"] + P["coriolis_span_deg"] * np.abs(np.sin(np.radians(g.lat))))
+    return rotate_about(g.xyz, flow, -np.sign(g.lat) * theta) * P["monsoon_speed_scale"]
+
+
+def _rate(per_1000km, spacing_km):
+    return 1.0 - np.exp(-np.maximum(per_1000km, 0.0) * spacing_km / 1000.0)
+
+
+def eddy_mixing(g, P):
+    """kappa · (nbr − I): the per-step eddy mixing with the neighbours — the same for every season (build once)."""
+    kappa = P["eddy_k_m2s"] / (P["eddy_wind_ms"] * g.spacing_km * 1000.0)
+    nbr = sparse.csr_matrix((1.0 / g.counts[g.src], (g.src, g.dst)), shape=(g.n, g.n))
+    return kappa * (nbr - sparse.identity(g.n, format="csr"))
+
+
+def precipitation(g, wind, T, z, land, itcz_lat, P, mixing=None):
+    """Steady-state moisture transport along the wind on the cell graph; returns rain (relative units).
+    mixing: eddy_mixing(g, P), when several seasons share it."""
+    t = g.edge_tangents
+    out = np.maximum(rowdot_at(wind, g.src, t), 0.0)
+    tot = np.bincount(g.src, weights=out, minlength=g.n)
+    frac = np.where(tot[g.src] > 0, out / np.maximum(tot[g.src], 1e-12), 0.0)
+    Tm = sparse.csr_matrix((frac, (g.dst, g.src)), shape=(g.n, g.n))
+    stay = (tot <= 0).astype(np.float64)
+    upslope = np.maximum(np.sum(wind * gradient(g, np.maximum(z, 0.0) / 1000.0), axis=1), 0.0)
+    conv = P["conv_rate"] * np.exp(-((g.lat - itcz_lat) / P["itcz_width_deg"]) ** 2) * np.clip((T - 10.0) / 20.0, 0, 1)
+    subs = P["subsidence"] * np.exp(-((np.abs(g.lat - itcz_lat) - P["hadley_edge_deg"]) / 6.0) ** 2)
+    front = P["front_rate"] * np.exp(-((np.abs(g.lat - itcz_lat) - P["front_lat_deg"]) / P["front_width_deg"]) ** 2)
+    per = np.maximum(P["base_rate"] + conv + front + P["oro_rate"] * upslope - subs, P["min_rate"])
+    r = _rate(per, g.spacing_km)
+    evap = np.where(land, 0.0, np.exp(P["cc_per_c"] * (np.clip(T, -2.0, 35.0) - 25.0)))
+    # per advection step (one cell, time h/U): rain out, move downwind, eddy-mix with neighbours
+    mix = eddy_mixing(g, P) if mixing is None else mixing
+    eye = sparse.identity(g.n, format="csr")
+    keep = sparse.diags(1.0 - r)
+    recyc = sparse.diags(np.where(land, P["recycle"] * r, 0.0))   # land evapotranspiration returns rain
+    system = (eye - (Tm @ keep + sparse.diags(stay) @ keep) - mix - recyc).tocsr()
+    W, info = bicgstab_jacobi(system, evap, evap / np.maximum(r, 1e-6), 1.0 / system.diagonal(), 1e-7, 5000)
+    if info != 0:
+        raise ValueError(f"precipitation: moisture solve did not converge (info={info})")
+    return r * np.maximum(W, 0.0)
+
+
+def run_locked(ctx, g, z, land, P) -> dict:
+    """One face always to the sun: temperature by sun angle, no seasons; surface wind from night to the sun point."""
+    sub = latlon_to_xyz(*P["lock_at"])
+    mu = np.maximum(g.xyz @ sub, 0.0)
+    t = P["lock_night_c"] + (P["lock_day_c"] - P["lock_night_c"]) * mu ** 0.25
+    lapse = -P["lapse_c_per_km"] * np.maximum(z, 0.0) / 1000.0
+    T = smooth_km(g, t, P["heat_transport_km"]) + lapse
+    dist_ocean = distance_to(g, ~land) if (~land).any() else np.full(g.n, 1e4)
+    wind = P["lock_wind_ms"] * tangent_dir(g.xyz, np.broadcast_to(sub, g.xyz.shape))
+    rain = np.exp(-((T - P["lock_melt_c"]) / P["lock_melt_width_c"]) ** 2)   # meltwater and frost in the twilight ring
+    k = P["global_mean_mm"] / max(np.sum(rain * g.area_km2) / g.area_km2.sum(), 1e-12)
+    out = {}
+    for s in SEASONS:
+        out[f"wind_{s}"] = wind.astype(np.float32)
+        out[f"P_{s}"] = rain * k
+        out[f"T_{s}"] = T
+    out["P_ann"] = rain * k
+    out["T_mean"] = T
+    out["T_range"] = np.zeros(g.n)
+    out["T_min"] = T
+    out["biotemp"] = biotemperature(T, out["T_range"])
+    out["PET"] = 58.93 * out["biotemp"]
+    out["dist_ocean_km"] = dist_ocean
+    return out
+
+
+def _still_ocean(g, t_mean):
+    """No circulation (ocean disabled, locked world): zero currents/upwelling/productivity, SST = T_mean."""
+    z = np.zeros(g.n, np.float32)
+    return {"current": np.zeros((g.n, 3), np.float32), "current_speed": z, "sst": np.asarray(t_mean, np.float32),
+            "upwelling": z.copy(), "productivity": z.copy()}
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "climate", DEFAULTS)
+    tilt = float(ctx.cfg["planet"]["tilt_deg"])
+    (z,) = ctx.need("elevation_eroded_m")
+    z = z.astype(np.float64)
+    water = ctx.data.get("open_water", ctx.data.get("ocean"))            # big inland basins are water to the air
+    land = ~np.asarray(water) if water is not None else z > 0
+    O = params(ctx.cfg, "ocean", OC.DEFAULTS)
+    if P["lock"]:
+        out = run_locked(ctx, g, z, land, P)
+        out.update(_still_ocean(g, out["T_mean"]))
+        return out
+    sea = ~land
+    coupled = bool(O["enabled"]) and bool(sea.any())
+    T, dist_ocean = temperatures(g, z, land, P, tilt, cur=np.zeros(g.n) if coupled else None)
+    decl = declinations(tilt)
+    def season_wind(s):
+        itcz = P["itcz_shift_deg"] * decl[s] / max(tilt, 1e-9)
+        wind = band_winds(g, itcz, P)
+        return wind + monsoon_winds(g, T[s], land, P) if s != "eq" else wind
+    winds = dict(zip(SEASONS, pmap(season_wind, SEASONS)))
+    if coupled:                     # one pass: winds from current-free temperatures, then currents carry heat
+        day = float(ctx.cfg["planet"]["day_hours"])
+        w_ann = (winds["jun"] + winds["dec"] + 2 * winds["eq"]) / 4
+        u = OC.currents(g, sea, w_ann, O, day)
+        T_eq = (T["jun"] + T["dec"] + 2 * T["eq"]) / 4
+        T_s = OC.sst(g, sea, u, T_eq, O)
+        cur = smooth_km(g, np.where(sea, T_s - T_eq, 0.0), P["current_leak_km"])
+        T, _ = temperatures(g, z, land, P, tilt, cur=cur)
+    out = {}
+    mixing = eddy_mixing(g, P)
+    rain = dict(zip(SEASONS, pmap(lambda s: precipitation(g, winds[s], T[s], z, land,
+                                                           P["itcz_shift_deg"] * decl[s] / max(tilt, 1e-9), P, mixing),
+                                  SEASONS)))
+    del mixing
+    for s in SEASONS:
+        out[f"wind_{s}"] = winds[s].astype(np.float32)
+    ann = (rain["jun"] + rain["dec"] + 2 * rain["eq"]) / 4
+    k = P["global_mean_mm"] / max(np.sum(ann * g.area_km2) / g.area_km2.sum(), 1e-12)
+    for s in SEASONS:
+        out[f"P_{s}"] = rain[s] * k
+        out[f"T_{s}"] = T[s]
+    out["P_ann"] = ann * k
+    out["T_mean"] = (T["jun"] + T["dec"] + 2 * T["eq"]) / 4
+    out["T_range"] = np.abs(T["jun"] - T["dec"])
+    out["T_min"] = np.minimum(T["jun"], T["dec"])
+    out["biotemp"] = biotemperature(out["T_mean"], out["T_range"])
+    out["PET"] = 58.93 * out["biotemp"]
+    out["dist_ocean_km"] = dist_ocean
+    if coupled:
+        sst_c = np.where(sea, np.maximum(T_s, SEA_FREEZE_C), out["T_mean"])
+        w_up = OC.upwelling(g, sea, w_ann, O, day)
+        out.update({"current": u.astype(np.float32),
+                    "current_speed": np.linalg.norm(u, axis=1).astype(np.float32),
+                    "sst": sst_c.astype(np.float32),
+                    "upwelling": w_up.astype(np.float32),
+                    "productivity": OC.productivity(g, sea, w_up, z, out["T_range"], sst_c, O).astype(np.float32)})
+    else:
+        out.update(_still_ocean(g, out["T_mean"]))
+    return out
diff --git a/mapgen/config.py b/mapgen/config.py
new file mode 100644
index 0000000..a52986a
--- /dev/null
+++ b/mapgen/config.py
@@ -0,0 +1,269 @@
+"""Load and validate map/config/*.toml."""
+from __future__ import annotations
+
+import tomllib
+from pathlib import Path
+
+
+class ConfigError(ValueError):
+    pass
+
+
+MASK_DEFAULTS = {"o2_range": 0.5, "gravity_range": 0.7}
+
+REQUIRED = {
+    "planet": {
+        "radius_km": (1000.0, 100000.0),
+        "gravity_g": (0.1, 5.0),
+        "day_hours": (1.0, 10000.0),
+        "year_days": (1.0, 100000.0),
+        "tilt_deg": (0.0, 90.0),
+        "sea_level_pressure_bar": (0.01, 100.0),
+        "scale_height_km": (1.0, 100.0),
+        "o2_fraction": (0.0, 1.0),
+    },
+    "build": {
+        "seed": (0, 2**31 - 1),
+        "res_dev": (0, 8),
+        "res_final": (0, 8),
+        "raster_width": (64, 32768),
+        "preview_width": (64, 8192),
+        "land_fraction": (0.01, 0.99),
+    },
+}
+
+
+def _check_ranges(cfg: dict, schema: dict, src: str) -> None:
+    for sec, keys in schema.items():
+        if sec not in cfg:
+            raise ConfigError(f"{src}: missing section [{sec}]")
+        for k, (lo, hi) in keys.items():
+            if k not in cfg[sec]:
+                raise ConfigError(f"{src}: missing [{sec}].{k}")
+            v = cfg[sec][k]
+            if isinstance(v, bool) or not isinstance(v, (int, float)):
+                raise ConfigError(f"{src}: [{sec}].{k} must be a number, got {v!r}")
+            if not lo <= v <= hi:
+                raise ConfigError(f"{src}: [{sec}].{k}={v} outside [{lo}, {hi}]")
+
+
+def _latlon(v, what: str, src: str = "tectonics.toml") -> None:
+    ok = isinstance(v, list) and len(v) == 2 and all(isinstance(x, (int, float)) for x in v)
+    if not ok or not (-90 <= v[0] <= 90 and -180 <= v[1] <= 360):
+        raise ConfigError(f"{src}: {what} must be [lat, lon], got {v!r}")
+
+
+def validate_tectonics(t: dict) -> None:
+    plates = t.get("plate", [])
+    if len(plates) < 2:
+        raise ConfigError("tectonics.toml: need at least 2 [[plate]] tables")
+    seen = set()
+    for p in plates:
+        pid = p.get("id", "?")
+        for k in ("id", "seed", "kind", "motion"):
+            if k not in p:
+                raise ConfigError(f"tectonics.toml: plate {pid} missing {k}")
+        if pid in seen:
+            raise ConfigError(f"tectonics.toml: duplicate plate id {pid}")
+        seen.add(pid)
+        if p["kind"] not in ("continental", "oceanic"):
+            raise ConfigError(f"tectonics.toml: plate {pid} kind must be continental|oceanic")
+        _latlon(p["seed"], f"plate {pid}.seed")
+        m = p["motion"]
+        if not (isinstance(m, list) and len(m) == 2 and 0 <= m[1] <= 20):
+            raise ConfigError(f"tectonics.toml: plate {pid}.motion must be [azimuth_deg, speed 0..20 cm/yr]")
+    for kind in ("lip", "volcano", "hotspot", "microcontinent"):
+        for x in t.get(kind, []):
+            for k in ("name", "center"):
+                if k not in x:
+                    raise ConfigError(f"tectonics.toml: [[{kind}]] missing {k}")
+            _latlon(x["center"], f"{kind} {x['name']}.center")
+            if kind in ("lip", "volcano", "microcontinent") and not x.get("radius_km", 0) > 0:
+                raise ConfigError(f"tectonics.toml: {kind} {x['name']} needs radius_km > 0")
+            if kind == "hotspot" and not x.get("length_km", 0) > 0:
+                raise ConfigError(f"tectonics.toml: hotspot {x['name']} needs length_km > 0")
+
+
+ERAS_FILE = "eras.toml"          # events and eras: not part of the base world's inputs key (pipeline.inputs_key)
+ZONE_FIELDS = {"gravity_g": (0.05, 5.0), "o2_fraction": (0.0, 1.0), "pressure_bar": (0.01, 10.0),
+               "fire_reactivity": (0.0, 5.0)}                  # render.CONTINUOUS holds every value in range
+PLATEAU_SPACING_KM = 2500.0
+RESERVED_ERA_KEYS = {"order", "default"}
+
+
+def _num(v, what: str, lo: float, hi: float, src: str = "tectonics.toml") -> None:
+    if isinstance(v, bool) or not isinstance(v, (int, float)) or not lo <= v <= hi:
+        raise ConfigError(f"{src}: {what} must be a number in [{lo}, {hi}], got {v!r}")
+
+
+def _pair(v, what: str, lo: float, hi: float, src: str = "tectonics.toml") -> None:
+    ok = isinstance(v, list) and len(v) == 2 and all(isinstance(x, (int, float)) and not isinstance(x, bool) for x in v)
+    if not ok or not lo <= v[0] <= v[1] <= hi:
+        raise ConfigError(f"{src}: {what} must be [low, high] within [{lo}, {hi}], got {v!r}")
+
+
+def _named(items: list, kind: str, src: str = "tectonics.toml") -> set:
+    seen = set()
+    for x in items:
+        if not isinstance(x.get("name"), str):
+            raise ConfigError(f"{src}: [[{kind}]] needs a name")
+        if x["name"] in seen:
+            raise ConfigError(f"{src}: duplicate {kind} name {x['name']!r}")
+        seen.add(x["name"])
+    return seen
+
+
+def _event(e: dict) -> None:
+    w, src, kind = f"event {e['name']}", ERAS_FILE, e.get("kind")
+    if kind == "disintegrate":
+        _latlon(e.get("center"), f"{w}.center", src)
+        _num(e.get("radius_km"), f"{w}.radius_km", 1.0, 20000.0, src)
+        _num(e.get("depth_m"), f"{w}.depth_m", 1.0, 20000.0, src)
+    elif kind == "volcano":
+        _latlon(e.get("center"), f"{w}.center", src)
+        _num(e.get("radius_km"), f"{w}.radius_km", 1.0, 5000.0, src)
+        _num(e.get("peak_m"), f"{w}.peak_m", -11000.0, 12000.0, src)
+        if e.get("shape", "cone") not in ("cone", "shield", "caldera"):
+            raise ConfigError(f"{src}: {w}.shape must be cone|shield|caldera, got {e.get('shape')!r}")
+        if e.get("peak_mode", "above") not in ("above", "absolute"):
+            raise ConfigError(f"{src}: {w}.peak_mode must be above|absolute, got {e.get('peak_mode')!r}")
+    elif kind == "zone":
+        f = e.get("fields")
+        if not isinstance(f, dict) or not f:
+            raise ConfigError(f"{src}: {w}.fields must set at least one of {sorted(ZONE_FIELDS)}")
+        for k, v in f.items():
+            if k not in ZONE_FIELDS:
+                raise ConfigError(f"{src}: {w}.fields: unknown field {k!r} (known: {sorted(ZONE_FIELDS)})")
+            _num(v, f"{w}.fields.{k}", *ZONE_FIELDS[k], src)
+        shape = e.get("shape")
+        if shape == "circle":
+            _latlon(e.get("center"), f"{w}.center", src)
+            _num(e.get("radius_km"), f"{w}.radius_km", 1.0, 20000.0, src)
+        elif shape == "landmass":
+            _latlon(e.get("seed"), f"{w}.seed", src)
+            _pair(e.get("reach_km"), f"{w}.reach_km", 0.0, 20000.0, src)
+        else:
+            raise ConfigError(f"{src}: {w}.shape must be circle|landmass, got {shape!r}")
+        if "edge_km" in e:
+            _num(e["edge_km"], f"{w}.edge_km", 0.001, 5000.0, src)
+        elif e.get("profile") != "smooth":
+            raise ConfigError(f'{src}: {w} needs edge_km (a sharp edge) or profile = "smooth"')
+    else:
+        raise ConfigError(f"{src}: {w}.kind must be disintegrate|zone|volcano, got {kind!r}")
+
+
+def validate_revision(t: dict, radius_km: float) -> None:
+    """Plateaus, land patches, zones, events and eras."""
+    from .sphere import gc_dist_km, latlon_to_xyz
+    plateaus = t.get("plateau", [])
+    _named(plateaus, "plateau")
+    for p in plateaus:
+        w = f"plateau {p['name']}"
+        _latlon(p.get("center"), f"{w}.center")
+        _num(p.get("area_km2"), f"{w}.area_km2", 1.0e4, 5.0e7)
+        _num(p.get("elongation", 1.0), f"{w}.elongation", 1.0, 10.0)
+        _num(p.get("azimuth_deg", 0.0), f"{w}.azimuth_deg", -360.0, 360.0)
+        _pair(p.get("top_m"), f"{w}.top_m (metres below sea level, shallowest first)", 1.0, 11000.0)
+        if not isinstance(p.get("islands", False), bool):
+            raise ConfigError(f"tectonics.toml: {w}.islands must be true or false")
+        _num(p.get("vent", 1.0), f"{w}.vent", 0.0, 1.0)
+    for i, a in enumerate(plateaus):
+        for b in plateaus[i + 1:]:
+            d = float(gc_dist_km(latlon_to_xyz(*a["center"]), latlon_to_xyz(*b["center"]), radius_km))
+            if d < PLATEAU_SPACING_KM:
+                raise ConfigError(f"tectonics.toml: plateaus {a['name']} and {b['name']} are {d:.0f} km apart "
+                                  f"(need ≥ {PLATEAU_SPACING_KM:.0f} km)")
+    _named(t.get("land_patch", []), "land_patch")
+    for p in t.get("land_patch", []):
+        w = f"land_patch {p['name']}"
+        _latlon(p.get("center"), f"{w}.center")
+        _num(p.get("radius_km"), f"{w}.radius_km", 1.0, 20000.0)
+        _num(p.get("strength", 1.0), f"{w}.strength", 0.0, 2.0)
+        _num(p.get("edge_noise", 0.0), f"{w}.edge_noise", 0.0, 1.0)
+    _named(t.get("zone", []), "zone")
+    for z in t.get("zone", []):
+        w = f"zone {z['name']}"
+        if z.get("field") not in ("o2", "gravity"):
+            raise ConfigError(f"tectonics.toml: {w}.field must be o2|gravity, got {z.get('field')!r}")
+        _latlon(z.get("center"), f"{w}.center")
+        _num(z.get("radius_km"), f"{w}.radius_km", 1.0, 20000.0)
+        _num(z.get("v"), f"{w}.v", -1.0, 1.0)
+        if z.get("profile", "smooth") != "smooth":
+            raise ConfigError(f"tectonics.toml: {w}.profile must be smooth")
+    names = _named(t.get("event", []), "event", ERAS_FILE)
+    for e in t.get("event", []):
+        _event(e)
+    eras = t.get("eras")
+    if eras is None:
+        return
+    order = eras.get("order")
+    if not (isinstance(order, list) and order and all(isinstance(n, str) for n in order)
+            and len(set(order)) == len(order)) or RESERVED_ERA_KEYS & set(order):
+        raise ConfigError(f"{ERAS_FILE}: [eras].order must list distinct era names (not order/default), got {order!r}")
+    if eras.get("default") not in order:
+        raise ConfigError(f"{ERAS_FILE}: [eras].default {eras.get('default')!r} is not in order {order}")
+    extra = set(eras) - set(order) - RESERVED_ERA_KEYS
+    if extra:
+        raise ConfigError(f"{ERAS_FILE}: [eras] tables {sorted(extra)} are not in order {order}")
+    for n in order:
+        e = eras.get(n)
+        if not isinstance(e, dict) or not isinstance(e.get("label", n), str):
+            raise ConfigError(f"{ERAS_FILE}: missing [eras.{n}] (label, events)")
+        evs = e.get("events", [])
+        if not isinstance(evs, list) or not all(isinstance(x, str) for x in evs):
+            raise ConfigError(f"{ERAS_FILE}: [eras.{n}].events must be a list of event names")
+        if "years" in e:
+            _num(e["years"], f"[eras.{n}].years (since the era before it)", 0.0, 1.0e9, ERAS_FILE)
+        unknown = [x for x in evs if x not in names]
+        if unknown:
+            raise ConfigError(f"{ERAS_FILE}: era {n} names unknown events {unknown}")
+
+
+def era_events(tect: dict, name: str) -> list:
+    """The events of era `name`: those of every era before it in [eras].order, then its own (cumulative)."""
+    eras = tect.get("eras") or {}
+    order = eras.get("order", [])
+    if name not in order:
+        raise ConfigError(f"{ERAS_FILE}: unknown era {name!r} (eras: {order})")
+    by_name = {e["name"]: e for e in tect.get("event", [])}
+    out = []
+    for n in order[: order.index(name) + 1]:
+        out += [by_name[e] for e in eras.get(n, {}).get("events", [])]
+    return out
+
+
+def load(root: Path) -> tuple[dict, dict]:
+    wpath = root / "config" / "world.toml"
+    tpath = root / "config" / "tectonics.toml"
+    for p in (wpath, tpath):
+        if not p.exists():
+            raise ConfigError(f"missing {p}")
+    with open(wpath, "rb") as f:
+        cfg = tomllib.load(f)
+    with open(tpath, "rb") as f:
+        tect = tomllib.load(f)
+    for k in ("event", "eras"):
+        if k in tect:
+            raise ConfigError(f"tectonics.toml: [{k}] belongs in config/{ERAS_FILE} (editing an era must not rebuild "
+                              f"the base world)")
+    epath = root / "config" / ERAS_FILE
+    if epath.exists():
+        with open(epath, "rb") as f:
+            e = tomllib.load(f)
+        unknown = set(e) - {"event", "eras"}
+        if unknown:
+            raise ConfigError(f"{ERAS_FILE}: unknown tables {sorted(unknown)} (only [[event]] and [eras])")
+        tect = {**tect, **e}
+    _check_ranges(cfg, REQUIRED, "world.toml")
+    validate_tectonics(tect)
+    validate_revision(tect, float(cfg["planet"]["radius_km"]))
+    return cfg, tect
+
+
+def params(cfg: dict, section: str, defaults: dict) -> dict:
+    """Module defaults overridden by world.toml [section]; unknown keys are errors (typo guard)."""
+    over = cfg.get(section, {})
+    unknown = set(over) - set(defaults)
+    if unknown:
+        raise ConfigError(f"world.toml [{section}]: unknown keys {sorted(unknown)}")
+    return {**defaults, **over}
diff --git a/mapgen/crust.py b/mapgen/crust.py
new file mode 100644
index 0000000..345a1aa
--- /dev/null
+++ b/mapgen/crust.py
@@ -0,0 +1,97 @@
+"""Stage `crust`: continental vs oceanic crust, age classes, oceanic age."""
+from __future__ import annotations
+
+import numpy as np
+
+from .config import params
+from . import plateaus as PL
+from .graph import distance_to, smooth_km
+from .noise import fbm, name_seed, ridged
+from .plates import CONV, DIV
+from .sphere import azimuth_deg, gc_dist_km, latlon_to_xyz
+
+OCEANIC, CRATON, PRE_OROGEN, POST_OROGEN, RIFT, LIP, SCAR, VOLCANO = range(8)
+AGE_NAMES = ["oceanic", "craton (3g era)", "pre-Lightening orogen", "post-Lightening orogen",
+             "rift", "Lightening basalt province", "collapse scar", "overshoot volcano"]
+DEFAULTS = {"continental_threshold": 0.25, "shelf_km": 450.0, "edge_noise": 0.5, "edge_noise_freq": 4.0, "edge_noise_gain": 0.65, "orogen_width_km": 700.0, "rift_width_km": 200.0,
+            "pre_orogen_threshold": 0.72, "half_spreading_km_myr": 30.0, "max_ocean_age_myr": 200.0,
+            "scar_inner": 0.4, "scar_outer": 1.6, "scar_halfwidth_deg": 25.0}
+
+
+def center_dist(g, center):
+    return gc_dist_km(g.xyz, latlon_to_xyz(*center), g.radius_km)
+
+
+def lip_profile(g, lip, seed):
+    r = lip["radius_km"] * (1.0 + 0.25 * fbm(g.xyz, name_seed(seed, lip["name"]), 3, 6.0))
+    x = center_dist(g, lip["center"]) / r
+    return np.clip((1.0 - x) / 0.25, 0.0, 1.0)
+
+
+def _scar_azimuths(v, seed):
+    if "scar_azimuths" in v:
+        return list(v["scar_azimuths"])
+    rng = np.random.default_rng(name_seed(seed, v["name"]))
+    return list(rng.uniform(0, 360, size=2))
+
+
+def patch_field(g, patch: dict, seed: int):
+    """[[land_patch]]: land hint (strength) over a disc whose radius varies ± 25 % · edge_noise (fractal edge)."""
+    d = center_dist(g, patch["center"])
+    r = float(patch["radius_km"])
+    out = np.zeros(g.n)
+    near = d < 1.3 * r
+    if near.any():
+        w = np.clip(2.0 * fbm(g.xyz[near], name_seed(seed, patch["name"]), 6, 12.0), -1.0, 1.0)
+        out[near] = patch.get("strength", 1.0) * (d[near] < r * (1.0 + 0.25 * patch.get("edge_noise", 0.0) * w))
+    return out
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "crust", DEFAULTS)
+    land, land_hint, btype = ctx.need("sk_land", "m_land_hint", "bnd_type")
+    rough = (fbm(g.xyz, ctx.seed + 23, 7, P["edge_noise_freq"], gain=P["edge_noise_gain"])
+             if P["edge_noise"] > 0 else None)
+
+    def continental_of(raw):
+        field = smooth_km(g, raw, P["shelf_km"])
+        if rough is not None:     # fractal margins: bays, peninsulas, offshore continental fragments
+            f = np.clip(field, 0.0, 1.0)
+            near = 4.0 * f * (1.0 - f)                         # strongest at the margin; none deep inland / far offshore
+            field = field + P["edge_noise"] * near * rough / max(float(rough.std()), 1e-12)
+        else:
+            field = np.maximum(field, (raw > 0.5).astype(float))
+        cont = field > P["continental_threshold"]
+        for m in ctx.tect.get("microcontinent", []):
+            cont |= center_dist(g, m["center"]) < m["radius_km"]
+        return cont
+
+    raw = np.clip(land + 0.5 * land_hint, 0.0, 1.0)
+    continental_base = continental_of(raw)                     # without land patches: its sea level is the world's
+    patches = sum((patch_field(g, p, ctx.seed) for p in ctx.tect.get("land_patch", [])), np.zeros(g.n))
+    continental = continental_of(np.clip(raw + patches, 0.0, 1.0)) if patches.any() else continental_base.copy()
+    plateau_id = PL.cell_ids(g.xyz, ctx.tect.get("plateau", []), ctx.seed, g.radius_km)
+    continental |= plateau_id >= 0                             # sunken plateaus: continental crust that never rose
+    d_conv = distance_to(g, btype == CONV)
+    d_div = distance_to(g, btype == DIV)
+    age = np.full(g.n, CRATON, np.int8)
+    age[continental & (ridged(g.xyz, ctx.seed + 21, 4, 4.0) > P["pre_orogen_threshold"])] = PRE_OROGEN
+    age[continental & (d_div < P["rift_width_km"])] = RIFT
+    age[continental & (d_conv < P["orogen_width_km"])] = POST_OROGEN
+    age[~continental] = OCEANIC
+    for lip in ctx.tect.get("lip", []):
+        age[lip_profile(g, lip, ctx.seed) > 0] = LIP
+    for v in ctx.tect.get("volcano", []):
+        d = center_dist(g, v["center"])
+        r = v["radius_km"]
+        age[d < r] = VOLCANO
+        az = azimuth_deg(latlon_to_xyz(*v["center"]), g.xyz)
+        for a in _scar_azimuths(v, ctx.seed):
+            diff = np.abs((az - a + 180.0) % 360.0 - 180.0)
+            age[(diff < P["scar_halfwidth_deg"]) & (d > P["scar_inner"] * r) & (d < P["scar_outer"] * r)] = SCAR
+    oceanic_age = np.minimum(d_div / P["half_spreading_km_myr"], P["max_ocean_age_myr"])
+    ocean_age = np.where(continental & (plateau_id < 0), 0.0, oceanic_age)   # plateaus: the floor around them
+    return {"continental": continental, "age_class": age, "ocean_age_myr": ocean_age.astype(np.float32),
+            "d_conv_km": d_conv, "d_div_km": d_div, "continental_base": continental_base, "plateau_id": plateau_id,
+            "land_patch": patches.astype(np.float32)}
diff --git a/mapgen/elevation.py b/mapgen/elevation.py
new file mode 100644
index 0000000..47f0060
--- /dev/null
+++ b/mapgen/elevation.py
@@ -0,0 +1,197 @@
+"""Stage `elevation`: tectonic relief + Lightening provinces + hotspots + hints; sea level solved."""
+from __future__ import annotations
+
+import numpy as np
+
+from .config import params
+from .crust import PRE_OROGEN, SCAR, center_dist, lip_profile, name_seed
+from .fields import gravity_mod
+from .graph import OCEAN_MIN_KM2, distance_to, nearest_source, ocean_mask, smooth_km
+from .noise import fbm, ridged
+from . import plateaus as PL
+from .pipeline import StageError
+from .plates import edge_convergence
+from .sphere import east_north, gc_dist_km, great_circle_point, latlon_to_xyz
+
+OVER, SUB, COLLISION = 1, 2, 3
+DEFAULTS = {
+    "continental_base_m": 400.0, "continental_noise_m": 350.0,
+    "margin_km": 500.0, "shelf_m": -200.0, "slope_km": 250.0,
+    "coast_noise_m": 900.0, "coast_band_km": 600.0, "coast_noise_freq": 6.0,
+    "ridge_depth_m": 2500.0, "age_depth_coeff": 350.0, "abyss_m": 6500.0,
+    "rate_full_m_yr": 0.05, "min_rate_m_yr": 0.005,
+    "trench_depth_m": 4000.0, "trench_width_km": 70.0,
+    "arc_cont_m": 5500.0, "arc_ocean_m": 3500.0, "arc_offset_km": 200.0, "arc_width_km": 110.0,
+    "collision_peak_m": 9500.0, "collision_width_km": 260.0, "plateau_m": 4500.0, "plateau_km": 800.0,
+    "rift_depth_m": 1200.0, "rift_shoulder_m": 800.0,
+    "pre_orogen_m": 1200.0, "lip_m": 1500.0, "lip_step_m": 300.0,
+    "apron_m": 500.0, "scar_drop_m": 1500.0,
+    "hotspot_m": 5500.0, "hotspot_spacing_km": 150.0, "hotspot_radius_km": 70.0,
+    "hint_m": 2500.0, "hint_km": 80.0, "land_hint_m": 800.0, "spire_m": 2500.0, "detail_m": 250.0,
+    "max_land_m": 12000.0, "min_ocean_m": -11000.0,
+}
+
+
+def solve_sea_level(z, area, land_fraction):
+    order = np.argsort(-z)
+    cum = np.cumsum(area[order]) / area.sum()
+    k = min(int(np.searchsorted(cum, land_fraction)), len(z) - 1)
+    return z - z[order[k]]
+
+
+def solve_sea_level_connected(g, z, land_fraction, min_sea_km2=OCEAN_MIN_KM2, iters=30):
+    """Shift z so land = everything outside the connected ocean covers land_fraction (interior pits stay land)."""
+    area, tot = g.area_km2, g.area_km2.sum()
+    if abs(area[~ocean_mask(g, z, min_sea_km2)].sum() / tot - land_fraction) <= 0.5 * area.min() / tot:
+        return z                                   # already there: a flat sea floor would pull the bisection onto it
+    s0 = float(np.asarray(z)[np.argsort(-z)][min(int(np.searchsorted(np.cumsum(area[np.argsort(-z)]) / tot,
+                                                                           land_fraction)), len(z) - 1)])
+    lo, hi = s0 - 3000.0, s0 + 3000.0
+    for _ in range(iters):
+        mid = 0.5 * (lo + hi)
+        if area[~ocean_mask(g, z - mid, min_sea_km2)].sum() / tot > land_fraction:
+            lo = mid
+        else:
+            hi = mid
+    return z - 0.5 * (lo + hi)
+
+
+def _smoothstep(a, b, x):
+    t = np.clip((x - a) / (b - a), 0.0, 1.0)
+    return t * t * (3 - 2 * t)
+
+
+def roles(g, plate, vel, continental, ocean_age, plate_continental, min_rate):
+    conv, _ = edge_convergence(g, vel)
+    s, d = g.src, g.dst
+    e = (plate[s] != plate[d]) & (conv > min_rate)
+    cs, cd = continental[s[e]], continental[d[e]]
+    ks, kd = plate_continental[plate[s[e]]], plate_continental[plate[d[e]]]
+    over_s = np.where(cs & ~cd, True, np.where(~cs & cd, False,
+                      np.where(ks & ~kd, True, np.where(~ks & kd, False, ocean_age[s[e]] < ocean_age[d[e]]))))
+    role_e = np.where(cs & cd, COLLISION, np.where(over_s, OVER, SUB)).astype(np.int8)
+    order = np.argsort(conv[e])
+    role = np.zeros(g.n, np.int8)
+    role[s[e][order]] = role_e[order]
+    return role
+
+
+def _on_plate(g, plate, role_mask, other=None):
+    d, src = nearest_source(g, np.flatnonzero(role_mask))
+    ok = src >= 0
+    same = ok & (plate == plate[np.maximum(src, 0)])
+    if other is not None:
+        same |= ok & (plate == other[np.maximum(src, 0)])
+    return np.where(same, d, np.inf), np.maximum(src, 0)
+
+
+def hotspot_track(g, h: dict, vel, spacing_km: float):
+    """(unit vector, km from the active end) along a hotspot chain; the chain follows the plate's motion."""
+    c = latlon_to_xyz(*h["center"])
+    v = vel[g.cell_index(*h["center"])]
+    sp = np.linalg.norm(v)
+    t = v / sp if sp > 0 else east_north(c[None])[0][0]
+    return [(great_circle_point(c, t, k * spacing_km, g.radius_km), k * spacing_km)
+            for k in range(int(h["length_km"] // spacing_km) + 1)]
+
+
+def _relief(ctx, g, P, cont):
+    """Heights before the sea-level solve for the continental mask `cont` → (z, role, d_over, d_sub, d_coll)."""
+    R, seed = g.radius_km, ctx.seed
+    (plate, vel, pk, brate, bother, age, ocean_age, d_div, land_hint, sk_mtn, mtn_hint) = ctx.need(
+        "plate", "vel", "plate_continental", "bnd_rate", "bnd_other", "age_class", "ocean_age_myr", "d_div_km",
+        "m_land_hint", "sk_mountains", "m_mountain_hint")
+    f = lambda r: np.clip(r / P["rate_full_m_yr"], 0.3, 1.0)
+
+    # passive-margin profile: interior plateau → coastal ramp → shelf → slope → abyssal floor
+    d_in = distance_to(g, ~cont) if (~cont).any() else np.full(g.n, np.inf)
+    d_out = distance_to(g, cont) if cont.any() else np.full(g.n, np.inf)
+    ramp = _smoothstep(0.0, P["margin_km"], d_in)
+    z_cont = P["shelf_m"] + (P["continental_base_m"] - P["shelf_m"]) * ramp
+    z_cont = z_cont + P["continental_noise_m"] * fbm(g.xyz, seed + 31, 5, 3.0)
+    z_ocean = -np.minimum(P["ridge_depth_m"] + P["age_depth_coeff"] * np.sqrt(ocean_age), P["abyss_m"])
+    z_ocean = P["shelf_m"] + (z_ocean - P["shelf_m"]) * _smoothstep(0.0, P["slope_km"], d_out)
+    z = np.where(cont, z_cont, z_ocean)
+    # multi-scale coastal noise: sea level cuts it → bays, headlands, drowned valleys, offshore islands
+    d_edge = np.where(cont, d_in, d_out)
+    z += P["coast_noise_m"] * np.exp(-d_edge / P["coast_band_km"]) * fbm(g.xyz, seed + 91, 7, P["coast_noise_freq"])
+
+    role = roles(g, plate, vel, cont, ocean_age, pk, P["min_rate_m_yr"])
+    d_over, s_over = _on_plate(g, plate, role == OVER)
+    d_sub, s_sub = _on_plate(g, plate, role == SUB)
+    d_coll, s_coll = _on_plate(g, plate, role == COLLISION, other=bother)
+
+    z += -P["trench_depth_m"] * f(brate[s_sub]) * np.exp(-(d_sub / P["trench_width_km"]) ** 2)
+    arc_h = np.where(cont, P["arc_cont_m"], P["arc_ocean_m"])
+    z += arc_h * f(brate[s_over]) * np.exp(-((d_over - P["arc_offset_km"]) / P["arc_width_km"]) ** 2)
+    w = P["collision_width_km"]
+    z += P["collision_peak_m"] * f(brate[s_coll]) * np.exp(-(d_coll / w) ** 2)
+    plateau = _smoothstep(0.5 * w, 1.5 * w, d_coll) * (1 - _smoothstep(0.7 * P["plateau_km"], P["plateau_km"], d_coll))
+    z += np.where(cont, P["plateau_m"] * f(brate[s_coll]) * plateau, 0.0)
+    rift = -P["rift_depth_m"] * np.exp(-(d_div / 60.0) ** 2) + P["rift_shoulder_m"] * np.exp(-((d_div - 120.0) / 60.0) ** 2)
+    z += np.where(cont, rift, 0.0)
+
+    z += np.where(age == PRE_OROGEN, P["pre_orogen_m"] * ridged(g.xyz, seed + 21, 4, 4.0), 0.0)
+    for lip in ctx.tect.get("lip", []):
+        prof = lip_profile(g, lip, seed)
+        z += np.floor(P["lip_m"] * prof / P["lip_step_m"]) * P["lip_step_m"]
+    for v in ctx.tect.get("volcano", []):
+        d = center_dist(g, v["center"])
+        r, H = v["radius_km"], v.get("height_m", 7000.0)
+        z += H * np.clip(1 - d / r, 0, 1) ** 1.5 - 0.35 * H * np.exp(-(d / (0.12 * r)) ** 2)
+        z += P["apron_m"] * np.exp(-((d - r) / (0.3 * r)) ** 2) * (0.5 + 0.5 * fbm(g.xyz, name_seed(seed, v["name"]), 3, 20.0))
+    z -= np.where(age == SCAR, P["scar_drop_m"], 0.0)
+
+    for h in ctx.tect.get("hotspot", []):
+        peaks = np.zeros(g.n)
+        for p, s in hotspot_track(g, h, vel, P["hotspot_spacing_km"]):
+            dk = gc_dist_km(g.xyz, p, R)
+            hk = P["hotspot_m"] * (1 - s / h["length_km"])
+            peaks = np.maximum(peaks, hk * np.clip(1 - dk / P["hotspot_radius_km"], 0, 1) ** 1.2)
+        z += peaks
+
+    hint = smooth_km(g, np.clip(sk_mtn + mtn_hint, 0, 1), P["hint_km"])
+    z += np.where(cont, P["hint_m"] * hint, 0.0)
+    z += P["land_hint_m"] * land_hint
+    z += P["detail_m"] * fbm(g.xyz, seed + 41, 6, 12.0)
+    return z, role, d_over, d_sub, d_coll
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "elevation", DEFAULTS)
+    continental, sk_land, land_hint, m_grav, m_lock = ctx.need("continental", "sk_land", "m_land_hint",
+                                                               "m_gravity_zones", "m_lock")
+    plateau_id = np.asarray(ctx.data.get("plateau_id", np.full(g.n, -1)))
+    cont = np.asarray(continental) & (plateau_id < 0)          # plateaus: ocean until their surfaces are set (§3)
+    base = np.asarray(ctx.data.get("continental_base", continental)) & (plateau_id < 0)
+    z, role, d_over, d_sub, d_coll = _relief(ctx, g, P, cont)
+
+    target = ctx.cfg["build"]["land_fraction"]
+    cont_area = g.area_km2[base].sum() / g.area_km2.sum()
+    if target > cont_area + 0.005:
+        raise StageError(f"land_fraction {target} exceeds continental crust area {cont_area:.3f}: sea level would "
+                         f"lift ocean ridges into land; lower [build] land_fraction or raise [crust] shelf_km")
+    min_sea = ctx.cfg.get("erosion", {}).get("min_sea_km2", OCEAN_MIN_KM2)
+    sea_ref = np.zeros(g.n, bool)                              # sea in the world without land patches
+    if np.array_equal(base, cont):
+        z = solve_sea_level_connected(g, z, target, min_sea)
+    else:   # land patches (the eastern continent made whole) must not move every other coast: the sea level of the
+        z_raw = _relief(ctx, g, P, base)[0]                  # world without them, applied to the world with them
+        z_ref = solve_sea_level_connected(g, z_raw, target, min_sea)
+        z = z - float(np.mean(z_raw - z_ref))
+        sea_ref = ocean_mask(g, z_ref, min_sea)
+        patch = smooth_km(g, np.asarray(ctx.data.get("land_patch", np.zeros(g.n)), dtype=np.float64), P["hint_km"])
+        z += P["land_hint_m"] * patch                        # a patch is a land hint in height too (bare crust
+                                                             # sits below the world's sea level: land)
+    mult = np.clip(1.0 / gravity_mod(ctx.cfg, m_grav), 1.0, 3.0)
+    spires = P["spire_m"] * (mult - 1) * ridged(g.xyz, ctx.seed + 71, 4, 40.0)
+    z = np.where(z > 0, z * mult + spires, z)
+    target_land = (sk_land + 0.5 * land_hint) > 0.5
+    lock = m_lock > 0.5
+    z = np.where(lock & target_land, np.maximum(z, 50.0), np.where(lock & ~target_land, np.minimum(z, -50.0), z))
+    z = PL.apply(g, z, ctx.tect.get("plateau", []), plateau_id, ctx.seed)     # sunken plateaus (after the solve)
+    z = np.clip(z, P["min_ocean_m"], P["max_land_m"] * mult)
+    added = ~ocean_mask(g, z, min_sea) & (sea_ref | (plateau_id >= 0))       # land the sea-level solve never saw
+    return {"elevation_m": z.astype(np.float32), "role": role, "d_over_km": d_over, "d_sub_km": d_sub,
+            "d_coll_km": d_coll, "land_added": added}
diff --git a/mapgen/environment.py b/mapgen/environment.py
new file mode 100644
index 0000000..fbd3c8e
--- /dev/null
+++ b/mapgen/environment.py
@@ -0,0 +1,158 @@
+"""Stage `environment`: Holdridge life zone + seasonality tag + lithology + landform + ground."""
+from __future__ import annotations
+
+import numpy as np
+
+from .config import params
+from .crust import CRATON, LIP, POST_OROGEN, PRE_OROGEN, RIFT, SCAR, VOLCANO
+from .graph import nbr_max, nbr_min, steepest_receivers
+from . import minerals as MN
+from .noise import fbm
+
+REGIONS = ["polar", "subpolar", "boreal", "cool temperate", "warm temperate", "subtropical", "tropical"]
+ZONES = {
+    "polar": ["desert"],
+    "subpolar": ["dry tundra", "moist tundra", "wet tundra", "rain tundra"],
+    "boreal": ["desert", "dry scrub", "moist forest", "wet forest", "rain forest"],
+    "cool temperate": ["desert", "desert scrub", "steppe", "moist forest", "wet forest", "rain forest"],
+    "warm temperate": ["desert", "desert scrub", "thorn steppe", "dry forest", "moist forest", "wet forest",
+                       "rain forest"],
+    "subtropical": ["desert", "desert scrub", "thorn woodland", "dry forest", "moist forest", "wet forest",
+                    "rain forest"],
+    "tropical": ["desert", "desert scrub", "thorn woodland", "very dry forest", "dry forest", "moist forest",
+                 "wet forest", "rain forest"],
+}
+THRESH = {  # annual precipitation (mm) upper bounds between successive zones
+    "polar": [], "subpolar": [125, 250, 500], "boreal": [125, 250, 500, 1000],
+    "cool temperate": [125, 250, 500, 1000, 2000], "warm temperate": [125, 250, 500, 1000, 2000, 4000],
+    "subtropical": [125, 250, 500, 1000, 2000, 4000], "tropical": [125, 250, 500, 1000, 2000, 4000, 8000],
+}
+HOLDRIDGE_NAMES = [f"{r} {z}" for r in REGIONS for z in ZONES[r]]
+
+SEAS_F, SEAS_S, SEAS_W, SEAS_M = 0, 1, 2, 3
+SEASONALITY_NAMES = ["even rainfall", "dry summer", "dry winter / summer monsoon", "tropical monsoon"]
+
+LI_OCEAN, LI_GRANITE, LI_LIMESTONE, LI_BASALT, LI_ANDESITE, LI_SANDSTONE, LI_METAMORPHIC = range(7)
+LITHOLOGY_NAMES = ["oceanic basalt", "granite/gneiss", "limestone", "basalt", "andesite", "sandstone/shale",
+                   "metamorphic"]
+
+(LF_OCEAN, LF_PLAIN, LF_HILLS, LF_MOUNTAINS, LF_PLATEAU, LF_RIFT, LF_ESCARPMENT, LF_ARC, LF_MASSIF,
+ LF_BASALT_PLATEAU, LF_DUNES, LF_BADLANDS) = range(12)
+LANDFORM_NAMES = ["ocean", "plain", "hills", "mountains", "plateau", "rift valley", "escarpment",
+                  "volcanic arc", "volcanic massif", "basalt plateau", "dunes", "badlands"]
+
+(GR_NONE, GR_WETLAND, GR_BOG, GR_FLOODPLAIN, GR_DELTA, GR_MANGROVE, GR_SALT_FLAT, GR_PERMAFROST,
+ GR_KARST) = range(9)
+GROUND_NAMES = ["—", "wetland", "bog", "floodplain", "delta", "mangrove", "salt flat", "permafrost", "karst"]
+
+LEGENDS = {"deposits": MN.LEGEND, "holdridge": HOLDRIDGE_NAMES, "seasonality": SEASONALITY_NAMES, "landform": LANDFORM_NAMES,
+           "lithology": LITHOLOGY_NAMES, "ground": GROUND_NAMES}
+
+DEFAULTS = {"relief_km": 100.0, "coastal_km": 60.0, "dry_ratio_s": 3.0, "dry_ratio_w": 4.0, "monsoon_min_mm": 1500.0, "flat_m_per_km": 1.0,
+            "wet_min_mm": 700.0, "karst_min_mm": 800.0, "small_river_km3_yr": 5.0, "life": True}
+
+
+def holdridge(bio, p, tmin):
+    region = np.select([bio < 1.5, bio < 3, bio < 6, bio < 12, (bio < 24) & (tmin < 0), bio < 24],
+                       [0, 1, 2, 3, 4, 5], default=6).astype(np.int8)
+    zone = np.zeros(len(bio), np.int16)
+    offset = 0
+    for ri, r in enumerate(REGIONS):
+        sel = region == ri
+        zone[sel] = offset + np.searchsorted(THRESH[r], p[sel], side="right")
+        offset += len(ZONES[r])
+    return zone, region
+
+
+def seasonality(p_jun, p_dec, lat, tmin, p_ann, P):
+    summer = np.where(lat >= 0, p_jun, p_dec)
+    winter = np.where(lat >= 0, p_dec, p_jun)
+    tag = np.full(len(lat), SEAS_F, np.int8)
+    tag[summer < winter / P["dry_ratio_s"]] = SEAS_S
+    w = winter < summer / P["dry_ratio_w"]
+    tag[w] = SEAS_W
+    tag[w & (tmin >= 18.0) & (p_ann >= P["monsoon_min_mm"])] = SEAS_M
+    return tag
+
+
+def lithology(g, z, continental, age, d_over, seed):
+    nz = fbm(g.xyz, seed + 51, 4, 5.0)
+    lit = np.full(g.n, LI_GRANITE, np.int8)
+    lit[continental & (z < 600) & (nz > 0.1)] = LI_SANDSTONE
+    lit[continental & (z < 300) & (nz < -0.1)] = LI_LIMESTONE
+    lit[np.isin(age, [POST_OROGEN, PRE_OROGEN]) & (z > 2000)] = LI_METAMORPHIC
+    lit[~continental] = LI_OCEAN
+    lit[(d_over > 100) & (d_over < 320)] = LI_ANDESITE
+    lit[np.isin(age, [LIP, VOLCANO, SCAR])] = LI_BASALT
+    return lit
+
+
+def landform(g, z, age, lit, d_over, p_ann, ocean=None, relief_km=100.0):
+    hi, lo = np.asarray(z, dtype=np.float64), np.asarray(z, dtype=np.float64)
+    for _ in range(max(1, int(round(relief_km / g.spacing_km)))):   # window ≈ relief_km at any resolution
+        hi, lo = nbr_max(g, hi), nbr_min(g, lo)
+    relief = hi - lo
+    lf = np.full(g.n, LF_PLAIN, np.int8)
+    lf[relief >= 300] = LF_HILLS
+    lf[(z > 1000) & (relief < 500)] = LF_PLATEAU
+    lf[(relief >= 1500) | (z > 2500)] = LF_MOUNTAINS
+    lf[(p_ann < 250) & (relief < 300) & (lit == LI_SANDSTONE)] = LF_DUNES
+    lf[(p_ann >= 250) & (p_ann < 500) & (relief >= 200) & (relief < 800) & (lit == LI_SANDSTONE)] = LF_BADLANDS
+    lf[age == RIFT] = LF_RIFT
+    lf[age == LIP] = LF_BASALT_PLATEAU
+    ocean = z <= 0 if ocean is None else ocean
+    lf[(d_over > 100) & (d_over < 320) & ~ocean] = LF_ARC
+    lf[age == VOLCANO] = LF_MASSIF
+    lf[age == SCAR] = LF_ESCARPMENT
+    lf[ocean] = LF_OCEAN
+    return lf, relief
+
+
+def ground(g, z, lit, t_mean, t_min, p_ann, dist_ocean, river, strahler, discharge, lake, salt_flat, P, ocean=None):
+    land = ~(z <= 0 if ocean is None else ocean)
+    _, slope, _ = steepest_receivers(g, z)
+    flat = slope < P["flat_m_per_km"]
+    coastal = land & (dist_ocean < P["coastal_km"])
+    wet = land & flat & ~lake & (p_ann > P["wet_min_mm"]) & (discharge > P["small_river_km3_yr"])
+    gr = np.zeros(g.n, np.int8)
+    gr[land & (t_mean < -2.0)] = GR_PERMAFROST
+    gr[land & (lit == LI_LIMESTONE) & (p_ann > P["karst_min_mm"])] = GR_KARST
+    gr[wet & (t_mean < 5.0)] = GR_BOG
+    gr[wet & (t_mean >= 5.0)] = GR_WETLAND
+    gr[river & (strahler >= 3) & flat] = GR_FLOODPLAIN
+    gr[coastal & flat & (t_min >= 20.0) & (p_ann > 1000.0)] = GR_MANGROVE
+    gr[river & (strahler >= 4) & coastal] = GR_DELTA
+    gr[salt_flat] = GR_SALT_FLAT
+    if not P["life"]:
+        gr[np.isin(gr, [GR_WETLAND, GR_BOG, GR_MANGROVE])] = GR_NONE
+    return gr
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "environment", DEFAULTS)
+    (z, cont, age, d_over, t_mean, t_min, p_ann, p_jun, p_dec, bio, dist_ocean, river, strahler, q, lake,
+     salt) = ctx.need("elevation_eroded_m", "continental", "age_class", "d_over_km", "T_mean", "T_min", "P_ann",
+                      "P_jun", "P_dec", "biotemp", "dist_ocean_km", "river", "strahler", "discharge_km3_yr",
+                      "lake", "salt_flat")
+    z = z.astype(np.float64)
+    bed = z
+    if "lake_level_m" in ctx.data:                           # the ground's shape: lakes at their water, not their beds
+        lev = np.asarray(ctx.data["lake_level_m"], dtype=np.float64)
+        z = np.where(np.isfinite(lev), np.maximum(z, lev), z)
+    ocean = np.asarray(ctx.data["ocean"]) if "ocean" in ctx.data else z <= 0
+    zone, region = holdridge(bio, p_ann, t_min)
+    lit = lithology(g, z, cont, age, d_over, ctx.seed)
+    lf, relief = landform(g, z, age, lit, d_over, p_ann, ocean, P["relief_km"])
+    gr = ground(g, z, lit, t_mean, t_min, p_ann, dist_ocean, river, strahler, q, lake, salt, P, ocean)
+    coal = (lit == LI_SANDSTONE) & (p_ann > 600) & ~ocean & bool(P["life"])
+    get = lambda k, v: np.asarray(ctx.data[k]) if k in ctx.data else np.full(g.n, v)
+    deposits, main = MN.place(g, {"land": ~ocean & ~lake, "z": bed, "lit": lit, "age": age, "lf": lf, "relief": relief,
+                                  "t_mean": t_mean, "p_ann": p_ann, "d_coll": get("d_coll_km", np.inf),
+                                  "salt": salt, "endo": get("endorheic", False), "dist_ocean": dist_ocean,
+                                  "coal": coal}, ctx.seed)
+    iron = np.uint32(MN.BIT["bog_iron"] | MN.BIT["ironstone"] | MN.BIT["iron_high"])
+    return {"deposits": deposits, "deposit_main": main,"holdridge": zone, "hold_region": region,
+            "seasonality": seasonality(p_jun, p_dec, g.lat, t_min, p_ann, P),
+            "landform": lf, "relief_m": relief, "lithology": lit, "ground": gr,
+            "coal_potential": coal, "iron_potential": (deposits & iron) > 0}
diff --git a/mapgen/eras.py b/mapgen/eras.py
new file mode 100644
index 0000000..9b0dbc4
--- /dev/null
+++ b/mapgen/eras.py
@@ -0,0 +1,219 @@
+"""Eras: the base world plus local events. build_era applies an era's events
+to the built base, re-runs the stages they affect inside an influence mask (changed cells + mask_km) and keeps every
+cell outside it exactly as the base. Results: out/r<res>/eras/<era>/ (as out/r<res>/), previews in
+previews/r<res>/eras/<era>/. An era's cache key is the base world's inputs key plus its events (config/eras.toml). Each era's events happen at the end of the era before it; its `years` erode what the events so far
+changed; zone events follow the land at the start of their own era."""
+from __future__ import annotations
+
+import hashlib
+import importlib
+import json
+import re
+import shutil
+from pathlib import Path
+
+import numpy as np
+
+from . import config as C, environment as EN, events as EV, fields as FL, pipeline as P
+from .config import params
+from .erosion import DEFAULTS as EROSION_DEFAULTS, K_MULT, erode
+from .graph import distance_to, ocean_mask
+
+DEFAULTS = {"mask_km": 1500.0, "climate_blend_km": 300.0, "erosion_steps": 1, "years": 10000.0}
+RERUN = ("climate", "hydrology", "seabed", "environment", "ice", "fields")
+
+
+def era_dir(root: Path, res: int, name: str) -> Path:
+    return Path(root) / "out" / f"r{res}" / "eras" / name
+
+
+def steps(tect: dict, name: str) -> list:
+    """[(era, its own events, years)] for every era up to and including `name` ([eras] order): an era's events happen
+    at the end of the era before it; `years` is the time from them to the era's map (None: the default)."""
+    C.era_events(tect, name)                                  # an unknown era: ConfigError
+    eras = tect.get("eras") or {}
+    order = eras["order"]
+    by_name = {e["name"]: e for e in tect.get("event", [])}
+    return [(n, [by_name[e] for e in eras.get(n, {}).get("events", [])], eras.get(n, {}).get("years"))
+            for n in order[: order.index(name) + 1]]
+
+
+def fingerprint(base_key: str, steps_: list) -> str:
+    content = [[own, yrs] for _, own, yrs in steps_]
+    return hashlib.sha256((base_key + json.dumps(content, sort_keys=True)).encode()).hexdigest()[:16]
+
+
+def _same(a, b) -> bool:
+    """Equal arrays; NaN equals NaN in float arrays (a field may hold NaN where it has no value)."""
+    a, b = np.asarray(a), np.asarray(b)
+    return np.array_equal(a, b, equal_nan=a.dtype.kind in "fc" and b.dtype.kind in "fc")
+
+
+def merge_lake_ids(base, new, mask):
+    """Lake ids of an era: the base's outside the mask; inside it a lake the era left as it was keeps its base id and
+    every other lake gets a fresh id after the base's, so one id never names two lakes."""
+    base, new, mask = np.asarray(base), np.asarray(new), np.asarray(mask, bool)
+    out = base.copy()
+    nxt = max(int(base.max()) + 1, 0) if len(base) else 0
+    counts = np.bincount(base[base >= 0]) if (base >= 0).any() else np.zeros(0, np.int64)
+    ks = np.unique(new[mask & (new >= 0)])
+    idx = np.flatnonzero(np.isin(new, ks))
+    idx = idx[np.argsort(new[idx], kind="stable")]
+    starts = np.searchsorted(new[idx], ks)
+    ends = np.append(starts[1:], len(idx))
+    for s, e in zip(starts.tolist(), ends.tolist()):
+        cells = idx[s:e]
+        b = base[cells]
+        kept = b[0] >= 0 and bool(np.all(b == b[0])) and counts[b[0]] == len(cells)
+        into = cells[mask[cells]]
+        if kept:
+            out[into] = b[0]
+        else:
+            out[into] = nxt
+            nxt += 1
+    out[mask & (new < 0)] = -1
+    return out.astype(base.dtype)
+
+
+def zone_key(name: str) -> str:
+    return "zone_" + re.sub(r"[^A-Za-z0-9]", "_", name)
+
+
+def _heights(g, z, events, min_sea, log, name, uncut=lambda z: z):
+    for ev in events:                                        # heights first, in era order
+        if ev["kind"] == "disintegrate":
+            z, z_ref = EV.disintegrate(g.xyz, z, ~ocean_mask(g, uncut(z), min_sea), ev, g.radius_km)
+            if z_ref is None:
+                log(f"era {name}: warning: {ev['name']} has no land in its rim ring (0.9–1.0 × radius): "
+                    f"its bowl hangs from 0 m")
+            else:
+                log(f"era {name}: {ev['name']} (ground at its rim {z_ref:.0f} m)")
+        elif ev["kind"] == "volcano":
+            z = EV.volcano(g.xyz, z, ev, g.radius_km)
+    return z
+
+
+def build_era(root: Path, res: int, name: str, log=print, base=None, low_memory: bool = False) -> Path:
+    root = Path(root)
+    cfg, tect = C.load(root)
+    events = C.era_events(tect, name)
+    if not events:
+        log(f"era {name}: the base world (out/r{res})")
+        return root / "out" / f"r{res}"
+    out = era_dir(root, res, name)
+    base_key = P.inputs_key(root, res)
+    plan = steps(tect, name)
+    key = fingerprint(base_key, plan)
+    label = tect["eras"][name].get("label", name)
+    meta_path = out / "cells_meta.json"
+    if meta_path.exists() and (out / "cells.npz").exists():
+        meta = json.loads(meta_path.read_text())
+        era = meta.get("era", {})
+        if era.get("fingerprint") == key:
+            if (era.get("label"), era.get("name")) != (label, name):   # a new label or name needs no rebuild
+                era["label"], era["name"] = label, name
+                meta_path.write_text(json.dumps(meta, indent=1))
+            log(f"era {name}: cached")
+            return out
+    ctx = base if base is not None else P.build(root, res, stop="fields", log=lambda m: None, low_memory=low_memory)
+    low_memory = low_memory or ctx.low_memory
+    g, d = ctx.grid, ctx.data
+    E = params(ctx.cfg, "eras", DEFAULTS)
+    EP = {**params(ctx.cfg, "erosion", EROSION_DEFAULTS), "steps": E["erosion_steps"]}
+    kmult = np.vectorize(K_MULT.get)(np.asarray(d["age_class"])).astype(np.float64)
+    z0 = np.asarray(d["elevation_eroded_m"], dtype=np.float64)
+    z, ocean = z0.copy(), np.asarray(d["ocean"])
+    water = np.asarray(d.get("open_water", ocean))
+    cut = np.asarray(d.get("lake_cut_m", np.zeros(g.n)), dtype=np.float64)
+    uncut = lambda zz: np.where(zz == z0, zz + cut, zz)       # sea masks: the ground before lake beds were carved
+    zones, lands = [], []
+    for _, own, yrs in plan:                                 # era by era: its events at the end of the era before
+        for ev in own:                                       # zones follow the land before their own era's events
+            if ev["kind"] == "zone":
+                zones.append(ev)
+                lands.append(EV.landmass(g, ~water, ev["seed"], ev["name"]) if ev["shape"] == "landmass" else None)
+        z = _heights(g, z, own, EP["min_sea_km2"], log, name, uncut)
+        hit = z != z0
+        if hit.any():                                        # the changed ground weathers for the era's years
+            years = float(E["years"] if yrs is None else yrs)
+            ze = erode(g, z, kmult, {**EP, "dt_myr": years / 1.0e6 / E["erosion_steps"]})
+            z = np.where(hit, np.maximum(np.minimum(z, ze), -11000.0), z)
+            ocean, water = ocean_mask(g, uncut(z), np.inf), ocean_mask(g, uncut(z), EP["min_sea_km2"])
+    hit = z != z0
+    weights = [EV.zone_weight(g.xyz, ev, g.radius_km, ctx.seed, None if land is None else g.xyz[land])
+               for ev, land in zip(zones, lands)]
+    changed = hit | (ocean != np.asarray(d["ocean"])) | (water != np.asarray(d.get("open_water", d["ocean"])))
+    for w in weights:
+        changed |= w > 0
+    mask = distance_to(g, changed) <= E["mask_km"] if changed.any() else np.zeros(g.n, bool)
+
+    ec = P.Ctx(ctx.root, ctx.cfg, ctx.tect, ctx.res, dict(d), g, low_memory=low_memory)
+    ec.data.update({"elevation_eroded_m": z.astype(np.asarray(d["elevation_eroded_m"]).dtype), "ocean": ocean,
+                    "open_water": water, "lake_carve": hit})   # lake beds: geology, carved where the ground moved
+    new_keys = {}
+    for s in RERUN:
+        o = importlib.import_module(f"mapgen.{s}").run(ec)
+        ec.data.update(o)
+        new_keys[s] = list(o)
+    blend = np.clip(distance_to(g, ~mask) / E["climate_blend_km"], 0.0, 1.0) if (~mask).any() else np.ones(g.n)
+    shape = lambda a, v: v.reshape(-1, *([1] * (a.ndim - 1)))
+    merged = dict(d)
+    for s, keys in new_keys.items():
+        for k in keys:
+            nv, bv = np.asarray(ec.data[k]), np.asarray(d[k])
+            if s == "climate" and nv.dtype.kind == "f":      # solved on the whole world, blended in over 300 km
+                nv = (bv + shape(nv, blend) * (nv - bv)).astype(nv.dtype)
+            merged[k] = np.where(shape(nv, mask), nv, bv)
+    for k in ("deposits", "deposit_main", "iron_potential"):  # ore is the base's geology: events move ground, the
+        if k in d:                                           # rank-by-share placement would shift it world-wide
+            merged[k] = d[k]
+    if "lake_id" in merged:                                  # one id never names two lakes
+        merged["lake_id"] = merge_lake_ids(d["lake_id"], ec.data["lake_id"], mask)
+    merged["elevation_eroded_m"] = np.where(mask, ec.data["elevation_eroded_m"], d["elevation_eroded_m"])   # events and
+    merged["ocean"] = ocean                                                                  # lake beds: inside only
+    if "open_water" in d:
+        merged["open_water"] = water
+    H = float(ctx.cfg["planet"]["scale_height_km"]) * 1000.0
+    for ev, w, land in zip(zones, weights, lands):          # zone events write their fields last
+        merged.update(EV.apply_zone(merged, w, ev, merged["z_surface_m"], H))
+        merged[zone_key(ev["name"])] = w.astype(np.float32)
+        if land is not None:
+            merged[zone_key(ev["name"]) + "_land"] = land
+    if zones:                                                # plants settle into the zone's air and gravity
+        pl = ctx.cfg["planet"]
+        merged["plant_height_x"] = FL.plant_height(pl["gravity_g"], merged["gravity_g"]).astype(
+            np.asarray(d["plant_height_x"]).dtype)
+        sea_p = np.asarray(merged["pressure_bar"], dtype=np.float64) / FL.pressure(1.0, merged["z_surface_m"], H)
+        wet = np.sqrt(sea_p / pl["sea_level_pressure_bar"])  # more CO₂ per breath: less water lost per growth
+        hz, hr = EN.holdridge(merged["biotemp"], np.asarray(merged["P_ann"], dtype=np.float64) * wet, merged["T_min"])
+        inside = np.logical_or.reduce([w > 0 for w in weights])
+        merged["holdridge"] = np.where(inside, hz, merged["holdridge"]).astype(np.asarray(d["holdridge"]).dtype)
+        merged["hold_region"] = np.where(inside, hr, merged["hold_region"]).astype(np.asarray(d["hold_region"]).dtype)
+    merged["era_mask"] = mask
+    for k, a in d.items():                                   # guard: nothing outside the mask may differ
+        a, b = np.asarray(a), np.asarray(merged[k])
+        if a.shape[:1] == (g.n,) and not _same(a[~mask], b[~mask]):
+            raise RuntimeError(f"era {name}: {k} changed outside the influence mask")
+
+    rc = P.Ctx(ctx.root, ctx.cfg, ctx.tect, ctx.res, merged, g, out_dir=out,
+               preview_dir=root / "previews" / f"r{res}" / "eras" / name, low_memory=low_memory)
+    importlib.import_module("mapgen.render").run(rc)
+    meta = json.loads(meta_path.read_text())
+    meta["era"] = {"name": name, "label": label, "events": events, "fingerprint": key, "base_key": base_key,
+                   "steps": [[n, [e["name"] for e in own], yrs] for n, own, yrs in plan],
+                   "mask_km": E["mask_km"], "mask_cells": int(mask.sum())}
+    meta_path.write_text(json.dumps(meta, indent=1))
+    log(f"era {name}: {int(mask.sum())} of {g.n} cells inside the influence mask")
+    return out
+
+
+def build_all(root: Path, res: int, log=print, base=None, low_memory: bool = False) -> list:
+    """Every configured era with events (the base era is the base build); era folders no longer configured go."""
+    _, tect = C.load(Path(root))
+    names = [n for n in (tect.get("eras") or {}).get("order", []) if C.era_events(tect, n)]
+    built = [build_era(root, res, n, log, base, low_memory) for n in names]
+    top = Path(root) / "out" / f"r{res}" / "eras"
+    for p in (top.iterdir() if top.is_dir() else []):
+        if p.is_dir() and p.name not in names:
+            shutil.rmtree(p, ignore_errors=True)
+    return built
diff --git a/mapgen/erosion.py b/mapgen/erosion.py
new file mode 100644
index 0000000..5384988
--- /dev/null
+++ b/mapgen/erosion.py
@@ -0,0 +1,54 @@
+"""Stage `erosion`: implicit stream-power incision (Braun & Willett 2013) + km-scale hillslope smoothing."""
+from __future__ import annotations
+
+import numpy as np
+
+from .config import params
+from .crust import CRATON, LIP, OCEANIC, POST_OROGEN, PRE_OROGEN, RIFT, SCAR, VOLCANO
+from .elevation import solve_sea_level_connected
+from .graph import (accumulate, drop_smooth_cache, ocean_mask, priority_flood, receiver_levels, smooth_km,
+                    steepest_receivers)
+
+DEFAULTS = {"min_sea_km2": 5.0e6, "steps": 12, "dt_myr": 5.0, "k": 0.02, "m": 0.5, "hillslope_km": 25.0, "sediment_fill": 0.15}
+K_MULT = {OCEANIC: 0.0, CRATON: 1.5, PRE_OROGEN: 1.2, POST_OROGEN: 1.0, RIFT: 1.0, LIP: 0.6, SCAR: 0.9, VOLCANO: 0.8}
+
+
+def erode(g, z, kmult, P):
+    z = np.asarray(z, dtype=np.float64).copy()
+    ar = np.arange(g.n)
+    for _ in range(int(P["steps"])):
+        ocean = ocean_mask(g, z, P["min_sea_km2"])
+        zb = np.where(ocean, 0.0, z)                     # base level = sea level, not the seafloor
+        zf = priority_flood(g, zb, ocean)
+        zb = np.where(ocean, zb, zb + P["sediment_fill"] * (zf - zb))   # sediment infills closed basins
+        recv, _, dist = steepest_receivers(g, zf)
+        recv = np.where(ocean, ar, recv)
+        levels = receiver_levels(recv)
+        area = accumulate(recv, levels, g.area_km2)
+        F = P["k"] * kmult * P["dt_myr"] * area ** P["m"] / np.where(np.isfinite(dist), dist, 1.0)
+        F[ocean] = 0.0
+        znew = zb.copy()
+        for lv in levels[1:]:
+            r = recv[lv]
+            znew[lv] = np.minimum(zb[lv], (zb[lv] + F[lv] * znew[r]) / (1.0 + F[lv]))
+        znew = smooth_km(g, znew, P["hillslope_km"], keep=True)   # hillslope/sub-grid smoothing, fixed length in km
+        z = np.where(ocean, z, znew)
+    drop_smooth_cache(g)
+    return z
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "erosion", DEFAULTS)
+    z0, age = ctx.need("elevation_m", "age_class")
+    kmult = np.vectorize(K_MULT.get)(age).astype(np.float64)
+    z0 = np.asarray(z0, dtype=np.float64)
+    added = np.asarray(ctx.data.get("land_added", np.zeros(g.n, bool)), bool)   # land the sea-level solve didn't
+    target = ctx.cfg["build"]["land_fraction"] + g.area_km2[added & (z0 > 0)].sum() / g.area_km2.sum()   # see
+    z = erode(g, z0, kmult, P)
+    z = solve_sea_level_connected(g, z, target, P["min_sea_km2"])
+    z = np.maximum(z, -11000.0)
+    # the sea is the connected world ocean; a separate basin ≥ min_sea_km2 shaped the coasts above as sea and stays
+    # open water for the climate (open_water), but its water is a lake (hydrology), not sea
+    return {"elevation_eroded_m": z.astype(np.float32), "ocean": ocean_mask(g, z, np.inf),
+            "open_water": ocean_mask(g, z, P["min_sea_km2"])}
diff --git a/mapgen/events.py b/mapgen/events.py
new file mode 100644
index 0000000..cba874e
--- /dev/null
+++ b/mapgen/events.py
@@ -0,0 +1,87 @@
+"""Era events: pure functions of positions and the event's config, so the era
+builder (world cells) and the refinement (fine cells) apply the same event."""
+from __future__ import annotations
+
+import numpy as np
+from scipy.spatial import cKDTree
+
+from .fields import pressure
+from .graph import components
+from .noise import fbm, name_seed
+from .pipeline import StageError
+from .sphere import gc_dist_km, latlon_to_xyz
+from .zones import smootherstep
+
+
+def disintegrate(xyz, z, land, ev: dict, radius_km: float):
+    """(new heights, z_ref): inside radius_km of the centre min(old, z_ref − depth · sqrt(1 − (d/r)²)); z_ref = mean
+    ground height of the land in the ring 0.9–1.0 r before the event, None when no land lies in that ring (the bowl
+    then hangs from 0 m). Material vanishes: no ejecta, no rim."""
+    z = np.asarray(z, dtype=np.float64)
+    d = gc_dist_km(np.asarray(xyz, dtype=np.float64), latlon_to_xyz(*ev["center"]), radius_km)
+    r = float(ev["radius_km"])
+    ring = np.asarray(land, bool) & (d >= 0.9 * r) & (d <= r)
+    z_ref = float(np.mean(z[ring])) if ring.any() else None   # None: no land on its rim (the bowl hangs from 0 m)
+    bowl = (0.0 if z_ref is None else z_ref) - ev["depth_m"] * np.sqrt(np.clip(1.0 - (d / r) ** 2, 0.0, 1.0))
+    return np.where(d < r, np.minimum(z, bowl), z), z_ref
+
+
+def volcano(xyz, z, ev: dict, radius_km: float):
+    """Heights after a volcano (only raises): cone, shield or caldera; peak_m above the old ground (peak_mode "above")
+    or as an absolute height ("absolute")."""
+    z = np.asarray(z, dtype=np.float64)
+    d = gc_dist_km(np.asarray(xyz, dtype=np.float64), latlon_to_xyz(*ev["center"]), radius_km)
+    x = np.clip(d / ev["radius_km"], 0.0, 1.0)
+    prof = {"cone": (1.0 - x) ** 1.5, "shield": (1.0 - x * x) ** 2,
+            "caldera": np.maximum((1.0 - x) ** 1.5 - 0.6 * np.exp(-(x / 0.15) ** 2), 0.0)}[ev.get("shape", "cone")]
+    peak = float(ev["peak_m"])
+    new = z + peak * prof if ev.get("peak_mode", "above") == "above" else z + np.maximum(peak - z, 0.0) * prof
+    return np.maximum(z, new)
+
+
+def landmass(g, land, seed_latlon, event_name: str):
+    """The connected land (grid cells) containing the seed point."""
+    lab = components(g, land)
+    i = g.cell_index(*seed_latlon)
+    if lab[i] < 0:
+        land = np.flatnonzero(np.asarray(land, bool))
+        hint = ""
+        if len(land):
+            j = land[np.argmax(g.xyz[land] @ g.xyz[i])]
+            hint = f"; the nearest land is at [{g.lat[j]:.1f}, {g.lon[j]:.1f}]"
+        raise StageError(f"event {event_name}: its seed {list(seed_latlon)} lies in the sea in this era{hint} "
+                         f"(move the seed in config/eras.toml)")
+    return lab == lab[i]
+
+
+def zone_weight(xyz, ev: dict, radius_km: float, seed: int, land_xyz=None):
+    """0–1 per point: 1 inside, 0 outside; a linear ramp edge_km wide at the edge (sharp), or the smooth profile.
+    circle: centre + radius_km. landmass: within reach of the landmass cells land_xyz (unit vectors), the reach
+    varying between reach_km by low-frequency noise."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    if ev["shape"] == "circle":
+        d = gc_dist_km(xyz, latlon_to_xyz(*ev["center"]), radius_km)
+        reach = np.full(len(xyz), float(ev["radius_km"]))
+    else:
+        chord, _ = cKDTree(np.asarray(land_xyz, dtype=np.float64)).query(xyz)
+        d = 2.0 * np.arcsin(np.clip(chord / 2.0, 0.0, 1.0)) * radius_km
+        lo, hi = ev["reach_km"]
+        reach = lo + (hi - lo) * np.clip(0.5 + 1.5 * fbm(xyz, name_seed(seed, ev["name"]), 3, 1.5), 0.0, 1.0)
+    if "edge_km" in ev:
+        return np.clip((reach - d) / ev["edge_km"], 0.0, 1.0)
+    return 1.0 - smootherstep(d / reach)
+
+
+def apply_zone(fields: dict, w, ev: dict, z_surface_m, scale_height_m: float) -> dict:
+    """Fields inside a zone event: absolute values blended by w (gravity_g, o2_fraction, fire_reactivity; pressure_bar
+    is the sea-level value, falling with height as elsewhere); po2_bar follows. Never height. Cells with w = 0 keep
+    their values exactly."""
+    w = np.asarray(w, dtype=np.float64)
+    decay = pressure(1.0, z_surface_m, scale_height_m)
+    new = {}
+    for k, v in ev["fields"].items():
+        base = np.asarray(fields[k], dtype=np.float64)
+        new[k] = (base / decay * (1.0 - w) + v * w) * decay if k == "pressure_bar" else base * (1.0 - w) + v * w
+    o2 = new.get("o2_fraction", np.asarray(fields["o2_fraction"], dtype=np.float64))
+    new["po2_bar"] = o2 * new.get("pressure_bar", np.asarray(fields["pressure_bar"], dtype=np.float64))
+    return {k: np.where(w > 0, v, fields[k]).astype(np.asarray(fields[k]).dtype) for k, v in new.items()}
diff --git a/mapgen/fields.py b/mapgen/fields.py
new file mode 100644
index 0000000..bf83d3f
--- /dev/null
+++ b/mapgen/fields.py
@@ -0,0 +1,36 @@
+"""Stage `fields` and shared zone modifiers."""
+from __future__ import annotations
+
+import numpy as np
+
+from .config import MASK_DEFAULTS, params
+
+
+def gravity_mod(cfg, m_gravity):
+    M = params(cfg, "masks", MASK_DEFAULTS)
+    return np.clip(1.0 + M["gravity_range"] * np.asarray(m_gravity, dtype=np.float64), 0.05, None)
+
+
+def o2_mod(cfg, m_o2):
+    M = params(cfg, "masks", MASK_DEFAULTS)
+    return np.clip(1.0 + M["o2_range"] * np.asarray(m_o2, dtype=np.float64), 0.0, None)
+
+
+def pressure(p0_bar, z_m, scale_height_m):
+    """Air pressure at ground height z (bar): sea-level pressure p0 falling with height (below 0 m: p0)."""
+    return p0_bar * np.exp(-np.maximum(np.asarray(z_m, dtype=np.float64), 0.0) / scale_height_m)
+
+
+def plant_height(g_planet, gravity_g):
+    """How much taller plants grow than under the planet's normal gravity (height limit ∝ 1 / g)."""
+    return g_planet / np.asarray(gravity_g, dtype=np.float64)
+
+
+def run(ctx) -> dict:
+    pl = ctx.cfg["planet"]
+    zs, m_o2, m_g = ctx.need("z_surface_m", "m_o2_zones", "m_gravity_zones")
+    p = pressure(pl["sea_level_pressure_bar"], zs, pl["scale_height_km"] * 1000.0)
+    o2 = pl["o2_fraction"] * o2_mod(ctx.cfg, m_o2)
+    return {"po2_bar": o2 * p, "gravity_g": pl["gravity_g"] * gravity_mod(ctx.cfg, m_g), "pressure_bar": p,
+            "o2_fraction": o2, "fire_reactivity": np.ones(len(p)),   # zone events (eras) set other values
+            "plant_height_x": plant_height(pl["gravity_g"], pl["gravity_g"] * gravity_mod(ctx.cfg, m_g))}
diff --git a/mapgen/geo.py b/mapgen/geo.py
new file mode 100644
index 0000000..152ae66
--- /dev/null
+++ b/mapgen/geo.py
@@ -0,0 +1,200 @@
+"""GeoJSON exports: rivers, coastlines/lakes (marching squares), plate boundaries."""
+from __future__ import annotations
+
+import json
+from pathlib import Path
+
+import h3.api.basic_int as h3
+import numpy as np
+
+BOUNDARY_NAMES = {1: "convergent", 2: "divergent", 3: "transform"}
+# marching-squares cases; corner bits TL=8 TR=4 BR=2 BL=1; edges T R B L
+_CASES = {1: [("L", "B")], 2: [("B", "R")], 3: [("L", "R")], 4: [("T", "R")], 5: [("T", "R"), ("L", "B")],
+          6: [("T", "B")], 7: [("L", "T")], 8: [("L", "T")], 9: [("T", "B")], 10: [("L", "T"), ("B", "R")],
+          11: [("T", "R")], 12: [("L", "R")], 13: [("B", "R")], 14: [("L", "B")]}
+_EDGE = {"T": (1, 0), "R": (2, 1), "B": (1, 2), "L": (0, 1)}   # doubled (dx, dy) inside a 2x2 block
+
+
+def split_antimeridian(coords):
+    parts, cur = [], [list(coords[0])]
+    for a, b in zip(coords, coords[1:]):
+        if abs(b[0] - a[0]) > 180.0:
+            parts.append(cur)
+            cur = []
+        cur.append(list(b))
+    parts.append(cur)
+    return [p for p in parts if len(p) >= 2]
+
+
+def _feature(geom_type, coords, props):
+    return {"type": "Feature", "geometry": {"type": geom_type, "coordinates": coords}, "properties": props}
+
+
+def river_lines(g, recv, river, strahler, discharge):
+    """One LineString per run of equal Strahler order (runs share their junction vertex)."""
+    n = g.n
+    ar = np.arange(n)
+    has_donor = np.zeros(n, bool)
+    nr = river & (recv != ar)
+    has_donor[recv[nr]] = True
+    visited = np.zeros(n, bool)
+    feats = []
+    for s in np.flatnonzero(river & ~has_donor):
+        cells, c = [int(s)], int(s)
+        while True:
+            visited[c] = True
+            r = int(recv[c])
+            if r == c:
+                break
+            cells.append(r)
+            if not river[r] or visited[r]:
+                break
+            c = r
+        body = cells[:-1] if len(cells) > 1 else cells
+        start = 0
+        for i in range(1, len(body) + 1):
+            if i == len(body) or strahler[body[i]] != strahler[body[i - 1]]:
+                run = body[start:i] + [cells[i] if i < len(cells) else body[-1]]
+                coords = [[round(float(g.lon[k]), 4), round(float(g.lat[k]), 4)] for k in run]
+                props = {"order": int(strahler[body[start]]), "discharge_km3_yr": round(float(discharge[body[i - 1]]), 3)}
+                feats += [_feature("LineString", p, props) for p in split_antimeridian(coords)]
+                start = i
+    return feats
+
+
+def boundary_features(g, plate, bnd_type):
+    s, d = g.src, g.dst
+    e = np.flatnonzero((plate[s] != plate[d]) & (s < d))
+    lines = {}
+    for k in e:
+        a, b = int(g.ids[s[k]]), int(g.ids[d[k]])
+        seg = [[round(lng, 4), round(lat, 4)] for lat, lng in h3.directed_edge_to_boundary(h3.cells_to_directed_edge(a, b))]
+        if abs(seg[0][0] - seg[-1][0]) > 180:
+            continue
+        lines.setdefault(BOUNDARY_NAMES.get(int(bnd_type[s[k]]), "transform"), []).append(seg)
+    return [_feature("MultiLineString", v, {"type": t}) for t, v in sorted(lines.items())]
+
+
+def contours(mask):
+    m = np.pad(np.asarray(mask, dtype=np.uint8), 1)
+    code = m[:-1, :-1] * 8 + m[:-1, 1:] * 4 + m[1:, 1:] * 2 + m[1:, :-1]
+    adj = {}
+    for case, segs in _CASES.items():
+        ys, xs = np.nonzero(code == case)
+        for (e1, e2) in segs:
+            for y, x in zip(ys.tolist(), xs.tolist()):
+                p = (2 * x + _EDGE[e1][0], 2 * y + _EDGE[e1][1])
+                q = (2 * x + _EDGE[e2][0], 2 * y + _EDGE[e2][1])
+                adj.setdefault(p, []).append(q)
+                adj.setdefault(q, []).append(p)
+    seen, rings = set(), []
+    for start in adj:
+        if start in seen:
+            continue
+        ring, prev, cur = [start], None, start
+        seen.add(start)
+        while True:
+            nxt = [q for q in adj[cur] if q != prev]
+            if not nxt:
+                break
+            prev, cur = cur, nxt[0]
+            ring.append(cur)
+            if cur == start:
+                break
+            seen.add(cur)
+        rings.append([(px / 2.0 - 1.0, py / 2.0 - 1.0) for px, py in ring])
+    return rings
+
+
+def _signed_area(ring):
+    a = np.asarray(ring, dtype=np.float64)
+    return 0.5 * float(np.sum(a[:-1, 0] * a[1:, 1] - a[1:, 0] * a[:-1, 1]))
+
+
+def _point_in_ring(pt, ring):
+    a = np.asarray(ring, dtype=np.float64)
+    x1, y1, x2, y2 = a[:-1, 0], a[:-1, 1], a[1:, 0], a[1:, 1]
+    cross = (y1 > pt[1]) != (y2 > pt[1])
+    xi = x1 + (pt[1] - y1) * (x2 - x1) / np.where(y2 != y1, y2 - y1, 1e-300)
+    return bool(np.count_nonzero(cross & (pt[0] < xi)) % 2)
+
+
+def _orient(ring, ccw):
+    return ring if (_signed_area(ring) > 0) == ccw else ring[::-1]
+
+
+def _clip(ring, left, x0=180.0):
+    """Sutherland–Hodgman clip of a closed ring to x ≤ x0 (left) or x ≥ x0."""
+    inside = (lambda p: p[0] <= x0) if left else (lambda p: p[0] >= x0)
+    cut = lambda p, q: [x0, p[1] + (x0 - p[0]) / (q[0] - p[0]) * (q[1] - p[1])]
+    pts, out = ring[:-1], []
+    for i in range(len(pts)):
+        cur, prev = pts[i], pts[i - 1]
+        if inside(cur):
+            if not inside(prev):
+                out.append(cut(prev, cur))
+            out.append(cur)
+        elif inside(prev):
+            out.append(cut(prev, cur))
+    return out + [out[0]] if len(out) >= 3 else []
+
+
+def _split_seam(rings):
+    """rings[0] exterior + holes in continuous longitude; split at +180 into ≤2 polygons in [−180, 180]."""
+    if max(p[0] for p in rings[0]) <= 180.0:
+        return [rings]
+    polys = []
+    for left in (True, False):
+        ext = _clip(rings[0], left)
+        if len(ext) < 4:
+            continue
+        holes = [h for h in (_clip(r, left) for r in rings[1:]) if len(h) >= 4]
+        part = [ext] + holes
+        if not left:
+            part = [[[x - 360.0, y] for x, y in r] for r in part]
+        polys.append(part)
+    return polys
+
+
+def _round(poly):
+    return [[[round(x, 4), round(y, 4)] for x, y in r] for r in poly]
+
+
+def contour_features(mask, kind):
+    """Land/lake outlines as RFC 7946 (Multi)Polygons: CCW exteriors, CW holes, split at the antimeridian."""
+    H, W = mask.shape
+    col = int(np.argmin(mask.sum(axis=0)))
+    rings = []
+    for ring in contours(np.roll(mask, -col, axis=1)):
+        pts = [[(x + col + 0.5) / W * 360.0 - 180.0, max(-90.0, min(90.0, 90.0 - (y + 0.5) / H * 180.0))]
+               for x, y in ring]
+        if len(pts) >= 4 and pts[0] == pts[-1]:
+            rings.append(pts)
+    depth = [sum(_point_in_ring(r[0], o) for j, o in enumerate(rings) if j != i) for i, r in enumerate(rings)]
+    feats = []
+    for i, r in enumerate(rings):
+        if depth[i] % 2:
+            continue
+        holes = [_orient(rings[j], False) for j in range(len(rings))
+                 if depth[j] == depth[i] + 1 and _point_in_ring(rings[j][0], r)]
+        polys = [_round(p) for p in _split_seam([_orient(r, True)] + holes)]
+        if len(polys) == 1:
+            feats.append(_feature("Polygon", polys[0], {"kind": kind}))
+        elif polys:
+            feats.append(_feature("MultiPolygon", polys, {"kind": kind}))
+    return feats
+
+
+def _write(path: Path, feats) -> None:
+    path.write_text(json.dumps({"type": "FeatureCollection", "features": feats}, separators=(",", ":")))
+
+
+def write_all(ctx, out_dir: Path, land_raster, lake_raster) -> None:
+    g, d = ctx.grid, ctx.data
+    gdir = out_dir / "geo"
+    gdir.mkdir(parents=True, exist_ok=True)
+    _write(gdir / "rivers.geojson", river_lines(g, d["recv"], d["river"], d["strahler"], d["discharge_km3_yr"]))
+    _write(gdir / "plate_boundaries.geojson", boundary_features(g, d["plate"], d["bnd_type"]))
+    step = max(1, land_raster.shape[1] // 2048)
+    _write(gdir / "coast.geojson", contour_features(land_raster[::step, ::step], "land"))
+    _write(gdir / "lakes.geojson", contour_features(lake_raster[::step, ::step] & land_raster[::step, ::step], "lake"))
diff --git a/mapgen/graph.py b/mapgen/graph.py
new file mode 100644
index 0000000..126c71f
--- /dev/null
+++ b/mapgen/graph.py
@@ -0,0 +1,477 @@
+"""Operations on the H3 cell graph (CSR neighbours)."""
+from __future__ import annotations
+
+import heapq
+
+import numpy as np
+from scipy import sparse
+from scipy.sparse import csgraph
+from scipy.sparse import linalg as splinalg
+
+
+def nbr_mean(g, f):
+    return np.bincount(g.src, weights=f[g.dst], minlength=g.n) / g.counts
+
+
+def nbr_max(g, f):
+    out = np.asarray(f, dtype=np.float64).copy()
+    np.maximum.at(out, g.src, f[g.dst])
+    return out
+
+
+def nbr_min(g, f):
+    out = np.asarray(f, dtype=np.float64).copy()
+    np.minimum.at(out, g.src, f[g.dst])
+    return out
+
+
+def diffuse(g, f, iters: int, alpha: float = 0.5, mask=None):
+    f = np.asarray(f, dtype=np.float64)
+    for _ in range(int(iters)):
+        new = (1.0 - alpha) * f + alpha * nbr_mean(g, f)
+        f = new if mask is None else np.where(mask, new, f)
+    return f
+
+
+def laplacian_matrix(g):
+    """Graph Laplacian (per km²): Δf_i ≈ (4/k_i) Σ_j (f_j − f_i)/d_ij² (exact for a regular hex lattice)."""
+    w = 4.0 / (g.counts[g.src] * g.edge_km**2)
+    lap = sparse.csr_matrix((w, (g.src, g.dst)), shape=(g.n, g.n))
+    return lap - sparse.diags(np.asarray(lap.sum(axis=1)).ravel())
+
+
+def workers() -> int:
+    """Threads for independent jobs (seasons, components): WORLDGEN_THREADS, default 3. Each job computes exactly
+    what it would alone, so results never depend on it; numpy and scipy's sparse kernels release the GIL."""
+    import os
+    try:
+        return max(1, int(os.environ.get("WORLDGEN_THREADS", "3")))
+    except ValueError:
+        return 3
+
+
+def pmap(fn, items) -> list:
+    """[fn(x) for x in items] on up to workers() threads, in order."""
+    items = list(items)
+    if workers() <= 1 or len(items) <= 1:
+        return [fn(x) for x in items]
+    from concurrent.futures import ThreadPoolExecutor
+    with ThreadPoolExecutor(min(workers(), len(items))) as ex:
+        return list(ex.map(fn, items))
+
+
+def smooth_km(g, f, length_km: float, rtol: float = 1e-6, keep: bool = False):
+    """Resolution-independent smoothing: solve (I − L²Δ) s = f (screened Poisson, decay length ≈ L km).
+    keep: keep the matrix on g for the next call (loops over one grid; drop_smooth_cache(g) frees it)."""
+    f = np.asarray(f, dtype=np.float64)
+    scale = float(np.max(np.abs(f))) if f.size else 0.0
+    if length_km <= 0 or scale == 0.0:
+        return f.copy()
+    f = f / scale                                    # linear system: solve at unit scale (avoids breakdown)
+    cache = g.__dict__.get("_screened", {})
+    if length_km in cache:
+        A, inv_diag = cache[length_km]
+    else:
+        A = (sparse.identity(g.n, format="csr") - length_km**2 * laplacian_matrix(g)).tocsr()
+        inv_diag = 1.0 / A.diagonal()
+        if keep:
+            g.__dict__.setdefault("_screened", {})[length_km] = (A, inv_diag)
+    s, info = bicgstab_jacobi(A, f, f.copy(), inv_diag, rtol, 5000)
+    if info != 0:
+        raise ValueError(f"smooth_km: solver did not converge (info={info})")
+    return s * scale
+
+
+def _make_bicg_jit():
+    """bicgstab's vector updates fused into single passes (numba optional; WORLDGEN_NO_JIT=1 turns it off). Each
+    element gets the same operations in the same order as scipy's numpy statements; dot products, norms and sparse
+    products stay the very calls scipy makes, so the iterates — and the answer — are the same floats. (Threads were
+    tried and dropped: on a busy machine they wait more than they work.)"""
+    import os
+    if os.environ.get("WORLDGEN_NO_JIT"):
+        return None
+    try:
+        import numba
+    except ImportError:
+        return None
+
+    @numba.njit(cache=True, nogil=True)
+    def p_update(p, v, r, omega, beta):          # p -= omega*v; p *= beta; p += r
+        for i in range(len(p)):
+            p[i] = (p[i] - omega * v[i]) * beta + r[i]
+
+    @numba.njit(cache=True, nogil=True)
+    def scale(out, d, x):                        # out = d * x (the Jacobi preconditioner)
+        for i in range(len(x)):
+            out[i] = d[i] * x[i]
+
+    @numba.njit(cache=True, nogil=True)
+    def axpy_neg(r, a, v):                       # r -= a*v
+        for i in range(len(r)):
+            r[i] = r[i] - a * v[i]
+
+    @numba.njit(cache=True, nogil=True)
+    def x_update(x, alpha, phat, omega, shat):   # x += alpha*phat; x += omega*shat
+        for i in range(len(x)):
+            x[i] = (x[i] + alpha * phat[i]) + omega * shat[i]
+
+    @numba.njit(cache=True, nogil=True)
+    def axpy(x, a, v):                           # x += a*v
+        for i in range(len(x)):
+            x[i] = x[i] + a * v[i]
+
+    return p_update, scale, axpy_neg, x_update, axpy
+
+
+_bicg_jit = _make_bicg_jit()
+
+
+def _csr_matvec_into(A):
+    """mv(x, y): y = A @ x, by scipy's own kernel into a reused buffer (A @ x zero-fills a fresh array and calls the
+    same csr_matvec; a fresh 16 MB array costs more in page faults than the product)."""
+    try:
+        from scipy.sparse import _sparsetools
+        fn = _sparsetools.csr_matvec
+    except (ImportError, AttributeError):
+        fn = None
+    M, N = A.shape
+
+    def mv(x, y):
+        if fn is None:
+            y[:] = A @ x
+            return
+        y.fill(0.0)
+        fn(M, N, A.indptr, A.indices, A.data, x, y)
+    return mv
+
+
+JIT_MIN_N = 50_000          # smaller systems: scipy (thread start-up outweighs the gain; same answer either way)
+
+
+def bicgstab_jacobi(A, b, x0, inv_diag, rtol, maxiter):
+    """scipy.sparse.linalg.bicgstab(A, b, x0=x0, rtol=rtol, maxiter=maxiter, M=diag(inv_diag)) for a CSR matrix and
+    float64 vectors: the same iterates, statement by statement (scipy 1.12+'s pure-Python loop), with the vector
+    updates fused and everything written into preallocated buffers (memory-bound; fresh 16 MB temporaries cost
+    more in page faults than in arithmetic). Falls back to scipy without numba."""
+    b = np.asarray(b, dtype=np.float64).ravel()
+    if (_bicg_jit is None or A.shape[0] < JIT_MIN_N or not sparse.isspmatrix_csr(A)
+            and not isinstance(A, sparse.csr_array) or A.dtype != np.float64):
+        M = splinalg.LinearOperator(A.shape, matvec=lambda x: inv_diag * x)
+        return splinalg.bicgstab(A, b, x0=x0, rtol=rtol, maxiter=maxiter, M=M)
+    p_update, scale, axpy_neg, x_update, axpy = _bicg_jit
+    mv = _csr_matvec_into(A)
+    inv_diag = np.ascontiguousarray(inv_diag, dtype=np.float64)
+    x = np.array(x0, dtype=np.float64)
+    bnrm2 = np.linalg.norm(b)
+    atol = max(0.0, float(rtol) * float(bnrm2))
+    if bnrm2 == 0:
+        return b, 0
+    rhotol = np.finfo(x.dtype.char).eps ** 2
+    omegatol = rhotol
+    rho_prev, omega, alpha, p, v = None, None, None, None, None
+    r = b - A @ x if x.any() else b.copy()
+    rtilde = r.copy()
+    phat, shat, v, t = np.empty_like(r), np.empty_like(r), np.empty_like(r), np.empty_like(r)
+    for iteration in range(maxiter):
+        if np.linalg.norm(r) < atol:
+            return x, 0
+        rho = np.dot(rtilde, r)
+        if np.abs(rho) < rhotol:
+            return x, -10
+        if iteration > 0:
+            if np.abs(omega) < omegatol:
+                return x, -11
+            beta = (rho / rho_prev) * (alpha / omega)
+            p_update(p, v, r, omega, beta)
+        else:
+            p = r.copy()
+        scale(phat, inv_diag, p)
+        mv(phat, v)
+        rv = np.dot(rtilde, v)
+        if rv == 0:
+            return x, -11
+        alpha = rho / rv
+        axpy_neg(r, alpha, v)
+        s = r                                     # scipy copies r into s here and reads both unchanged until r -= omega*t
+        if np.linalg.norm(s) < atol:
+            axpy(x, alpha, phat)
+            return x, 0
+        scale(shat, inv_diag, s)
+        mv(shat, t)
+        omega = np.dot(t, s) / np.dot(t, t)
+        x_update(x, alpha, phat, omega, shat)
+        axpy_neg(r, omega, t)
+        rho_prev = rho
+    return x, maxiter
+
+
+def drop_smooth_cache(g) -> None:
+    """Free the matrices smooth_km keeps on g."""
+    g.__dict__.pop("_screened", None)
+
+
+def gradient(g, f):
+    """Tangent-plane gradient (f per km): (2/k) Σ_j (f_j − f_i)/d_ij · t_ij."""
+    w = (f[g.dst] - f[g.src]) / g.edge_km
+    s = np.stack([np.bincount(g.src, weights=w * g.edge_tangents[:, c], minlength=g.n) for c in range(3)], axis=1)
+    return s * (2.0 / g.counts)[:, None]
+
+
+def _csr(g, weights):
+    return sparse.csr_matrix((weights, (g.src, g.dst)), shape=(g.n, g.n))
+
+
+def nearest_source(g, sources, weights=None):
+    sources = np.asarray(sources, dtype=np.int64)
+    if len(sources) == 0:
+        return np.full(g.n, np.inf), np.full(g.n, -9999, dtype=np.int64)
+    w = g.edge_km if weights is None else weights
+    dist, _, src = csgraph.dijkstra(_csr(g, w), directed=True, indices=sources,
+                                    min_only=True, return_predecessors=True)
+    return dist, src.astype(np.int64)
+
+
+def distance_to(g, mask):
+    return nearest_source(g, np.flatnonzero(mask))[0]
+
+
+def priority_flood(g, z, sink_mask, eps: float = 0.01):
+    """Barnes (2014) priority-flood + ε: every non-sink cell gets a strictly descending path to a sink."""
+    sink_mask = np.asarray(sink_mask, dtype=bool)
+    if not sink_mask.any():
+        raise ValueError("priority_flood: no sink cells")
+    has_open = np.bincount(g.src, weights=(~sink_mask)[g.dst].astype(np.float64), minlength=g.n) > 0
+    seeds = np.flatnonzero(sink_mask & has_open)
+    z = np.asarray(z, dtype=np.float64)
+    if _flood_jit is not None:
+        return _flood_jit(z.copy(), sink_mask.copy(), np.asarray(g.nbr_ptr, np.int64), np.asarray(g.nbr_idx, np.int64),
+                          seeds.astype(np.int64), float(eps))
+    return _flood_py(z, sink_mask, g.nbr_ptr, g.nbr_idx, seeds, eps)
+
+
+def _flood_py(z, sink_mask, nbr_ptr, nbr_idx, seeds, eps):
+    zf = np.asarray(z, dtype=np.float64).tolist()
+    done = np.asarray(sink_mask).tolist()
+    ptr = np.asarray(nbr_ptr).tolist()
+    idx = np.asarray(nbr_idx).tolist()
+    heap = [(zf[i], i) for i in np.asarray(seeds).tolist()]
+    heapq.heapify(heap)
+    while heap:
+        zc, c = heapq.heappop(heap)
+        for k in range(ptr[c], ptr[c + 1]):
+            n = idx[k]
+            if not done[n]:
+                done[n] = True
+                if zf[n] < zc + eps:
+                    zf[n] = zc + eps
+                heapq.heappush(heap, (zf[n], n))
+    return np.array(zf)
+
+
+def _make_flood_jit():
+    """_flood_py compiled with numba when it is installed (optional: same heap order, same float steps, same result;
+    WORLDGEN_NO_JIT=1 turns it off)."""
+    import os
+    if os.environ.get("WORLDGEN_NO_JIT"):
+        return None
+    try:
+        import numba
+    except ImportError:
+        return None
+
+    @numba.njit(cache=True)
+    def flood(zf, done, ptr, idx, seeds, eps):
+        heap = [(zf[i], i) for i in seeds]
+        heapq.heapify(heap)
+        while len(heap):
+            zc, c = heapq.heappop(heap)
+            for k in range(ptr[c], ptr[c + 1]):
+                n = idx[k]
+                if not done[n]:
+                    done[n] = True
+                    if zf[n] < zc + eps:
+                        zf[n] = zc + eps
+                    heapq.heappush(heap, (zf[n], n))
+        return zf
+    return flood
+
+
+_flood_jit = _make_flood_jit()
+
+
+def _make_steep_jit():
+    """steepest_receivers' per-row maximum, compiled (numba optional, as the flood): the first edge in row order
+    with the largest slope, slopes computed as numpy does. ok=False (empty row, NaN): use the numpy path."""
+    import os
+    if os.environ.get("WORLDGEN_NO_JIT"):
+        return None
+    try:
+        import numba
+    except ImportError:
+        return None
+
+    @numba.njit(cache=True)
+    def steep(z, ptr, dst, edge):
+        n = len(ptr) - 1
+        first = np.empty(n, np.int64)
+        best = np.empty(n)
+        for i in range(n):
+            a, b = ptr[i], ptr[i + 1]
+            if a == b:
+                return first, best, False
+            bk = a
+            bs = (z[i] - z[dst[a]]) / edge[a]
+            if bs != bs:
+                return first, best, False
+            for k in range(a + 1, b):
+                sk = (z[i] - z[dst[k]]) / edge[k]
+                if sk != sk:
+                    return first, best, False
+                if sk > bs:
+                    bs, bk = sk, k
+            first[i], best[i] = bk, bs
+        return first, best, True
+    return steep
+
+
+_steep_jit = _make_steep_jit()
+
+
+def steepest_receivers(g, z):
+    """Per cell: the neighbour of steepest descent (first in neighbour order on ties), the slope and the edge length;
+    no way down → itself, 0, inf."""
+    ptr = np.asarray(g.nbr_ptr, np.int64)
+    if _steep_jit is not None and g.n:
+        first, s, ok = _steep_jit(np.asarray(z), ptr, np.asarray(g.nbr_idx, np.int64), g.edge_km)
+        if ok:
+            down = s > 0
+            recv = np.where(down, g.dst[first], np.arange(g.n))
+            return recv, np.where(down, s, 0.0), np.where(down, g.edge_km[first], np.inf)
+    slope = (z[g.src] - z[g.dst]) / g.edge_km
+    if g.n == 0 or not (np.diff(ptr) > 0).all() or np.isnan(slope).any():
+        order = np.lexsort((-slope, g.src))                  # general case (empty rows, NaN)
+        first = order[ptr[:-1]]
+    else:                                                    # same edge as the stable lexsort, without sorting
+        smax = np.maximum.reduceat(slope, ptr[:-1])
+        cand = np.flatnonzero(slope == smax[g.src])
+        rows = g.src[cand]
+        first = cand[np.concatenate([[True], rows[1:] != rows[:-1]])]
+    s = slope[first]
+    down = s > 0
+    ar = np.arange(g.n)
+    recv = np.where(down, g.dst[first], ar)
+    return recv, np.where(down, s, 0.0), np.where(down, g.edge_km[first], np.inf)
+
+
+def _gather(ptr, arr, sel):
+    counts = ptr[sel + 1] - ptr[sel]
+    tot = int(counts.sum())
+    if tot == 0:
+        return arr[:0]
+    starts = np.repeat(ptr[sel] - np.concatenate([[0], np.cumsum(counts)[:-1]]), counts)
+    return arr[starts + np.arange(tot)]
+
+
+def receiver_levels(recv):
+    recv = np.asarray(recv, dtype=np.int64)
+    n = len(recv)
+    ar = np.arange(n)
+    root = recv == ar
+    donors = ar[~root]
+    donors = donors[np.argsort(recv[donors], kind="stable")]
+    dptr = np.concatenate([[0], np.cumsum(np.bincount(recv[donors], minlength=n))])
+    levels, frontier, seen = [], ar[root], 0
+    while len(frontier):
+        levels.append(frontier)
+        seen += len(frontier)
+        frontier = _gather(dptr, donors, frontier)
+    if seen != n:
+        raise ValueError("receiver_levels: cycle in receivers")
+    return levels
+
+
+def accumulate(recv, levels, w):
+    """Sum w down the receiver tree. Per level only the receivers are touched (a full bincount per level costs
+    levels × cells); the sums are added in donor order from 0, as bincount does: the same floats."""
+    acc = np.asarray(w, dtype=np.float64).copy()
+    buf = np.zeros(len(acc))
+    if len(levels) > 1:
+        acc += 0.0                                    # as the first full-length add did: −0 becomes +0
+    for lv in reversed(levels[1:]):
+        r = recv[lv]
+        buf[r] = 0.0
+        np.add.at(buf, r, acc[lv])
+        acc[r] = acc[r] + buf[r]
+    return acc
+
+
+def _make_components_jit():
+    """Connected-component labels as scipy's connected_components numbers them — each component (every node not in
+    the mask is one by itself) by the order of its lowest node — by union-find over the edges, with no sparse matrix
+    (numba optional; WORLDGEN_NO_JIT=1 turns it off)."""
+    import os
+    if os.environ.get("WORLDGEN_NO_JIT"):
+        return None
+    try:
+        import numba
+    except ImportError:
+        return None
+
+    @numba.njit(cache=True, nogil=True)
+    def labels(n, src, dst, mask):
+        parent = np.arange(n)
+        for e in range(src.shape[0]):
+            a, b = src[e], dst[e]
+            if mask[a] and mask[b]:
+                while parent[a] != a:
+                    parent[a] = parent[parent[a]]
+                    a = parent[a]
+                while parent[b] != b:
+                    parent[b] = parent[parent[b]]
+                    b = parent[b]
+                if a != b:
+                    if a < b:
+                        parent[b] = a
+                    else:
+                        parent[a] = b
+        root_lab = np.full(n, -1, np.int32)
+        out = np.empty(n, np.int32)
+        count = 0
+        for v in range(n):
+            r = v
+            while parent[r] != r:
+                r = parent[r]
+            if root_lab[r] < 0:
+                root_lab[r] = count
+                count += 1
+            out[v] = root_lab[r] if mask[v] else -1
+        return out
+    return labels
+
+
+_components_jit = _make_components_jit()
+
+
+def components(g, mask):
+    mask = np.asarray(mask, dtype=bool)
+    if _components_jit is not None:
+        return _components_jit(g.n, g.src, g.dst, mask)
+    e = mask[g.src] & mask[g.dst]
+    m = sparse.csr_matrix((np.ones(int(e.sum())), (g.src[e], g.dst[e])), shape=(g.n, g.n))
+    _, lab = csgraph.connected_components(m, directed=False)
+    return np.where(mask, lab, -1)
+
+
+OCEAN_MIN_KM2 = 5.0e6
+
+
+def ocean_mask(g, z, min_area_km2: float = OCEAN_MIN_KM2):
+    """The connected world ocean plus any separate basin ≥ min_area_km2; smaller interior lows count as land."""
+    wet = np.asarray(z) <= 0
+    if not wet.any():
+        return wet
+    lab = components(g, wet)
+    area = np.bincount(lab[wet], weights=g.area_km2[wet])
+    keep = area >= min_area_km2
+    keep[np.argmax(area)] = True
+    return wet & keep[np.maximum(lab, 0)]
diff --git a/mapgen/grid.py b/mapgen/grid.py
new file mode 100644
index 0000000..8138b5f
--- /dev/null
+++ b/mapgen/grid.py
@@ -0,0 +1,121 @@
+"""H3 hexagonal grid as a CSR cell graph. Stage `grid`."""
+from __future__ import annotations
+
+from dataclasses import dataclass
+from functools import cached_property
+
+import h3.api.basic_int as h3
+import numpy as np
+
+from .sphere import latlon_to_xyz, tangent_dir
+
+
+EDGE_BLOCK = 1 << 20
+
+
+def edge_blocks(n: int, block: int | None = None):
+    """Slices over n edges, EDGE_BLOCK at a time: per-edge (row-wise) formulas give the same values block by block,
+    with (block, 3) float64 temporaries instead of (edges, 3) ones (290 MB each at r5, 2 GB at r6)."""
+    block = block or EDGE_BLOCK
+    return [slice(a, min(a + block, n)) for a in range(0, n, block)]
+
+
+def rowdot_at(A, idx, t) -> np.ndarray:
+    """np.sum(A[idx] * t, axis=1) — per edge, A's row at an end dotted with the edge's vector — in edge blocks."""
+    out = np.empty(len(idx))
+    for s in edge_blocks(len(idx)):
+        out[s] = np.sum(A[idx[s]] * t[s], axis=1)
+    return out
+
+
+def cell_parents(cells, res: int) -> np.ndarray:
+    """h3.cell_to_parent for an array of cells (all of resolution ≥ res), by the index bits: the resolution field
+    set to res, the digits below it set to 7 (unused). Same ids as h3, without a Python call per cell."""
+    cells = np.asarray(cells, dtype=np.uint64)
+    if len(cells) and int(((cells >> np.uint64(52)) & np.uint64(15)).min()) < res:
+        raise ValueError(f"cell_parents: a cell is coarser than resolution {res}")
+    unused = 0
+    for r in range(res + 1, 16):
+        unused |= 7 << ((15 - r) * 3)
+    out = (cells & np.uint64(~(15 << 52) & (2**64 - 1))) | np.uint64(res << 52)
+    return out | np.uint64(unused)
+
+
+@dataclass
+class Grid:
+    res: int
+    radius_km: float
+    ids: np.ndarray
+    lat: np.ndarray
+    lon: np.ndarray
+    xyz: np.ndarray
+    area_km2: np.ndarray
+    nbr_ptr: np.ndarray
+    nbr_idx: np.ndarray
+
+    @property
+    def n(self) -> int:
+        return len(self.ids)
+
+    @cached_property
+    def counts(self):
+        return np.diff(self.nbr_ptr)
+
+    @cached_property
+    def src(self):
+        return np.repeat(np.arange(self.n), self.counts)
+
+    @property
+    def dst(self):
+        return self.nbr_idx
+
+    @cached_property
+    def edge_km(self):
+        out = np.empty(len(self.dst))
+        for s in edge_blocks(len(self.dst)):
+            d = np.sum(self.xyz[self.src[s]] * self.xyz[self.dst[s]], axis=1)
+            out[s] = self.radius_km * np.arccos(np.clip(d, -1.0, 1.0))
+        return out
+
+    @cached_property
+    def edge_tangents(self):
+        out = np.empty((len(self.dst), 3))
+        for s in edge_blocks(len(self.dst)):
+            out[s] = tangent_dir(self.xyz[self.src[s]], self.xyz[self.dst[s]])
+        return out
+
+    @cached_property
+    def spacing_km(self) -> float:
+        return float(self.edge_km.mean())
+
+    def cell_index(self, lat: float, lon: float) -> int:
+        c = np.uint64(h3.latlng_to_cell(float(lat), float(lon), self.res))
+        return int(np.searchsorted(self.ids, c))
+
+    def to_arrays(self) -> dict:
+        return {"g_ids": self.ids, "g_lat": self.lat, "g_lon": self.lon, "g_xyz": self.xyz,
+                "g_area_km2": self.area_km2, "g_nbr_ptr": self.nbr_ptr, "g_nbr_idx": self.nbr_idx}
+
+    @classmethod
+    def from_arrays(cls, d: dict, res: int, radius_km: float) -> "Grid":
+        return cls(res, radius_km, d["g_ids"], d["g_lat"], d["g_lon"], d["g_xyz"], d["g_area_km2"],
+                   d["g_nbr_ptr"], d["g_nbr_idx"])
+
+
+def build_grid(res: int, radius_km: float) -> Grid:
+    cells: list[int] = []
+    for r0 in h3.get_res0_cells():
+        cells.extend(h3.cell_to_children(r0, res))
+    ids = np.array(sorted(cells), dtype=np.uint64)
+    ll = np.array([h3.cell_to_latlng(int(c)) for c in ids], dtype=np.float64)
+    area = np.array([h3.cell_area(int(c), unit="rads^2") for c in ids]) * radius_km**2
+    rings = [h3.grid_ring(int(c), 1) for c in ids]
+    counts = np.array([len(r) for r in rings], dtype=np.int64)
+    ptr = np.concatenate([[0], np.cumsum(counts)]).astype(np.int64)
+    flat = np.array([c for r in rings for c in r], dtype=np.uint64)
+    idx = np.searchsorted(ids, flat).astype(np.int64)
+    return Grid(res, radius_km, ids, ll[:, 0], ll[:, 1], latlon_to_xyz(ll[:, 0], ll[:, 1]), area, ptr, idx)
+
+
+def run(ctx) -> dict:
+    return build_grid(ctx.res, float(ctx.cfg["planet"]["radius_km"])).to_arrays()
diff --git a/mapgen/hydrology.py b/mapgen/hydrology.py
new file mode 100644
index 0000000..d7adde1
--- /dev/null
+++ b/mapgen/hydrology.py
@@ -0,0 +1,231 @@
+"""Stage `hydrology`: depression filling, lakes vs endorheic basins, discharge, rivers, Strahler order."""
+from __future__ import annotations
+
+import numpy as np
+from scipy import sparse
+from scipy.sparse import csgraph
+
+from .config import params
+from .crust import RIFT
+from .graph import accumulate, components, distance_to, priority_flood, receiver_levels, steepest_receivers
+from .noise import fbm
+from .pipeline import StageError
+
+DEFAULTS = {"fill_eps_m": 0.01, "min_depth_m": 1.0, "river_min_km3_yr": 2.0,
+            "salt_flat_max_p_mm": 300.0, "salt_flat_fraction": 0.2, "dry_lake_fraction": 0.05,
+            # lake beds: deepest point = k × area^exp m (Earth-like: ≈95 m at
+            # 1,000 km², ≈230 m at 20,000 km², ≈610 m at 500,000 km²) × 0.5–2 (noise at the lake), × rift_x in
+            # rifts, × arid_x for dry terminal lakes; never shallower than lake_min_m (refinement keeps ≥ 15 m)
+            "lake_depth_k": 12.0, "lake_depth_exp": 0.3, "lake_min_m": 25.0, "lake_max_m": 1800.0,
+            "lake_rift_x": 2.5, "lake_arid_x": 0.4, "lake_arid_p_mm": 400.0, "lake_shore": 0.25}
+
+
+def strahler(recv, levels, river):
+    n = len(recv)
+    order = np.zeros(n, np.int8)
+    mx = np.zeros(n, np.int8)
+    cnt = np.zeros(n, np.int16)
+    for lv in reversed(levels):
+        cells = lv[river[lv]]
+        if len(cells) == 0:
+            continue
+        order[cells] = np.where(cnt[cells] >= 2, mx[cells] + 1, np.maximum(mx[cells], 1))
+        nonroot = cells[recv[cells] != cells]
+        r, o = recv[nonroot], order[nonroot]
+        np.maximum.at(mx, r, o)
+        np.add.at(cnt, r, (o == mx[r]).astype(np.int16))
+    return order
+
+
+def _groups(lab):
+    idx = np.argsort(lab, kind="stable")
+    ls = lab[idx]
+    start = np.searchsorted(ls, 0)
+    idx, ls = idx[start:], ls[start:]
+    cuts = np.flatnonzero(np.diff(ls)) + 1
+    return np.split(idx, cuts) if len(idx) else []
+
+
+def _leaves(recv, lab, x, limit=100000):
+    """Does the flow from depression cell x leave its depression for good (not back in over shallow ground)?"""
+    own, y = lab[x], recv[x]
+    for _ in range(limit):
+        if lab[y] == own:
+            return False
+        if lab[y] >= 0 or recv[y] == y:
+            return True
+        y = recv[y]
+    return True
+
+
+def _leaves_all(recv, lab, levels, limit=100000):
+    """_leaves for every cell at once: walking down from recv[x], the first cell that is in a depression or a root
+    (stop) and how many steps away it is; x leaves unless that cell is in x's own depression (or the walk is longer
+    than limit steps)."""
+    n = len(recv)
+    stop, d = np.arange(n), np.zeros(n, np.int64)
+    for lv in levels[1:]:                                    # roots first: a cell's receiver is done before it
+        free = lv[lab[lv] < 0]
+        stop[free] = stop[recv[free]]
+        d[free] = d[recv[free]] + 1
+    y = recv
+    return (d[y] >= limit) | (lab[stop[y]] != lab)
+
+
+def _route_inside(g, dep, lab, recv, targets):
+    """Within each depression, point every cell along the shortest intra-depression path to its target cell."""
+    e = dep[g.src] & dep[g.dst] & (lab[g.src] == lab[g.dst])
+    m = sparse.csr_matrix((g.edge_km[e], (g.src[e], g.dst[e])), shape=(g.n, g.n))
+    _, pred, _ = csgraph.dijkstra(m, directed=False, indices=targets, min_only=True, return_predecessors=True)
+    out = recv.copy()
+    inside = dep & (pred >= 0)
+    out[inside] = pred[inside]
+    return out
+
+
+def _sweep(recv, levels, water, outlets, cap):
+    """Accumulate flow upstream→downstream; at each spilling-lake outlet remove up to `cap` (lake evaporation)."""
+    acc = np.asarray(water, dtype=np.float64).copy()
+    loss = np.zeros(len(acc))
+    is_out = np.zeros(len(acc), bool)
+    is_out[outlets] = True
+    cap_cell = np.zeros(len(acc))
+    cap_cell[outlets] = cap
+    buf = np.zeros(len(acc))                    # per level only the receivers change (graph.accumulate): same floats
+    if len(levels) > 1:
+        acc += 0.0
+    for lv in reversed(levels[1:]):
+        push = acc[lv].copy()
+        o = is_out[lv]
+        if o.any():
+            cells = lv[o]
+            lost = np.minimum(acc[cells], cap_cell[cells])
+            loss[cells] = lost
+            push[o] = acc[cells] - lost
+        r = recv[lv]
+        buf[r] = 0.0
+        np.add.at(buf, r, push)
+        acc[r] = acc[r] + buf[r]
+    return acc, loss
+
+
+def lake_levels(z, zf, lab, lake, lake_id, endo):
+    """Each lake cell's water level (NaN elsewhere): a spilling lake stands at its spill height; a terminal lake at the
+    lowest ground of its basin it does not cover (the next cell to flood). Carving the beds leaves both unchanged."""
+    level = np.full(len(z), np.nan)
+    m = int(lab.max()) + 1 if len(lab) and lab.max() >= 0 else 0
+    dry = np.full(m, np.inf)                                    # per depression: its lowest uncovered ground
+    sel = (lab >= 0) & ~lake
+    np.minimum.at(dry, lab[sel], z[sel])
+    for c in _groups(np.where(lake, lake_id, -1)):
+        k = lab[c[0]]
+        level[c] = dry[k] if endo[c[0]] and k >= 0 and np.isfinite(dry[k]) else zf[c].min()
+    return level
+
+
+def lake_beds(g, z, level, lake, lake_id, endo, rift, p_ann, seed, P, only=None):
+    """Lake beds carved below their level (heights unchanged elsewhere): each lake's deepest point from its area
+    (P lake_*), noise keyed by where the lake lies (not its id), the depth rising from lake_shore × that at the shore
+    to all of it at the cell farthest from shore; never above the ground (min), never shallower than lake_min_m.
+    only: carve just the lakes touching these cells (eras: where events changed the ground)."""
+    z = np.asarray(z, dtype=np.float64).copy()
+    if not lake.any():
+        return z
+    shore = distance_to(g, ~lake) if (~lake).any() else np.full(g.n, g.spacing_km)
+    for c in _groups(np.where(lake, lake_id, -1)):
+        if only is not None and not only[c].any():
+            continue
+        area = g.area_km2[c].sum()
+        ctr = g.xyz[c].mean(0)
+        ctr = ctr / max(np.linalg.norm(ctr), 1e-12)
+        d = P["lake_depth_k"] * area ** P["lake_depth_exp"] * 2.0 ** float(fbm(ctr[None], seed + 4421, 3, 6.0)[0] * 1.4)
+        if rift[c].mean() > 0.3:
+            d *= P["lake_rift_x"]
+        if endo[c[0]] and p_ann[c].mean() < P["lake_arid_p_mm"]:
+            d *= P["lake_arid_x"]
+        d = float(np.clip(d, P["lake_min_m"], P["lake_max_m"]))
+        t = shore[c] / max(shore[c].max(), 1e-9)
+        prof = np.maximum(d * (P["lake_shore"] + (1 - P["lake_shore"]) * t ** 0.6), P["lake_min_m"])
+        z[c] = np.minimum(z[c], level[c] - prof)
+    return z
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "hydrology", DEFAULTS)
+    z, p_ann, pet = ctx.need("elevation_eroded_m", "P_ann", "PET")
+    z = z.astype(np.float64)
+    ar = np.arange(g.n)
+    ocean = np.asarray(ctx.data["ocean"]) if "ocean" in ctx.data else z <= 0
+    aet = p_ann / np.sqrt(1.0 + (p_ann / np.maximum(pet, 1e-6)) ** 2)      # Pike (1964)
+    runoff = np.where(ocean, 0.0, np.maximum(p_ann - aet, 0.0))
+    water = runoff * g.area_km2 * 1e-6                                       # km³/yr per cell
+    zf = priority_flood(g, z, ocean, P["fill_eps_m"])
+    depth = zf - z
+    dep = ~ocean & (depth > P["min_depth_m"])
+    lab = components(g, dep)
+    recv1, _, _ = steepest_receivers(g, zf)
+    recv1 = np.where(ocean, ar, recv1)
+    levels1 = receiver_levels(recv1)
+    q1 = accumulate(recv1, levels1, water)
+    evap_net = np.maximum(pet - p_ann, 0.0) * g.area_km2 * 1e-6              # full-lake evaporation, km³/yr
+
+    groups = _groups(lab)
+    leaves = _leaves_all(recv1, lab, levels1) if len(groups) else None
+    outlets = np.zeros(len(groups), np.int64)
+    terminals = np.zeros(len(groups), np.int64)
+    cap = np.zeros(len(groups))
+    for k, c in enumerate(groups):
+        ext = c[lab[recv1[c]] != lab[c[0]]]
+        ext = ext[leaves[ext]]
+        outlets[k] = ext[np.argmax(q1[ext])] if len(ext) else c[np.argmax(q1[c])]
+        terminals[k] = c[np.argmin(z[c])]
+        cap[k] = evap_net[c].sum()
+
+    # 1) every depression drains to its spill outlet; one upstream-first sweep decides spill vs endorheic
+    recv_a = _route_inside(g, dep, lab, recv1, outlets) if len(groups) else recv1
+    acc_a, _ = _sweep(recv_a, receiver_levels(recv_a), water, outlets, cap)
+    inflow = acc_a[outlets]
+    spill = inflow >= cap
+    # 2) endorheic depressions drain to their single lowest cell instead
+    targets = np.where(spill, outlets, terminals)
+    recv = _route_inside(g, dep, lab, recv1, targets) if len(groups) else recv1
+    recv[terminals[~spill]] = terminals[~spill]
+    try:
+        levels = receiver_levels(recv)
+    except ValueError as e:
+        raise StageError(f"hydrology: {e}") from e
+    q, loss = _sweep(recv, levels, water, outlets[spill], cap[spill])
+
+    lake = np.zeros(g.n, bool)
+    endo = np.zeros(g.n, bool)
+    salt = np.zeros(g.n, bool)
+    lake_id = np.full(g.n, -1, np.int32)
+    for k, c in enumerate(groups):
+        if spill[k]:
+            lake[c] = True
+            lake_id[c] = k
+            continue
+        endo[c] = True
+        cs = c[np.argsort(z[c])]
+        nl = int(np.searchsorted(np.cumsum(evap_net[cs]), inflow[k]))
+        if inflow[k] > 0:
+            nl = max(nl, 1)                         # the terminal always holds some water
+        lake[cs[:nl]] = True
+        lake_id[cs[:nl]] = k
+        if nl == 0 or (nl < P["dry_lake_fraction"] * len(cs) and p_ann[c].mean() < P["salt_flat_max_p_mm"]):
+            a = np.cumsum(g.area_km2[cs])
+            ns = max(nl + 1, int(np.searchsorted(a, P["salt_flat_fraction"] * a[-1])))
+            salt[cs[nl:ns]] = True
+    river = ~ocean & ~lake & (q >= P["river_min_km3_yr"])
+    level = lake_levels(z, zf, lab, lake, lake_id, endo)
+    rift = np.asarray(ctx.data["age_class"]) == RIFT if "age_class" in ctx.data else np.zeros(g.n, bool)
+    zb = lake_beds(g, z, level, lake, lake_id, endo, rift, p_ann, ctx.seed, P, ctx.data.get("lake_carve"))
+    cut = np.asarray(ctx.data.get("lake_cut_m", 0.0), dtype=np.float64) + (z - zb)   # sea masks see the uncut ground
+    z = zb
+    depth = zf - z
+    dtype = np.asarray(ctx.data["elevation_eroded_m"]).dtype
+    return {"elevation_eroded_m": z.astype(dtype), "lake_level_m": level.astype(np.float32),
+            "lake_cut_m": np.broadcast_to(cut, z.shape).astype(np.float32), "z_filled_m": zf, "recv": recv, "discharge_km3_yr": q, "runoff_mm": runoff, "aet_mm": aet,
+            "lake": lake, "lake_id": lake_id, "endorheic": endo, "salt_flat": salt, "river": river,
+            "strahler": strahler(recv, levels, river), "depression_depth_m": depth, "lake_loss_km3_yr": loss}
diff --git a/mapgen/ice.py b/mapgen/ice.py
new file mode 100644
index 0000000..2aef33f
--- /dev/null
+++ b/mapgen/ice.py
@@ -0,0 +1,37 @@
+"""Stage `ice`: ice sheets, mountain glaciers, seasonal/perennial sea ice."""
+from __future__ import annotations
+
+import numpy as np
+
+from .config import params
+from .graph import components, distance_to
+
+ICE_NONE, ICE_SHEET, ICE_GLACIER, ICE_SEA_SEASONAL, ICE_SEA_PERENNIAL = range(5)
+ICE_NAMES = ["none", "ice sheet", "glacier", "seasonal sea ice", "perennial sea ice"]
+DEFAULTS = {"melt_summer_c": 0.0, "sheet_min_km2": 250000.0, "sheet_min_p_mm": 100.0, "sea_ice_t_c": -1.8, "sheet_max_m": 3000.0,
+            "sheet_edge_m": 200.0, "sheet_growth_m_per_km": 3.0}
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "ice", DEFAULTS)
+    z, t_mean, t_jun, t_dec, p_ann = ctx.need("elevation_eroded_m", "T_mean", "T_jun", "T_dec", "P_ann")
+    z = z.astype(np.float64)
+    land = ~np.asarray(ctx.data["ocean"]) if "ocean" in ctx.data else z > 0
+    summer = np.where(g.lat >= 0, t_jun, t_dec)
+    winter = np.where(g.lat >= 0, t_dec, t_jun)
+    cold = land & (summer < P["melt_summer_c"])                  # snow survives the summer → permanent ice
+    lab = components(g, cold)
+    area = np.bincount(lab[cold], weights=g.area_km2[cold]) if cold.any() else np.zeros(1)
+    big = cold & (area[np.maximum(lab, 0)] >= P["sheet_min_km2"])
+    sheet = big & (p_ann > P["sheet_min_p_mm"])
+    thick = np.zeros(g.n)
+    if sheet.any():
+        d_edge = distance_to(g, ~sheet)
+        thick = np.where(sheet, np.minimum(P["sheet_max_m"], P["sheet_edge_m"] + P["sheet_growth_m_per_km"] * d_edge), 0.0)
+    ice = np.zeros(g.n, np.int8)
+    ice[cold & ~sheet] = ICE_GLACIER
+    ice[sheet] = ICE_SHEET
+    ice[~land & (winter < P["sea_ice_t_c"])] = ICE_SEA_SEASONAL
+    ice[~land & (summer < P["sea_ice_t_c"])] = ICE_SEA_PERENNIAL
+    return {"ice": ice, "ice_thickness_m": thick, "z_surface_m": z + thick}
diff --git a/mapgen/minerals.py b/mapgen/minerals.py
new file mode 100644
index 0000000..ab173c9
--- /dev/null
+++ b/mapgen/minerals.py
@@ -0,0 +1,172 @@
+"""Mineral deposits: where each ore, fuel and industrial mineral can be mined, by
+geological rules, placed at an Earth-like share of the land. Everything a society needs to reach steam, rifled guns
+and aluminium airships is somewhere; the good stuff (high-grade iron ore and the alloy metals) lies mostly deep in
+mountains. Each deposit is a bit of `deposits` (uint32); `deposit_main` names the rarest one in a cell (a map layer).
+
+A deposit's cells: the land cells (not sea, not under lakes) where its rule allows it, scored by the rule × a district
+noise of its own (deposits cluster into mining districts). Half the `share` is taken province by province (≈ 1500 km
+cells, each its share of its own land, best-scored first): every region has its local bog iron, coal pits and small
+mines where its rocks allow. The rest goes to the best-scored land world-wide (the great districts)."""
+from __future__ import annotations
+
+import numpy as np
+
+from scipy.spatial import cKDTree
+
+from .noise import fbm, name_seed
+
+# name, label, share of land (fraction). Shares are rough Earth-like extents of workable districts, not grades.
+DEPOSITS = [
+    ("bog_iron", "bog / laterite iron (low grade)", 0.030),
+    ("ironstone", "ironstone (sedimentary iron)", 0.025),
+    ("iron_high", "high-grade iron (magnetite, banded)", 0.006),
+    ("coal", "coal (incl. coking)", 0.050),
+    ("copper", "copper", 0.012),
+    ("molybdenum", "molybdenum", 0.002),
+    ("tin", "tin", 0.003),
+    ("tungsten", "tungsten", 0.002),
+    ("lead_zinc", "lead & zinc", 0.010),
+    ("silver", "silver", 0.005),
+    ("gold", "gold", 0.008),
+    ("manganese", "manganese", 0.004),
+    ("chromium", "chromium", 0.0015),
+    ("nickel", "nickel", 0.003),
+    ("cobalt", "cobalt", 0.001),
+    ("vanadium", "vanadium & titanium", 0.002),
+    ("platinum", "platinum metals", 0.0005),
+    ("mercury", "mercury (cinnabar)", 0.0012),
+    ("sulfur", "sulfur", 0.004),
+    ("saltpetre", "saltpetre (nitrates)", 0.004),
+    ("bauxite", "bauxite (aluminium)", 0.010),
+    ("fluorite", "fluorite / cryolite (flux)", 0.0025),
+    ("rock_salt", "rock salt & evaporites", 0.018),
+    ("potash", "potash", 0.004),
+    ("phosphate", "phosphate", 0.005),
+    ("oil_gas", "oil & gas", 0.025),
+    ("helium", "helium (with gas)", 0.0025),
+    ("diamond", "diamonds (kimberlite)", 0.0004),
+    ("magnesite", "magnesite / dolomite (magnesium)", 0.003),
+    ("kaolin", "kaolin / fire clay", 0.010),
+    ("glass_sand", "glass sand (quartz)", 0.015),
+    ("graphite", "graphite", 0.002),
+]
+NAMES = [d[0] for d in DEPOSITS]
+BIT = {n: 1 << i for i, n in enumerate(NAMES)}
+LEGEND = ["—"] + [d[1] for d in DEPOSITS]
+# the good stuff: mostly in deep mountains
+DEEP = ("iron_high", "tungsten", "molybdenum", "chromium", "cobalt", "vanadium", "platinum", "manganese", "tin")
+
+
+PROVINCE_KM = 1500.0   # spacing of the provinces that each take their local share
+LOCAL = 0.5            # part of a deposit's share taken province by province
+LOCAL_MIN = 0.25       # ... only on ground at least this good relative to the world's marginal district
+
+
+def provinces(g, seed):
+    """Province index of each cell: the nearest of random centres ≈ PROVINCE_KM apart."""
+    n = max(1, int(round(4 * np.pi * g.radius_km ** 2 / PROVINCE_KM ** 2)))
+    c = np.random.default_rng(name_seed(seed, "deposit-provinces")).normal(size=(n, 3))
+    c /= np.linalg.norm(c, axis=1)[:, None]
+    return cKDTree(c).query(g.xyz / np.linalg.norm(g.xyz, axis=1)[:, None])[1]
+
+
+def _rank_select(score, share, land, prov):
+    """`share` × land cells where score > 0: LOCAL of it the best of each province (its share of its own land, on
+    ground ≥ LOCAL_MIN × the world's k-th best score), the rest the best world-wide."""
+    ok = np.flatnonzero((score > 0) & land)
+    k = min(len(ok), int(round(share * land.sum())))
+    out = np.zeros(len(score), bool)
+    if k == 0:
+        return out
+    cut = LOCAL_MIN * np.partition(score[ok], len(ok) - k)[len(ok) - k]     # the k-th best world-wide
+    ok_l = ok[score[ok] >= cut]
+    quota = np.floor(LOCAL * share * np.bincount(prov[land], minlength=prov.max() + 1) + 0.5).astype(np.int64)
+    o = ok_l[np.lexsort((-score[ok_l], prov[ok_l]))]           # by province, best first
+    p = prov[o]
+    first = np.r_[0, np.flatnonzero(p[1:] != p[:-1]) + 1]
+    rank = np.arange(len(o)) - np.repeat(first, np.diff(np.r_[first, len(o)]))
+    local = o[rank < quota[p]]
+    out[local[np.argsort(-score[local], kind="stable")[:k]]] = True
+    rest = ok[~out[ok]]
+    m = k - int(out.sum())
+    if m > 0:
+        out[rest[np.argpartition(-score[rest], m - 1)[:m]]] = True
+    return out
+
+
+def rules(f):
+    """Each deposit's suitability (≥ 0, 0 = impossible) from the geology/climate fields in f."""
+    from .crust import CRATON, LIP, POST_OROGEN, PRE_OROGEN, RIFT, VOLCANO
+    from .environment import (LF_ARC, LF_MASSIF, LF_MOUNTAINS, LI_ANDESITE, LI_BASALT, LI_GRANITE, LI_LIMESTONE,
+                              LI_METAMORPHIC, LI_SANDSTONE, LF_DUNES, LF_PLAIN)
+    z, lit, age, lf, rel = f["z"], f["lit"], f["age"], f["lf"], f["relief"]
+    t, p = f["t_mean"], f["p_ann"]
+    b = lambda m: np.asarray(m, dtype=np.float64)
+    deep = np.clip((z - 800.0) / 2000.0, 0, 1) * np.clip(rel / 1500.0, 0.3, 1.0)   # deep-mountain mining ground
+    deep = np.maximum(deep, 0.6 * b(np.isin(lf, [LF_MOUNTAINS, LF_MASSIF])) * np.clip(z / 2500.0, 0, 1))
+    oro = b(np.isin(age, [PRE_OROGEN, POST_OROGEN]))
+    craton = b(age == CRATON)
+    rift = b(age == RIFT)
+    arc = b((lit == LI_ANDESITE) | (lf == LF_ARC))
+    volc = b(age == VOLCANO) + arc
+    gran = b(lit == LI_GRANITE)
+    meta = b(lit == LI_METAMORPHIC)
+    lime = b(lit == LI_LIMESTONE)
+    sand = b(lit == LI_SANDSTONE)
+    basalt = b((lit == LI_BASALT) | (age == LIP))
+    sed = lime + sand
+    low = b(z < 800)
+    hot, wet = b(t > 20), b(p > 1200)
+    arid = b(p < 250) * b(t > 10)
+    suture = b(f["d_coll"] < 400) * (oro + meta)
+    return {
+        "bog_iron": b(p > 700) * low * (b(lf == LF_PLAIN) + 0.5) + hot * wet * low,
+        "ironstone": sed * b(z < 1200),
+        "iron_high": deep * (1.0 + craton + meta + oro),
+        "coal": b(f["coal"]),
+        "copper": arc * (1 + deep) + 0.4 * sand * rift + 0.3 * basalt,
+        "molybdenum": arc * deep,
+        "tin": (gran + meta) * oro * (0.3 + deep),
+        "tungsten": (gran + meta) * (oro + 0.3) * deep,
+        "lead_zinc": lime + 0.3 * sand * (1 + rift),
+        "silver": arc + 0.3 * lime + 0.3 * deep * oro,
+        "gold": (oro + craton * 0.7 + meta) * (0.3 + deep) + 0.3 * arc,
+        "manganese": (sed + craton) * (0.2 + deep),
+        "chromium": suture * deep + 0.3 * craton * basalt * deep,
+        "nickel": basalt * hot * wet + (craton + suture) * deep,
+        "cobalt": (basalt * hot * wet + craton * deep + arc * 0.3) * deep,
+        "vanadium": (basalt + craton) * deep,
+        "platinum": craton * (basalt + meta + 0.3) * deep,
+        "mercury": volc * (0.5 + 0.5 * b(f["t_mean"] > -30)),
+        "sulfur": volc + 0.5 * b(f["salt"]) + 0.3 * sed * arid,
+        "saltpetre": arid * b(lf != LF_MOUNTAINS) + 0.3 * lime * b(p > 800),
+        "bauxite": hot * wet * (gran + basalt + 0.3) * b(rel < 800) + 0.4 * lime * b(t > 14) * b(p > 600),
+        "fluorite": (gran + lime * 0.4) * (rift + oro + 0.2),
+        "rock_salt": b(f["salt"]) + 0.6 * b(f["endo"]) + 0.3 * sed * b(p < 500),
+        "potash": b(f["salt"]) + 0.3 * b(f["endo"]) + 0.2 * sed * b(p < 400),
+        "phosphate": sed * low * b(f["dist_ocean"] < 400) + 0.3 * lime,
+        "oil_gas": sed * low * (1 + b(f["d_coll"] < 1200) + rift),
+        "helium": sed * low * (craton + gran),
+        "diamond": craton * (1 + deep),
+        "magnesite": (lime + meta * 0.5) * (0.3 + deep) + 0.3 * b(f["salt"]),
+        "kaolin": (gran + sand * 0.4) * b(p > 900) * b(t > 8),
+        "glass_sand": sand * (1 + b(lf == LF_DUNES) + b(f["dist_ocean"] < 150)),
+        "graphite": meta * (0.3 + deep),
+    }
+
+
+def place(g, f, seed):
+    """(deposits uint32 bitmask, deposit_main uint8 legend index) over the grid."""
+    land = np.asarray(f["land"], bool)
+    R = rules(f)
+    prov = provinces(g, seed)
+    bits = np.zeros(g.n, np.uint32)
+    main = np.zeros(g.n, np.uint8)
+    best = np.full(g.n, np.inf)
+    for i, (name, _, share) in enumerate(DEPOSITS):
+        nz = 0.5 + 0.5 * fbm(g.xyz, name_seed(seed, "deposit:" + name), 3, 6.0)       # mining districts
+        sel = _rank_select(R[name] * (0.25 + nz) ** 2, share, land, prov)
+        bits[sel] |= np.uint32(1 << i)
+        rarer = sel & (share < best)
+        main[rarer], best[rarer] = i + 1, share
+    return bits, main
diff --git a/mapgen/newworld.py b/mapgen/newworld.py
new file mode 100644
index 0000000..80109b0
--- /dev/null
+++ b/mapgen/newworld.py
@@ -0,0 +1,117 @@
+"""`mapgen.py --world DIR new-world`: a random starting world (sketch/*.png continents, config/tectonics.toml plates),
+so a build runs without drawing anything. Used for the moons; edit or redraw afterwards."""
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+from scipy import ndimage
+
+from .noise import fbm, ridged
+from .sphere import latlon_to_xyz
+
+W, H = 2000, 1000
+
+
+def _spread(rng, n, min_deg, avoid=(), max_lat=70.0, tries=4000):
+    pts = list(avoid)
+    out = []
+    for _ in range(tries):
+        if len(out) == n:
+            break
+        lat = np.degrees(np.arcsin(rng.uniform(np.sin(np.radians(-max_lat)), np.sin(np.radians(max_lat)))))
+        lon = rng.uniform(-180, 180)
+        p = latlon_to_xyz(lat, lon)
+        if all(np.degrees(np.arccos(np.clip(p @ latlon_to_xyz(*q), -1, 1))) >= min_deg for q in pts):
+            pts.append((lat, lon))
+            out.append((float(lat), float(lon)))
+    return out
+
+
+def _centroid(mask, LAT, LON):
+    w = np.cos(np.radians(LAT)) * mask
+    p = (latlon_to_xyz(LAT, LON) * w[..., None]).sum(axis=(0, 1))
+    p /= np.linalg.norm(p)
+    return float(np.degrees(np.arcsin(p[2]))), float(np.degrees(np.arctan2(p[1], p[0])))
+
+
+def make(root: Path, seed: int, continents: int = 6, land_share: float = 0.27, toward=None, ridge: float = 0.86,
+         force: bool = False) -> None:
+    root = Path(root)
+    sk, cfg = root / "sketch", root / "config"
+    targets = [sk / "land.png", cfg / "tectonics.toml"]
+    if not force and any(p.exists() for p in targets):
+        raise SystemExit("new-world: sketch/ or config/tectonics.toml already exist (use --force to overwrite)")
+    rng = np.random.default_rng(seed)
+    w, h = W // 2, H // 2
+    LAT, LON = np.meshgrid(90.0 - (np.arange(h) + 0.5) * 180.0 / h, -180.0 + (np.arange(w) + 0.5) * 360.0 / w,
+                           indexing="ij")
+    xyz = latlon_to_xyz(LAT, LON).reshape(-1, 3)
+
+    centres = _spread(rng, continents, 35.0, max_lat=55.0)
+    field = np.zeros(len(xyz))
+    for lat, lon in centres:
+        r = np.radians(rng.uniform(18.0, 34.0))
+        d = np.arccos(np.clip(xyz @ latlon_to_xyz(lat, lon), -1, 1))
+        field = np.maximum(field, np.clip(1.0 - d / r, 0.0, None) * rng.uniform(0.8, 1.2))
+    field += 0.55 * fbm(xyz, seed + 7, 6, 2.2)
+    if toward is not None:
+        field += 1.2 * (xyz @ latlon_to_xyz(*toward))
+    area = np.cos(np.radians(LAT)).ravel()
+    order = np.argsort(-field)
+    cut = field[order][np.searchsorted(np.cumsum(area[order]) / area.sum(), land_share)]
+    land = (field > cut).reshape(h, w)
+
+    mount = (ridged(xyz, seed + 11, 5, 3.0) > ridge).reshape(h, w) & ndimage.binary_erosion(land, iterations=6)
+
+    def save(name, m):
+        Image.fromarray((m * 255).astype(np.uint8), "L").resize((W, H), Image.BILINEAR).point(
+            lambda v: 255 if v > 127 else 0).save(sk / f"{name}.png")
+
+    sk.mkdir(parents=True, exist_ok=True)
+    cfg.mkdir(parents=True, exist_ok=True)
+    save("land", land)
+    save("mountains", mount)
+    for name in ("desert", "rainforest", "trench"):
+        save(name, np.zeros_like(land))
+
+    lab, n = ndimage.label(land)
+    for a, b in zip(lab[:, 0], lab[:, -1]):
+        if a and b and a != b:
+            lab[lab == b] = a
+    sizes = sorted(((int((lab == i).sum()), i) for i in np.unique(lab) if i), reverse=True)
+    plates, seeds = [], []
+    for k, (px, i) in enumerate(sizes[:continents]):
+        m = lab == i
+        if px > 0.035 * m.size:                                    # a big landmass: two plates, a collision belt between
+            lat, lon = _centroid(m, LAT, LON)
+            ang = rng.uniform(0, np.pi)
+            side = (np.sin(ang) * (LAT - lat) + np.cos(ang) * ((LON - lon + 180) % 360 - 180) * np.cos(np.radians(lat))) > 0
+            parts = [m & side, m & ~side]
+        else:
+            parts = [m]
+        for j, part in enumerate(parts):
+            if part.sum() < 50:
+                continue
+            lat, lon = _centroid(part, LAT, LON)
+            seeds.append((lat, lon))
+            plates.append((f"continent-{k + 1}{'ab'[j] if len(parts) > 1 else ''}", lat, lon, "continental",
+                           rng.uniform(0, 360), rng.uniform(2.0, 5.0)))
+    for k, (lat, lon) in enumerate(_spread(rng, 14 - len(plates) // 2, 28.0, seeds, max_lat=85.0)):
+        plates.append((f"ocean-{k + 1}", lat, lon, "oceanic", rng.uniform(0, 360), rng.uniform(3.0, 8.0)))
+
+    ocean = np.argwhere(~ndimage.binary_dilation(land, iterations=10))
+    def sea_point():
+        y, x = ocean[rng.integers(len(ocean))]
+        return float(LAT[y, x]), float(LON[y, x])
+    f = lambda v: f"[{v[0]:.1f}, {v[1]:.1f}]"
+    t = [f"# Generated by `mapgen.py new-world --seed {seed}`. Edit freely (see README.md → Tectonics).",
+         "# seed = [lat, lon] where the plate grows from; motion = [azimuth° clockwise from north, speed cm/yr].", ""]
+    for pid, lat, lon, kind, az, sp in plates:
+        t += ["[[plate]]", f'id = "{pid}"', f"seed = {f((lat, lon))}", f'kind = "{kind}"', f"motion = [{az:.0f}.0, {sp:.1f}]", ""]
+    if len(ocean):
+        t += ["[[hotspot]]", 'name = "island-chain"', f"center = {f(sea_point())}", "length_km = 600.0", ""]
+    (cfg / "tectonics.toml").write_text("\n".join(t))
+
+    print(f"new-world: {len(plates)} plates, sketch/ and config/tectonics.toml written (seed {seed})")
diff --git a/mapgen/noise.py b/mapgen/noise.py
new file mode 100644
index 0000000..e52b7f8
--- /dev/null
+++ b/mapgen/noise.py
@@ -0,0 +1,111 @@
+"""Seeded, vectorized 3D value noise and fractal sums on points (N, 3)."""
+from __future__ import annotations
+
+import zlib
+
+import numpy as np
+
+_MASK = (1 << 64) - 1
+_K1 = np.uint64(0x9E3779B185EBCA87)
+_K2 = np.uint64(0xC2B2AE3D27D4EB4F)
+_K3 = np.uint64(0x165667B19E3779F9)
+_M1 = np.uint64(0x94D049BB133111EB)
+
+
+def _hash(ix, iy, iz, seed: int):
+    s = np.uint64((seed * 0x27D4EB2F165667C5 + 0x632BE59BD9B4E019) & _MASK)
+    h = ix.astype(np.uint64) * _K1 ^ iy.astype(np.uint64) * _K2 ^ iz.astype(np.uint64) * _K3 ^ s
+    h ^= h >> np.uint64(31)
+    h *= _M1
+    h ^= h >> np.uint64(29)
+    return (h >> np.uint64(11)).astype(np.float64) / float(1 << 53) * 2.0 - 1.0
+
+
+def value_noise(p, seed: int):
+    if _noise_jit is not None:
+        p = np.asarray(p)
+        if p.dtype == np.float64 and p.ndim == 2 and p.shape[1] == 3:
+            s = np.uint64((seed * 0x27D4EB2F165667C5 + 0x632BE59BD9B4E019) & _MASK)
+            return _noise_jit(np.ascontiguousarray(p), s)
+    return _value_noise(p, seed)
+
+
+def _value_noise(p, seed: int):
+    f = np.floor(p)
+    i = f.astype(np.int64)
+    t = p - f
+    t = t * t * (3.0 - 2.0 * t)
+    out = np.zeros(len(p))
+    for dx in (0, 1):
+        wx = t[:, 0] if dx else 1.0 - t[:, 0]
+        for dy in (0, 1):
+            wy = t[:, 1] if dy else 1.0 - t[:, 1]
+            for dz in (0, 1):
+                wz = t[:, 2] if dz else 1.0 - t[:, 2]
+                out += wx * wy * wz * _hash(i[:, 0] + dx, i[:, 1] + dy, i[:, 2] + dz, seed)
+    return out
+
+
+def _make_noise_jit():
+    """value_noise compiled (numba optional): per point the same integer hash and the same float steps in the same
+    order, so the same values; one call costs microseconds instead of a dozen numpy passes (river meanders call it
+    on a few points at a time). WORLDGEN_NO_JIT=1 turns it off."""
+    import os
+    if os.environ.get("WORLDGEN_NO_JIT"):
+        return None
+    try:
+        import numba
+    except ImportError:
+        return None
+    K1, K2, K3, M1 = _K1, _K2, _K3, _M1
+
+    @numba.njit(cache=True)
+    def noise(p, s):
+        n = p.shape[0]
+        out = np.empty(n)
+        for k in range(n):
+            fx, fy, fz = np.floor(p[k, 0]), np.floor(p[k, 1]), np.floor(p[k, 2])
+            ix, iy, iz = np.int64(fx), np.int64(fy), np.int64(fz)
+            tx, ty, tz = p[k, 0] - fx, p[k, 1] - fy, p[k, 2] - fz
+            tx = tx * tx * (3.0 - 2.0 * tx)
+            ty = ty * ty * (3.0 - 2.0 * ty)
+            tz = tz * tz * (3.0 - 2.0 * tz)
+            acc = 0.0
+            for dx in range(2):
+                wx = tx if dx else 1.0 - tx
+                for dy in range(2):
+                    wy = ty if dy else 1.0 - ty
+                    for dz in range(2):
+                        wz = tz if dz else 1.0 - tz
+                        h = (np.uint64(ix + dx) * K1) ^ (np.uint64(iy + dy) * K2) ^ (np.uint64(iz + dz) * K3) ^ s
+                        h ^= h >> np.uint64(31)
+                        h *= M1
+                        h ^= h >> np.uint64(29)
+                        v = np.float64(h >> np.uint64(11)) / 9007199254740992.0 * 2.0 - 1.0
+                        acc += wx * wy * wz * v
+            out[k] = acc
+        return out
+    return noise
+
+
+_noise_jit = _make_noise_jit()
+
+
+def fbm(xyz, seed: int, octaves: int = 5, freq: float = 2.0, lacunarity: float = 2.0, gain: float = 0.5):
+    total = np.zeros(len(xyz))
+    amp, norm, f = 1.0, 0.0, freq
+    for o in range(octaves):
+        total += amp * value_noise(xyz * f, seed + 1013 * o)
+        norm += amp
+        amp *= gain
+        f *= lacunarity
+    return total / norm
+
+
+def ridged(xyz, seed: int, octaves: int = 5, freq: float = 2.0):
+    return 1.0 - np.abs(fbm(xyz, seed, octaves, freq))
+
+
+def name_seed(seed: int, name: str) -> int:
+    """A per-feature noise seed: the build seed plus a hash of the feature's name (stable across runs)."""
+    return seed + zlib.crc32(name.encode()) % 100000
diff --git a/mapgen/ocean.py b/mapgen/ocean.py
new file mode 100644
index 0000000..67dfaf7
--- /dev/null
+++ b/mapgen/ocean.py
@@ -0,0 +1,166 @@
+"""Ocean circulation for stage `climate`: wind-driven surface currents, sea-surface temperature, Ekman upwelling
+and a productivity index.
+
+Currents: the Stommel model on the sphere, r∇²ψ + βψ_x = k·curl(τ/ρH), solved for the stream function ψ over
+every cell; land is the same fluid with `land_friction`× the friction (Brinkman penalisation), so coasts block the
+flow and islands need no special treatment. u = k × ∇ψ. Western boundary currents come out ≈ r/β wide.
+"""
+from __future__ import annotations
+
+import numpy as np
+from scipy import sparse
+from scipy.sparse import linalg as splinalg
+
+from .graph import bicgstab_jacobi, distance_to, gradient, pmap, smooth_km
+from .grid import rowdot_at
+from .sphere import east_north
+
+DEFAULTS = {
+    "enabled": True, "friction_days": 5.0, "layer_m": 150.0, "stress_k": 1.0, "land_friction": 1000.0,
+    "rho_air": 1.2, "drag": 1.3e-3, "direct_max_cells": 500000,
+    "relax_days": 300.0, "kappa_m2s": 1000.0,
+    "ekman_min_lat": 3.0, "upwell_coast_km": 100.0,
+    "prod_upwell": 0.6, "prod_shelf": 0.4, "prod_mix": 0.3, "upwell_ref_m_yr": 100.0,
+    "prod_front": 0.4, "front_min": 0.5, "front_ref": 1.0,
+}
+RHO_W = 1025.0
+YEAR_S = 3.15576e7
+
+
+def omega(day_hours):
+    return 2.0 * np.pi / (float(day_hours) * 3600.0)
+
+
+def divergence(g, V):
+    """Graph divergence of a tangent field (V per km): (2/k_i) Σ_j V_j·t_ij / d_ij (exact for linear V on a hex
+    lattice; the V_i terms cancel)."""
+    w = rowdot_at(V, g.dst, g.edge_tangents) / g.edge_km
+    return np.bincount(g.src, weights=w, minlength=g.n) * (2.0 / g.counts)
+
+
+def _edge_operator(g, w):
+    """Sparse operator f ↦ Σ_j w_ij (f_j − f_i)."""
+    A = sparse.csr_matrix((w, (g.src, g.dst)), shape=(g.n, g.n))
+    return A - sparse.diags(np.asarray(A.sum(axis=1)).ravel())
+
+
+def wind_stress(wind, P):
+    """Bulk formula τ = k·ρ_air·C_d·|w|·w (N/m²)."""
+    wind = np.asarray(wind, dtype=np.float64)
+    return P["stress_k"] * P["rho_air"] * P["drag"] * np.linalg.norm(wind, axis=1, keepdims=True) * wind
+
+
+def streamfunction(g, ocean, wind, P, day_hours):
+    """ψ (m²/s) of the wind-driven surface flow; one direct sparse solve over every cell."""
+    Om = omega(day_hours)
+    R = g.radius_km * 1e3
+    beta = 2.0 * Om * np.cos(np.radians(g.lat)) / R                    # 1/(m·s)
+    beta30 = 2.0 * Om * np.cos(np.radians(30.0)) / R
+    r = max(1.0 / (P["friction_days"] * 86400.0), beta30 * g.spacing_km * 1e3)   # boundary layer ≥ one cell
+    F = np.where(ocean[:, None], wind_stress(wind, P), 0.0) / (RHO_W * P["layer_m"])
+    curl = divergence(g, np.cross(F, g.xyz)) * 1e-3                    # k·curl F = ∇·(F×k), 1/s²
+    fr = np.where(ocean, 1.0, P["land_friction"])
+    fe = 2.0 / (1.0 / fr[g.src] + 1.0 / fr[g.dst])                    # harmonic mean: coasts count as land
+    L = _edge_operator(g, 4.0 / (g.counts[g.src] * g.edge_km ** 2) * fe)          # ∇·(fr∇), per km²
+    e, _ = east_north(g.xyz)
+    Dx = _edge_operator(g, (2.0 / g.counts[g.src]) * rowdot_at(e, g.src, g.edge_tangents) / g.edge_km)
+    A = (L + sparse.diags(beta / r * 1e3) @ Dx).tolil()                # β/r·1e3: per km
+    b = 1e6 * curl / r
+    # gauge: ψ is defined up to a constant; pin it deep inland, where land friction soaks up the solve's
+    # compatibility residual (pinned at sea it would leave a point vortex there)
+    k = int(np.argmax(distance_to(g, ocean))) if (~ocean).any() else 0
+    A[k, :] = 0.0
+    A[k, k] = 1.0
+    b[k] = 0.0
+    return splinalg.spsolve(A.tocsc(), b)
+
+
+def velocity(g, psi):
+    """u = k × ∇ψ (m/s), tangent 3-vectors."""
+    return np.cross(g.xyz, gradient(g, psi) * 1e-3)
+
+
+def _coarse_currents(g, ocean, wind, P, day_hours):
+    """Grids too big for a direct solve: solve on the parent H3 resolution, carry u down, smooth."""
+    from .grid import build_grid, cell_parents
+    cg = build_grid(g.res - 1, g.radius_km)
+    parent = np.searchsorted(cg.ids, cell_parents(g.ids, g.res - 1))
+    cnt = np.maximum(np.bincount(parent, minlength=cg.n), 1)
+    c_ocean = np.bincount(parent, weights=ocean.astype(float), minlength=cg.n) / cnt > 0.5
+    c_wind = np.stack([np.bincount(parent, weights=wind[:, k], minlength=cg.n) / cnt for k in range(3)], axis=1)
+    cu = currents(cg, c_ocean, c_wind, P, day_hours)
+    u = np.stack(pmap(lambda k: smooth_km(g, cu[parent, k], cg.spacing_km), range(3)), axis=1)
+    return u - np.sum(u * g.xyz, axis=1, keepdims=True) * g.xyz       # back into the tangent plane
+
+
+def currents(g, ocean, wind, P, day_hours):
+    """Surface current (n,3) m/s; 0 on land."""
+    ocean = np.asarray(ocean, bool)
+    wind = np.asarray(wind, dtype=np.float64)
+    if not ocean.any():
+        return np.zeros((g.n, 3))
+    if g.n > P["direct_max_cells"]:
+        u = _coarse_currents(g, ocean, wind, P, day_hours)
+    else:
+        u = velocity(g, streamfunction(g, ocean, wind, P, day_hours))
+    return np.where(ocean[:, None], u, 0.0)
+
+
+def _masked_edges(g, mask, w):
+    return _edge_operator(g, np.where(mask[g.src] & mask[g.dst], w, 0.0))
+
+
+def sst(g, ocean, u, T_eq, P):
+    """Annual-mean sea-surface temperature (°C): steady u·∇T − κ∇²T = λ(T_eq − T) over the ocean (upwind
+    advection, no flux into land); land keeps T_eq."""
+    ocean = np.asarray(ocean, bool)
+    T_eq = np.asarray(T_eq, dtype=np.float64)
+    lam = 1.0 / (P["relax_days"] * 86400.0)
+    up = (4.0 / g.counts[g.src]) * np.maximum(-rowdot_at(u, g.src, g.edge_tangents), 0.0) / (g.edge_km * 1e3)
+    adv = -_masked_edges(g, ocean, up)                                 # Σ_upstream c_ij (T_i − T_j), 1/s
+    dif = _masked_edges(g, ocean, 4.0 / (g.counts[g.src] * g.edge_km ** 2)) * (P["kappa_m2s"] * 1e-6)
+    A = (adv - dif) / lam + sparse.identity(g.n)
+    A = (sparse.diags(ocean.astype(float)) @ A + sparse.diags((~ocean).astype(float))).tocsr()
+    T, info = bicgstab_jacobi(A, T_eq, T_eq.copy(), 1.0 / A.diagonal(), 1e-9, 5000)
+    if info != 0:
+        raise ValueError(f"sst: solver did not converge (info={info})")
+    return T
+
+
+def upwelling(g, ocean, wind, P, day_hours):
+    """Ekman pumping (m/yr, + up): w = ∇·M, M = τ×k/(ρf), |f| floored at `ekman_min_lat`. Transport pointing off a
+    coast leaves the coast cell (land carries none), so coastal upwelling needs no separate rule. Smoothed over
+    `upwell_coast_km`; 0 on land."""
+    ocean = np.asarray(ocean, bool)
+    Om = omega(day_hours)
+    fmin = 2.0 * Om * np.sin(np.radians(P["ekman_min_lat"]))
+    f = np.where(g.lat >= 0, 1.0, -1.0) * np.maximum(np.abs(2.0 * Om * np.sin(np.radians(g.lat))), fmin)
+    tau = np.where(ocean[:, None], wind_stress(wind, P), 0.0)
+    M = np.cross(tau, g.xyz) / (RHO_W * f[:, None])                    # m²/s
+    w = np.where(ocean, divergence(g, M) * 1e-3 * YEAR_S, 0.0)
+    return np.where(ocean, smooth_km(g, w, P["upwell_coast_km"]), 0.0)
+
+
+def sst_front(g, ocean, sst_c):
+    """|∇SST| over the sea (°C per 100 km); edges to land count as flat, so coasts are no front."""
+    sst_c = np.asarray(sst_c, dtype=np.float64)
+    both = ocean[g.src] & ocean[g.dst]
+    w = np.where(both, sst_c[g.dst] - sst_c[g.src], 0.0) / g.edge_km
+    grad = np.stack([np.bincount(g.src, weights=w * g.edge_tangents[:, c], minlength=g.n) for c in range(3)], axis=1)
+    return np.where(ocean, np.linalg.norm(grad, axis=1) * (2.0 / g.counts) * 100.0, 0.0)
+
+
+def productivity(g, ocean, w, z, t_range, sst_c, P):
+    """0–1 sea productivity: upwelling (saturating), shallow shelf, winter mixing and SST fronts (where warm and cold
+    currents meet, e.g. a Brazil–Malvinas confluence; gradients under `front_min` °C/100 km add nothing), dimmed
+    toward the poles."""
+    ocean = np.asarray(ocean, bool)
+    depth = np.maximum(-np.asarray(z, dtype=np.float64), 0.0)
+    x = np.maximum(w, 0.0) / P["upwell_ref_m_yr"]
+    shelf = np.clip((1000.0 - depth) / 800.0, 0.0, 1.0)
+    mix = np.clip(np.asarray(t_range) / 20.0, 0.0, 1.0) * np.clip((20.0 - np.asarray(sst_c)) / 20.0, 0.0, 1.0)
+    light = 0.3 + 0.7 * np.cos(np.radians(g.lat))
+    xf = np.maximum(sst_front(g, ocean, sst_c) - P["front_min"], 0.0) / P["front_ref"]
+    p = (P["prod_upwell"] * x / (1.0 + x) + P["prod_shelf"] * shelf + P["prod_mix"] * mix
+         + P["prod_front"] * xf / (1.0 + xf)) * light
+    return np.where(ocean, np.clip(p, 0.0, 1.0), 0.0)
diff --git a/mapgen/pipeline.py b/mapgen/pipeline.py
new file mode 100644
index 0000000..fca2fa1
--- /dev/null
+++ b/mapgen/pipeline.py
@@ -0,0 +1,180 @@
+"""Stage runner with an input-hash cache."""
+from __future__ import annotations
+
+import hashlib
+import importlib
+import os
+import time
+from dataclasses import dataclass, field
+from pathlib import Path
+
+import numpy as np
+
+from . import config as C
+from .grid import Grid
+
+STAGES = ["grid", "sketch", "plates", "crust", "elevation", "erosion", "climate",
+          "hydrology", "seabed", "environment", "ice", "fields", "render"]
+
+
+class StageError(RuntimeError):
+    pass
+
+
+@dataclass
+class Ctx:
+    root: Path
+    cfg: dict
+    tect: dict
+    res: int
+    data: dict = field(default_factory=dict)
+    grid: Grid | None = None
+    out_dir: Path | None = None        # render: out/r<res> unless set (an era writes out/r<res>/eras/<era>)
+    preview_dir: Path | None = None
+    low_memory: bool = False           # trade speed for a lower memory peak; never changes the results
+
+    @property
+    def seed(self) -> int:
+        return int(self.cfg["build"]["seed"])
+
+    def need(self, *keys):
+        missing = [k for k in keys if k not in self.data]
+        if missing:
+            raise StageError(f"missing inputs {missing}: run the stage that produces them first")
+        return [self.data[k] for k in keys]
+
+
+def inputs_key(root: Path, res: int) -> str:
+    h = hashlib.sha256(str(res).encode())
+    files = []
+    for sub, pat in (("config", "*.toml"), ("masks", "*.png"), ("sketch", "*.png")):
+        files += sorted(p for p in (root / sub).glob(pat) if p.name != C.ERAS_FILE)   # eras: their own key
+    files += sorted(Path(__file__).resolve().parent.glob("*.py"))
+    for p in files:
+        h.update(p.name.encode())
+        h.update(p.read_bytes())
+    return h.hexdigest()[:16]
+
+
+ALIGN = 64
+
+
+def save_npz_aligned(path, /, **arrays) -> None:
+    """np.savez(path, **arrays), but each array's data starts at a multiple of ALIGN bytes in the file (the local zip
+    header gets a padding extra field, as zipalign does): an ordinary .npz for np.load, whose arrays npz_maps can
+    map in place. Alignment matters for more than speed: numpy sums 8-byte-misaligned data in buffered chunks,
+    i.e. in another order — mapped misaligned inputs would change the last bits of results."""
+    import io
+    import struct
+    import zipfile
+    from numpy.lib import format as F
+    with zipfile.ZipFile(path, "w", compression=zipfile.ZIP_STORED, allowZip64=True) as zf:
+        for name, v in arrays.items():
+            v = np.asarray(v)
+            if v.dtype.hasobject:
+                raise ValueError(f"save_npz_aligned: {name} has object dtype")
+            head = io.BytesIO()
+            d = F.header_data_from_array_1_0(v)
+            try:
+                F.write_array_header_1_0(head, d)
+            except ValueError:
+                head = io.BytesIO()
+                F.write_array_header_2_0(head, d)
+            info = zipfile.ZipInfo(f"{name}.npy", date_time=(1980, 1, 1, 0, 0, 0))
+            info.compress_type = zipfile.ZIP_STORED
+            start = zf.fp.tell() + 30 + len(info.filename.encode()) + 20 + len(head.getvalue())   # 20: zip64 field
+            pad = -start % ALIGN
+            if 0 < pad < 4:                                          # an extra field is at least its 4-byte header
+                pad += ALIGN
+            if pad:
+                info.extra = struct.pack("<HH", 0xA1A1, pad - 4) + bytes(pad - 4)
+            with zf.open(info, "w", force_zip64=True) as m:
+                m.write(head.getvalue())
+                _write_data(m, v)
+
+
+def _write_data(m, v) -> None:
+    """The array's bytes (C order, or Fortran order for an F-contiguous array, as np.save) in 16 MB pieces."""
+    flat = v.T.reshape(-1) if (v.flags.f_contiguous and not v.flags.c_contiguous) else np.ascontiguousarray(v).reshape(-1)
+    step = max(1, (16 << 20) // max(v.itemsize, 1))
+    for i in range(0, flat.size, step):
+        m.write(flat[i:i + step].tobytes())
+
+
+def npz_maps(f: Path) -> dict:
+    """The arrays of an uncompressed .npz (np.savez) mapped from the file, copy-on-write: plain writable arrays with
+    the same values, whose pages the OS reads on use and can drop again (writes stay private, the file never
+    changes). Members that can't be mapped (compressed, object dtype) are read into memory as np.load would."""
+    import mmap
+    import struct
+    import zipfile
+    out = {}
+    with open(f, "rb") as fh, zipfile.ZipFile(fh) as zf:
+        mm = mmap.mmap(fh.fileno(), 0, access=mmap.ACCESS_COPY)
+        for info in zf.infolist():
+            name = info.filename[:-4] if info.filename.endswith(".npy") else info.filename
+            ok = info.compress_type == zipfile.ZIP_STORED
+            if ok:
+                fh.seek(info.header_offset)
+                local = fh.read(30)
+                n_name, n_extra = struct.unpack("<HH", local[26:30])
+                fh.seek(info.header_offset + 30 + n_name + n_extra)
+                version = np.lib.format.read_magic(fh)
+                read_header = {(1, 0): np.lib.format.read_array_header_1_0,
+                               (2, 0): np.lib.format.read_array_header_2_0}.get(version)
+                ok = read_header is not None
+            if ok:
+                shape, fortran, dtype = read_header(fh, max_header_size=1 << 20)
+                ok = not dtype.hasobject
+                ok = fh.tell() % ALIGN == 0                         # misaligned: read (see save_npz_aligned)
+            if ok:
+                count = int(np.prod(shape, dtype=np.int64))
+                a = np.frombuffer(mm, dtype=dtype, count=count, offset=fh.tell()) if count else np.empty(0, dtype)
+                out[name] = a.reshape(shape, order="F" if fortran else "C")
+            else:
+                with zf.open(info) as m:
+                    out[name] = np.lib.format.read_array(m, allow_pickle=False)
+    return out
+
+
+def _after(ctx: Ctx, name: str) -> None:
+    if name == "grid":
+        ctx.grid = Grid.from_arrays(ctx.data, ctx.res, float(ctx.cfg["planet"]["radius_km"]))
+
+
+def build(root: Path, res: int, start: str | None = None, stop: str | None = None, log=print,
+          low_memory: bool = False) -> Ctx:
+    cfg, tect = C.load(root)
+    ctx = Ctx(root, cfg, tect, res, low_memory=low_memory)
+    cache = root / "out" / "cache" / f"r{res}"
+    cache.mkdir(parents=True, exist_ok=True)
+    key = inputs_key(root, res)
+    stages = STAGES[: STAGES.index(stop) + 1] if stop else STAGES
+    forced = False
+    for name in stages:
+        forced = forced or name == start
+        f = cache / f"{name}.npz"
+        t0 = time.time()
+        if not forced and f.exists():
+            with np.load(f, allow_pickle=False) as z:
+                hit = str(z["_key"]) == key
+                if hit and not ctx.low_memory:
+                    ctx.data.update({k: z[k] for k in z.files if k != "_key"})
+            if hit:
+                if ctx.low_memory:
+                    ctx.data.update({k: v for k, v in npz_maps(f).items() if k != "_key"})
+                _after(ctx, name)
+                log(f"{name}: cached")
+                continue
+        forced = True
+        out = importlib.import_module(f"mapgen.{name}").run(ctx)
+        tmp = f.with_name(f"{f.stem}.tmp-{os.getpid()}.npz")
+        save_npz_aligned(tmp, _key=np.array(key), **out)
+        os.replace(tmp, f)                  # never truncated in place: maps of the old file stay valid
+        if ctx.low_memory:                  # fields kept on disk from here on (the cache file just written)
+            maps = npz_maps(f)
+            out = {k: maps.get(k, v) if isinstance(v, np.ndarray) else v for k, v in out.items()}
+        ctx.data.update(out)
+        _after(ctx, name)
+        log(f"{name}: {time.time() - t0:.1f}s")
+    return ctx
diff --git a/mapgen/plateaus.py b/mapgen/plateaus.py
new file mode 100644
index 0000000..9647d04
--- /dev/null
+++ b/mapgen/plateaus.py
@@ -0,0 +1,179 @@
+"""Sunken plateaus: Kerguelen-type continental crust that never rose. Outline,
+surface and volcanic field are functions of position and the plateau's config (deterministic per name), so the world
+build, the seabed pass and the viewer agree at any resolution."""
+from __future__ import annotations
+
+import numpy as np
+
+from .graph import distance_to
+from .noise import fbm, name_seed
+from .sphere import east_north, gc_dist_km, great_circle_point, latlon_to_xyz, tangent_dir, xyz_to_latlon
+from .zones import smootherstep
+
+EDGE_WARP = 0.15        # outline radius ± 15 % (fractal noise): bays and lobes, like a small continent
+MARGIN_KM = 150.0       # the surface blends down to the surrounding sea floor over this distance inside the outline
+RELIEF_M = 400.0        # internal relief (ridges, basins)
+HIDDEN_MAX_M = -550.0   # hidden plateaus: every cell at least this deep (checked at −500 m after erosion's re-solve)
+SITE_RULES = (("land", 400.0), ("a plate boundary", 150.0), ("the trench", 200.0))   # km from the outline
+
+
+def semi_axes(p: dict):
+    """(long, short) semi-axes (km) of the plateau's ellipse: π·a·b = area_km2, a / b = elongation."""
+    e = float(p.get("elongation", 1.0))
+    b = float(np.sqrt(p["area_km2"] / (np.pi * e)))
+    return e * b, b
+
+
+def _near(xyz, p, radius_km, pad_km=0.0):
+    a, _ = semi_axes(p)
+    lim = min(np.pi, (a * (1.0 + EDGE_WARP) / (1.0 - EDGE_WARP) + pad_km) / radius_km)
+    return xyz @ latlon_to_xyz(*p["center"]) > np.cos(lim)
+
+
+def rho(xyz, p: dict, seed: int, radius_km: float):
+    """Normalised radius of points: 0 at the centre, < 1 inside the noise-warped elliptical outline."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    c = latlon_to_xyz(*p["center"])
+    e, n = east_north(c[None])
+    ang = np.arccos(np.clip(xyz @ c, -1.0, 1.0)) * radius_km            # km from the centre along the surface
+    t = tangent_dir(np.broadcast_to(c, xyz.shape), xyz)
+    x, y = ang * (t @ e[0]), ang * (t @ n[0])                            # east, north (azimuthal equidistant)
+    az = np.radians(p.get("azimuth_deg", 0.0))
+    along, across = x * np.sin(az) + y * np.cos(az), x * np.cos(az) - y * np.sin(az)
+    a, b = semi_axes(p)
+    r = np.sqrt((along / a) ** 2 + (across / b) ** 2)
+    w = np.clip(2.0 * fbm(xyz, name_seed(seed, p["name"]), 4, 25.0), -1.0, 1.0)
+    return r / (1.0 + EDGE_WARP * w)
+
+
+def cell_ids(xyz, plateaus: list, seed: int, radius_km: float):
+    """Plateau index per point (−1 outside every plateau)."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    ids = np.full(len(xyz), -1, np.int16)
+    for k, p in enumerate(plateaus):
+        idx = np.flatnonzero(_near(xyz, p, radius_km))
+        if len(idx):
+            ids[idx[rho(xyz[idx], p, seed, radius_km) < 1.0]] = k
+    return ids
+
+
+def _place(rng, p, seed, radius_km, origin, max_km, inside=0.85, tries=30):
+    """A random point within max_km of origin inside the outline (rho < inside); None if none is found."""
+    e, n = east_north(origin[None])
+    for _ in range(tries):
+        az, dist = rng.uniform(0.0, 2.0 * np.pi), rng.uniform(0.0, max_km)
+        q = great_circle_point(origin, np.sin(az) * e[0] + np.cos(az) * n[0], dist, radius_km)
+        if rho(q[None], p, seed, radius_km)[0] < inside:
+            return q
+    return None
+
+
+def _ll(q):
+    lat, lon = xyz_to_latlon(np.asarray(q, dtype=np.float64))
+    return round(float(lat), 4), round(float(lon), 4)
+
+
+def features(p: dict, seed: int, radius_km: float) -> dict:
+    """The plateau's volcanic field at world scale (deterministic per name): its own hotspot point, 6–20 cones and
+    1–3 calderas (none when vent = 0); island plateaus lift 3–6 of the cones to +0.6…+1.5 km."""
+    rng = np.random.default_rng(name_seed(seed, p["name"]) + 7)
+    c = latlon_to_xyz(*p["center"])
+    _, b = semi_axes(p)
+    vent = float(p.get("vent", 1.0))
+    hot = _place(rng, p, seed, radius_km, c, 0.3 * b)
+    hot = c if hot is None else hot
+    cones, calderas = [], []
+    if vent > 0:
+        for _ in range(int(rng.integers(6, 21))):
+            q = _place(rng, p, seed, radius_km, hot, 0.6 * b)
+            r_km, h = float(rng.uniform(15.0, 40.0)), float(rng.uniform(500.0, 2000.0)) * vent
+            if q is not None:
+                lat, lon = _ll(q)
+                cones.append({"lat": lat, "lon": lon, "radius_km": r_km, "height_m": h, "island": False})
+        for _ in range(int(rng.integers(1, 4))):
+            q = _place(rng, p, seed, radius_km, hot, 0.5 * b)
+            r_km = float(rng.uniform(20.0, 50.0))
+            if q is not None:
+                lat, lon = _ll(q)
+                calderas.append({"lat": lat, "lon": lon, "radius_km": r_km, "rim_m": 300.0 * vent,
+                                 "floor_m": -400.0 * vent})
+    if p.get("islands", False) and cones:
+        k = min(len(cones), int(rng.integers(3, 7)))
+        for i in rng.choice(len(cones), size=k, replace=False):
+            cones[int(i)].update(island=True, peak_m=float(rng.uniform(600.0, 1500.0)),
+                                 radius_km=float(rng.uniform(40.0, 70.0)))
+    return {"name": p["name"], "center": [float(v) for v in p["center"]], "hotspot": list(_ll(hot)),
+            "cones": cones, "calderas": calderas}
+
+
+def surface(xyz, p: dict, feat: dict, seed: int, radius_km: float, z_floor):
+    """Heights of points inside the outline: the plateau top (depth across top_m by low-frequency noise, ± RELIEF_M)
+    plus its cones and calderas, blended down to z_floor over MARGIN_KM inside the edge."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    s = name_seed(seed, p["name"])
+    t = np.clip(0.5 + 1.5 * fbm(xyz, s + 1, 3, 8.0), 0.0, 1.0)
+    top0, top1 = p["top_m"]
+    z = -(top0 + (top1 - top0) * t) + RELIEF_M * np.clip(2.5 * fbm(xyz, s + 2, 5, 40.0), -1.0, 1.0)
+    for c in feat["cones"]:
+        if not c["island"]:
+            d = gc_dist_km(xyz, latlon_to_xyz(c["lat"], c["lon"]), radius_km)
+            z = z + c["height_m"] * np.clip(1.0 - d / c["radius_km"], 0.0, 1.0) ** 1.5
+    for c in feat["calderas"]:
+        d = gc_dist_km(xyz, latlon_to_xyz(c["lat"], c["lon"]), radius_km)
+        r = c["radius_km"]
+        z = z + c["rim_m"] * np.exp(-((d - r) / (0.25 * r)) ** 2) + c["floor_m"] * (1.0 - smootherstep(d / r))
+    _, b = semi_axes(p)
+    w = smootherstep((1.0 - rho(xyz, p, seed, radius_km)) * b / MARGIN_KM)
+    return np.asarray(z_floor, dtype=np.float64) + (z - z_floor) * w
+
+
+def islands(xyz, feat: dict, radius_km: float, z):
+    """Island cones lift the ground to their peak_m (small volcanic islands); only raises."""
+    z = np.asarray(z, dtype=np.float64)
+    for c in feat["cones"]:
+        if c["island"]:
+            d = gc_dist_km(np.asarray(xyz, dtype=np.float64), latlon_to_xyz(c["lat"], c["lon"]), radius_km)
+            f = np.clip(d / c["radius_km"], 0.0, 1.0)
+            z = np.where(d < c["radius_km"], np.maximum(z, c["peak_m"] - (c["peak_m"] - z) * f ** 1.2), z)
+    return z
+
+
+def apply(g, z, plateaus: list, plateau_id, seed: int):
+    """Grid heights with the plateaus set (after the world's sea-level solve): surfaces and volcanic fields; hidden
+    plateaus clamped to HIDDEN_MAX_M; islands, each island's nearest cell raised to its peak (so every resolution
+    keeps it)."""
+    z = np.asarray(z, dtype=np.float64).copy()
+    R = g.radius_km
+    for k, p in enumerate(plateaus):
+        idx = np.flatnonzero(np.asarray(plateau_id) == k)
+        if len(idx) == 0:
+            continue
+        feat = features(p, seed, R)
+        zk = surface(g.xyz[idx], p, feat, seed, R, z[idx])
+        if p.get("islands", False):
+            z[idx] = islands(g.xyz[idx], feat, R, zk)
+            for c in feat["cones"]:
+                if c["island"]:
+                    i = g.cell_index(c["lat"], c["lon"])
+                    z[i] = max(z[i], c["peak_m"])
+        else:
+            z[idx] = np.minimum(zk, HIDDEN_MAX_M)
+    return z
+
+
+def site_report(g, data: dict, plateaus: list) -> list:
+    """Plateaus that no longer fit their site on this world, as warnings: outline ≥ 400 km from land,
+    ≥ 150 km from plate boundaries, ≥ 200 km from the sketch's trench. Reported, never moved."""
+    ids = np.asarray(data["plateau_id"])
+    masks = {"land": ~np.asarray(data["ocean"]) & (ids < 0), "a plate boundary": np.asarray(data["bnd_type"]) > 0,
+             "the trench": np.asarray(data.get("sk_trench", np.zeros(g.n))) > 0.5}
+    out = []
+    for what, lim in SITE_RULES:
+        if not masks[what].any():
+            continue
+        d = distance_to(g, masks[what])
+        for k, p in enumerate(plateaus):
+            m = ids == k
+            if m.any() and float(d[m].min()) < lim:
+                out.append(f"plateau {p['name']}: {float(d[m].min()):.0f} km from {what} (rule ≥ {lim:.0f} km)")
+    return out
diff --git a/mapgen/plates.py b/mapgen/plates.py
new file mode 100644
index 0000000..af798bf
--- /dev/null
+++ b/mapgen/plates.py
@@ -0,0 +1,86 @@
+"""Stage `plates`: grow plates from seeds, assign Euler motions, classify boundaries."""
+from __future__ import annotations
+
+import numpy as np
+from scipy.spatial import cKDTree
+
+from .config import params
+from .graph import nearest_source
+from .grid import edge_blocks
+from .noise import fbm
+from .pipeline import StageError
+from .sphere import motion_to_omega, velocity
+
+NONE, CONV, DIV, TRANS = 0, 1, 2, 3
+DEFAULTS = {"land_cost": 0.25, "noise": 0.35, "noise_freq": 3.0, "warp_km": 1000.0, "warp_freq": 1.5,
+            "min_rate_m_yr": 0.005}
+
+
+def edge_convergence(g, vel):
+    """Per directed edge s→d: closing speed (m/yr, >0 converging) and tangential slip speed."""
+    t, src, dst = g.edge_tangents, g.src, g.dst
+    along, tang = np.empty(len(dst)), np.empty(len(dst))
+    for s in edge_blocks(len(dst)):                      # per edge: the same values, small temporaries
+        rel = vel[dst[s]] - vel[src[s]]
+        along[s] = np.sum(rel * t[s], axis=1)
+        tang[s] = np.linalg.norm(rel - along[s][:, None] * t[s], axis=1)
+    return -along, tang
+
+
+def classify_boundaries(g, plate, vel, min_rate):
+    conv, tang = edge_convergence(g, vel)
+    s, d = g.src, g.dst
+    b = plate[s] != plate[d]
+    cnt = np.bincount(s[b], minlength=g.n)
+    c_mean = np.bincount(s[b], weights=conv[b], minlength=g.n) / np.maximum(cnt, 1)
+    t_mean = np.bincount(s[b], weights=tang[b], minlength=g.n) / np.maximum(cnt, 1)
+    other = np.full(g.n, -1, np.int16)
+    other[s[b]] = plate[d[b]]
+    typ = np.zeros(g.n, np.int8)
+    isb = cnt > 0
+    typ[isb] = TRANS
+    typ[isb & (c_mean > min_rate)] = CONV
+    typ[isb & (c_mean < -min_rate)] = DIV
+    rate = np.where(typ == CONV, c_mean, np.where(typ == DIV, -c_mean, np.where(isb, t_mean, 0.0)))
+    return typ, rate.astype(np.float32), other
+
+
+def _warp_labels(g, plate, seeds, P, seed, continental_kind):
+    """Domain warp: each cell takes the label found at a noise-displaced point → meandering boundaries.
+    Only same-kind swaps (cont↔cont, ocean↔ocean): ocean–continent margins stay where growth put them."""
+    if P["warp_km"] <= 0:
+        return plate
+    amp = P["warp_km"] / g.radius_km
+    disp = np.stack([fbm(g.xyz, seed + 101 + k, 3, P["warp_freq"]) for k in range(3)], axis=1)
+    disp /= max(float(disp.std()), 1e-12)              # warp_km = RMS displacement per component
+    q = g.xyz + amp * disp
+    q /= np.linalg.norm(q, axis=1, keepdims=True)
+    _, j = cKDTree(g.xyz).query(q)
+    out = plate[j].astype(np.int16)
+    keep = continental_kind[out] != continental_kind[plate]
+    out[keep] = plate[keep]
+    out[seeds] = np.arange(len(seeds), dtype=np.int16)
+    return out
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "plates", DEFAULTS)
+    land, land_hint = ctx.need("sk_land", "m_land_hint")
+    plates = ctx.tect["plate"]
+    seeds = np.array([g.cell_index(*p["seed"]) for p in plates], dtype=np.int64)
+    if len(np.unique(seeds)) < len(seeds):
+        raise StageError("two plate seeds fall in the same cell; move one in tectonics.toml")
+    hint = np.clip(land + 0.5 * land_hint, 0.0, 1.0)
+    nz = fbm(g.xyz, ctx.seed + 11, octaves=4, freq=P["noise_freq"])
+    mult = np.where(hint[g.dst] > 0.5, P["land_cost"], 1.0) * (1.0 + P["noise"] * nz[g.dst])
+    _, src = nearest_source(g, seeds, g.edge_km * np.maximum(mult, 0.05))
+    order = np.argsort(seeds)
+    plate = order[np.searchsorted(seeds[order], src)].astype(np.int16)
+    kinds = np.array([p["kind"] == "continental" for p in plates])
+    plate = _warp_labels(g, plate, seeds, P, ctx.seed, kinds)
+    omegas = np.array([motion_to_omega(*p["seed"], *p["motion"], g.radius_km) for p in plates])
+    vel = velocity(g.xyz, omegas[plate], g.radius_km)
+    btype, brate, bother = classify_boundaries(g, plate, vel, P["min_rate_m_yr"])
+    return {"plate": plate, "vel": vel, "plate_continental": kinds,
+            "bnd_type": btype, "bnd_rate": brate, "bnd_other": bother}
diff --git a/mapgen/projections.py b/mapgen/projections.py
new file mode 100644
index 0000000..67efb19
--- /dev/null
+++ b/mapgen/projections.py
@@ -0,0 +1,150 @@
+"""Map projections of the equirectangular relief: Equal Earth, Mollweide (equal-area) and orthographic globes."""
+from __future__ import annotations
+
+import numpy as np
+from PIL import Image, ImageDraw
+
+from .sketch import sample_equirect
+from .sphere import east_north, latlon_to_xyz, xyz_to_latlon
+
+BACKGROUND = (16, 18, 24)
+_A1, _A2, _A3, _A4 = 1.340264, -0.081106, 0.000893, 0.003796
+_M = np.sqrt(3.0) / 2.0
+EE_XMAX = 2.0 * np.sqrt(3.0) * np.pi / (3.0 * _A1)
+EE_YMAX = _A1 * np.pi / 3 + _A2 * (np.pi / 3) ** 3 + _A3 * (np.pi / 3) ** 7 + _A4 * (np.pi / 3) ** 9
+MW_XMAX, MW_YMAX = 2.0 * np.sqrt(2.0), np.sqrt(2.0)
+
+
+def equal_earth_forward(lat, lon):
+    t = np.arcsin(_M * np.sin(np.radians(lat)))
+    d = _A1 + 3 * _A2 * t**2 + 7 * _A3 * t**6 + 9 * _A4 * t**8
+    x = 2 * np.sqrt(3.0) * np.radians(lon) * np.cos(t) / (3 * d)
+    y = _A1 * t + _A2 * t**3 + _A3 * t**7 + _A4 * t**9
+    return x, y
+
+
+def equal_earth_inverse(x, y):
+    t = np.asarray(y, dtype=np.float64) / _A1
+    for _ in range(12):
+        f = _A1 * t + _A2 * t**3 + _A3 * t**7 + _A4 * t**9 - y
+        t = t - f / (_A1 + 3 * _A2 * t**2 + 7 * _A3 * t**6 + 9 * _A4 * t**8)
+    d = _A1 + 3 * _A2 * t**2 + 7 * _A3 * t**6 + 9 * _A4 * t**8
+    lat = np.degrees(np.arcsin(np.clip(np.sin(t) / _M, -1, 1)))
+    lon = np.degrees(3 * np.asarray(x) * d / (2 * np.sqrt(3.0) * np.cos(t)))
+    return lat, lon
+
+
+def mollweide_forward(lat, lon):
+    phi = np.radians(lat)
+    t = phi.copy() if isinstance(phi, np.ndarray) else np.array(phi, dtype=np.float64)
+    for _ in range(30):
+        f = 2 * t + np.sin(2 * t) - np.pi * np.sin(phi)
+        t = t - f / np.maximum(2 + 2 * np.cos(2 * t), 1e-12)
+    return MW_XMAX / np.pi * np.radians(lon) * np.cos(t), MW_YMAX * np.sin(t)
+
+
+def mollweide_inverse(x, y):
+    t = np.arcsin(np.clip(np.asarray(y) / MW_YMAX, -1, 1))
+    lat = np.degrees(np.arcsin(np.clip((2 * t + np.sin(2 * t)) / np.pi, -1, 1)))
+    lon = np.degrees(np.pi * np.asarray(x) / (MW_XMAX * np.maximum(np.cos(t), 1e-12)))
+    return lat, lon
+
+
+def orthographic_inverse(x, y, lat0, lon0):
+    """Unit-disc view coords → lat/lon on the visible hemisphere centred on (lat0, lon0); ok = inside the disc."""
+    x, y = np.asarray(x, dtype=np.float64), np.asarray(y, dtype=np.float64)
+    rho2 = x**2 + y**2
+    ok = rho2 <= 1.0
+    z = np.sqrt(np.clip(1.0 - rho2, 0.0, 1.0))
+    c = latlon_to_xyz(np.array([lat0]), np.array([lon0]))
+    e, n = east_north(c)
+    p = x[..., None] * e[0] + y[..., None] * n[0] + z[..., None] * c[0]
+    lat, lon = xyz_to_latlon(p)
+    return lat, lon, ok
+
+
+def _sample_rgb(img, lat, lon):
+    """Bilinear RGB samples, read straight from the uint8 channels (each sampled value promotes to float64 exactly,
+    so no full-size float copy of a big image is needed)."""
+    return np.stack([sample_equirect(img[..., k], lat, lon) for k in range(3)], axis=-1)
+
+
+def reproject(img, proj, width):
+    """Equirectangular RGB → equal-area world map ('equal_earth' | 'mollweide')."""
+    xmax, ymax, inv = {"equal_earth": (EE_XMAX, EE_YMAX, equal_earth_inverse),
+                       "mollweide": (MW_XMAX, MW_YMAX, mollweide_inverse)}[proj]
+    height = int(round(width * ymax / xmax))
+    xs = ((np.arange(width) + 0.5) / width * 2 - 1) * xmax
+    ys = (1 - (np.arange(height) + 0.5) / height * 2) * ymax
+    X, Y = np.meshgrid(xs, ys)
+    lat, lon = inv(X, Y)
+    ok = np.isfinite(lat) & np.isfinite(lon) & (np.abs(lon) <= 180.0)
+    rgb = _sample_rgb(img, np.where(ok, lat, 0.0), np.where(ok, lon, 0.0))
+    out = np.where(ok[..., None], rgb, np.array(BACKGROUND, dtype=np.float64))
+    return np.clip(np.round(out), 0, 255).astype(np.uint8)
+
+
+def globe(img, lat0, lon0, size):
+    """Orthographic view centred on (lat0, lon0), gentle limb darkening."""
+    v = (np.arange(size) + 0.5) / size * 2 - 1
+    X, Y = np.meshgrid(v, -v)
+    lat, lon, ok = orthographic_inverse(X, Y, lat0, lon0)
+    rgb = _sample_rgb(img, lat, lon) * (0.72 + 0.28 * np.sqrt(np.clip(1 - X**2 - Y**2, 0, 1)))[..., None]
+    out = np.where(ok[..., None], rgb, np.array(BACKGROUND, dtype=np.float64))
+    return np.clip(np.round(out), 0, 255).astype(np.uint8)
+
+
+def graticule(arr, proj, step=30):
+    """Draw a light lat/lon grid (every `step`°) onto an equal-area world map produced by `reproject`."""
+    fwd, xmax, ymax = {"equal_earth": (equal_earth_forward, EE_XMAX, EE_YMAX),
+                       "mollweide": (mollweide_forward, MW_XMAX, MW_YMAX)}[proj]
+    h, w = arr.shape[:2]
+    im = Image.fromarray(arr)
+    draw = ImageDraw.Draw(im)
+    to_px = lambda x, y: list(zip(((x / xmax + 1) / 2 * w).tolist(), ((1 - y / ymax) / 2 * h).tolist()))
+    t = np.linspace(-90, 90, 181)
+    for lon in range(-180, 181, step):
+        draw.line(to_px(*fwd(t, np.full_like(t, float(lon)))), fill=(200, 210, 225), width=1)
+    s = np.linspace(-180, 180, 361)
+    for lat in range(-90 + step, 90, step):
+        draw.line(to_px(*fwd(np.full_like(s, float(lat)), s)), fill=(200, 210, 225), width=1)
+    return np.array(im)
+
+
+def continent_centres(g, ocean, min_share=0.03):
+    """Centroids (lat, lon) of land bodies holding ≥ min_share of all land, largest first."""
+    from .graph import components
+    land = ~np.asarray(ocean)
+    lab = components(g, land)
+    area = np.bincount(lab[land], weights=g.area_km2[land])
+    out = []
+    for k in np.argsort(-area):
+        if area[k] < min_share * area.sum():
+            break
+        m = lab == k
+        c = (g.xyz[m] * g.area_km2[m, None]).sum(axis=0)
+        la, lo = xyz_to_latlon(c[None])
+        out.append((float(la[0]), float(lo[0])))
+    return out
+
+
+def write_all(relief, pdir, g, ocean, width, globe_size=1024):
+    """Write proj_equal_earth.png, proj_mollweide.png, globe_*.png and globes_sheet.png into pdir."""
+    for proj in ("equal_earth", "mollweide"):
+        Image.fromarray(graticule(reproject(relief, proj, width), proj)).save(pdir / f"proj_{proj}.png")
+    views = [("north_pole", 90.0, 0.0), ("south_pole", -90.0, 0.0)]
+    views = [(f"continent_{i + 1}", la, lo) for i, (la, lo) in enumerate(continent_centres(g, ocean))] + views
+    tiles = []
+    for name, la, lo in views:
+        im = Image.fromarray(globe(relief, la, lo, globe_size))
+        im.save(pdir / f"globe_{name}.png")
+        tiles.append((f"{name} ({la:.0f}°, {lo:.0f}°)", im))
+    cols, t = 4, globe_size // 2
+    rows = (len(tiles) + cols - 1) // cols
+    sheet = Image.new("RGB", (cols * t, rows * (t + 16)), BACKGROUND)
+    draw = ImageDraw.Draw(sheet)
+    for k, (label, im) in enumerate(tiles):
+        x, y = (k % cols) * t, (k // cols) * (t + 16)
+        sheet.paste(im.resize((t, t)), (x, y + 16))
+        draw.text((x + 4, y + 2), label, fill=(230, 230, 230))
+    sheet.save(pdir / "globes_sheet.png")
diff --git a/mapgen/render.py b/mapgen/render.py
new file mode 100644
index 0000000..2a408e6
--- /dev/null
+++ b/mapgen/render.py
@@ -0,0 +1,526 @@
+"""Stage `render`: equirectangular rasters, cells.npz, metadata, previews, contact sheet."""
+from __future__ import annotations
+
+import colorsys
+import json
+
+import numpy as np
+from PIL import Image, ImageDraw
+from scipy.spatial import cKDTree
+
+from . import geo, projections, viewer_export
+from . import plateaus as PL
+from .crust import AGE_NAMES
+from .environment import LEGENDS, REGIONS, ZONES
+from .ice import ICE_NAMES
+from .noise import fbm
+from .seabed import MINERAL_NAMES, SEABED_NAMES
+from .sphere import east_north, latlon_to_xyz
+
+CONTINUOUS = {
+    "elevation": ("z_surface_m", 0.5, -12000.0, "m"),
+    "T_mean": ("T_mean", 0.01, -100.0, "degC"),
+    "T_range": ("T_range", 0.01, 0.0, "degC"),
+    "P_ann": ("P_ann", 0.5, 0.0, "mm/yr"),
+    "P_jun": ("P_jun", 0.5, 0.0, "mm/yr"),
+    "P_dec": ("P_dec", 0.5, 0.0, "mm/yr"),
+    "po2": ("po2_bar", 2e-4, 0.0, "bar"),
+    "gravity": ("gravity_g", 1e-4, 0.0, "g"),
+    "pressure": ("pressure_bar", 2e-4, 0.0, "bar"),
+    "o2_fraction": ("o2_fraction", 2e-5, 0.0, "fraction"),
+    "fire": ("fire_reactivity", 1e-4, 0.0, "x"),
+    "vent_potential": ("vent_potential", 1e-4, 0.0, "0-1"),
+    "bottom_temp": ("bottom_temp_c", 0.01, -10.0, "degC"),
+    "sediment": ("sediment_m", 0.5, 0.0, "m"),
+    "plant_height": ("plant_height_x", 2e-3, 0.0, "x"),
+    "sst": ("sst", 0.01, -100.0, "degC"),
+    "productivity": ("productivity", 1e-4, 0.0, "0-1"),
+    "current_speed": ("current_speed", 1e-4, 0.0, "m/s"),
+    "upwelling": ("upwelling", 0.05, -1600.0, "m/yr"),
+}
+CATEGORICAL = {"plates": "plate", "age_class": "age_class", "holdridge": "holdridge",
+               "seasonality": "seasonality", "landform": "landform", "lithology": "lithology",
+               "ground": "ground", "ice": "ice",
+               "seabed_type": "seabed_type", "seabed_mineral": "seabed_mineral", "deposits": "deposit_main"}
+DETAIL_M = np.array([60, 40, 150, 400, 80, 150, 300, 300, 300, 60, 30, 120], dtype=np.float64)  # by landform
+CHUNK = 128
+COAST_DETAIL_M = 250.0
+
+
+def _by_rows(fn, v, dtype, tail=()):
+    """fn applied to CHUNK-row slices of v (fn elementwise): the same values as fn(v), without fn's full-size
+    float64 scratch (a raster ramp builds five H×W×3 float64 temporaries: ~4 GB at 8192 px)."""
+    v = np.asarray(v)
+    if v.ndim < 2 or v.shape[0] <= CHUNK:
+        return fn(v)
+    out = np.empty(v.shape + tuple(tail), dtype)
+    for y0 in range(0, v.shape[0], CHUNK):
+        out[y0:y0 + CHUNK] = fn(v[y0:y0 + CHUNK])
+    return out
+
+
+def encode(v, scale, offset):
+    return _by_rows(lambda x: np.clip(np.round((np.asarray(x, dtype=np.float64) - offset) / scale), 0, 65535)
+                    .astype(np.uint16), v, np.uint16)
+
+
+def decode(raw, scale, offset):
+    return np.asarray(raw, dtype=np.float64) * scale + offset
+
+
+def _rows_xyz(rows, W, H):
+    lat = 90.0 - (rows + 0.5) / H * 180.0
+    lon = (np.arange(W) + 0.5) / W * 360.0 - 180.0
+    LA, LO = np.meshgrid(lat, lon, indexing="ij")
+    return latlon_to_xyz(LA.ravel(), LO.ravel())
+
+
+def _make_render_jit():
+    """Compiled ramp and sampling loops (numba optional; WORLDGEN_NO_JIT=1 turns it off): per pixel the same float
+    steps in the same order as the numpy statements they replace, so the same bytes, without the temporaries."""
+    import os
+    if os.environ.get("WORLDGEN_NO_JIT"):
+        return None
+    try:
+        import numba
+    except ImportError:
+        return None
+
+    @numba.njit(cache=True, nogil=True)
+    def ramp(x, lo, span, stops, out):          # x: float64 (n,), out: uint8 (n, c)
+        m = stops.shape[0]
+        for j in range(x.shape[0]):
+            t = (x[j] - lo) / span
+            t = 0.0 if t < 0.0 else (1.0 if t > 1.0 else t)
+            t = t * (m - 1)
+            i = min(np.int64(t), m - 2)
+            f = t - i
+            for c in range(stops.shape[1]):
+                out[j, c] = np.uint8(np.int64(stops[i, c] * (1 - f) + stops[i + 1, c] * f))
+
+    @numba.njit(cache=True, nogil=True)
+    def sample(v, idx, wts, out):               # out[p] = Σ_k v[idx[p,k]] * w[p,k], k left to right (as np.sum, k < 8)
+        for p in range(idx.shape[0]):
+            s = v[idx[p, 0]] * np.float64(wts[p, 0])
+            for k in range(1, idx.shape[1]):
+                s += v[idx[p, k]] * np.float64(wts[p, k])
+            out[p] = s
+    return ramp, sample
+
+
+_render_jit = _make_render_jit()
+
+
+def pixel_neighbours(g, W, H, k=3):
+    from .graph import workers
+    tree = cKDTree(g.xyz)
+    idx = np.empty((H, W, k), np.int32)
+    wts = np.empty((H, W, k), np.float32)
+    for y0 in range(0, H, CHUNK):
+        rows = np.arange(y0, min(H, y0 + CHUNK))
+        d, i = tree.query(_rows_xyz(rows, W, H), k=k, workers=workers())   # per point: the same answer
+        w = 1.0 / np.maximum(d, 1e-9) ** 2
+        w /= w.sum(axis=1, keepdims=True)
+        idx[rows] = i.reshape(len(rows), W, k)
+        wts[rows] = w.reshape(len(rows), W, k)
+    return idx, wts
+
+
+def sample_cont(v, idx, wts):
+    out = np.empty(idx.shape[:2])
+    if (_render_jit is not None and v.dtype == np.float64 and v.ndim == 1 and wts.dtype == np.float32
+            and 1 <= idx.shape[-1] < 8 and idx.shape == wts.shape):
+        k = idx.shape[-1]
+        _render_jit[1](np.ascontiguousarray(v), np.ascontiguousarray(idx).reshape(-1, k),
+                       np.ascontiguousarray(wts).reshape(-1, k), out.reshape(-1))
+        return out
+    for y0 in range(0, idx.shape[0], CHUNK):
+        s = slice(y0, y0 + CHUNK)
+        out[s] = np.sum(v[idx[s]] * wts[s], axis=-1)
+    return out
+
+
+def sample_cat(v, idx):
+    return v[idx[..., 0]]
+
+
+def pixel_land(ocean_k, z):
+    """Pixel land mask: unanimous neighbour cells decide; at the coast (mixed) the sub-cell height does."""
+    all_sea = ocean_k.all(axis=-1)
+    all_land = ~ocean_k.any(axis=-1)
+    return all_land | (~all_sea & ~all_land & (z > 0))
+
+
+def hillshade(z, radius_km, az=315.0, alt=45.0, exag=4.0, low_memory=False):
+    """low_memory: the same values, CHUNK rows at a time (one halo row each side keeps the central differences)."""
+    if low_memory:
+        H = z.shape[0]
+        out = np.empty(z.shape)
+        for y0 in range(0, H, CHUNK):
+            a, b = max(y0 - 1, 0), min(y0 + CHUNK + 1, H)
+            part = _hillshade_rows(z[a:b], a, H, radius_km, az, alt, exag)
+            out[y0:min(y0 + CHUNK, H)] = part[y0 - a: y0 - a + min(CHUNK, H - y0)]
+        return out
+    return _hillshade_rows(z, 0, z.shape[0], radius_km, az, alt, exag)
+
+
+def _hillshade_rows(z, row0, H, radius_km, az, alt, exag):
+    """Hillshade of rows row0.. of an H-row raster (z holds those rows; edges of z use one-sided differences)."""
+    W = z.shape[1]
+    lat = 90.0 - (np.arange(row0, row0 + z.shape[0]) + 0.5) / H * 180.0
+    dy = np.pi * radius_km * 1000.0 / H
+    dx = 2 * np.pi * radius_km * 1000.0 * np.maximum(np.cos(np.radians(lat)), 0.01) / W
+    gy, gx = np.gradient(z)
+    dzdx = gx / dx[:, None] * exag
+    dzdn = -gy / dy * exag
+    norm = np.sqrt(dzdx**2 + dzdn**2 + 1.0)
+    a, b = np.radians(az), np.radians(alt)
+    L = (np.sin(a) * np.cos(b), np.cos(a) * np.cos(b), np.sin(b))
+    return np.clip((-dzdx * L[0] - dzdn * L[1] + L[2]) / norm, 0.0, 1.0)
+
+
+def _ramp(v, lo, hi, stops):
+    stops = np.asarray(stops, dtype=np.float64)
+
+    def part(x):
+        if (_render_jit is not None and isinstance(x, np.ndarray) and x.dtype == np.float64 and stops.ndim == 2
+                and len(stops) >= 2 and not np.isnan(x).any()):
+            out = np.empty(x.shape + (stops.shape[1],), np.uint8)
+            _render_jit[0](np.ascontiguousarray(x).reshape(-1), float(lo), float(hi - lo), stops,
+                           out.reshape(-1, stops.shape[1]))
+            return out
+        t = np.clip((x - lo) / (hi - lo), 0, 1) * (len(stops) - 1)
+        i = np.minimum(t.astype(np.int64), len(stops) - 2)
+        f = (t - i)[..., None]
+        return (stops[i] * (1 - f) + stops[i + 1] * f).astype(np.uint8)
+    return _by_rows(part, v, np.uint8, (stops.shape[-1],))
+
+
+def holdridge_palette():
+    ramp = np.array([[216, 200, 160], [208, 196, 140], [200, 200, 120], [152, 168, 96], [106, 150, 80],
+                     [70, 125, 68], [50, 105, 62], [37, 90, 56]], dtype=np.float64)
+    tint = {"polar": ((232, 236, 239), 0.9), "subpolar": ((170, 176, 160), 0.55), "boreal": ((70, 100, 80), 0.35),
+            "tropical": ((20, 90, 40), 0.15)}
+    cols = []
+    for r in REGIONS:
+        k = len(ZONES[r])
+        for j in range(k):
+            c = ramp[int(round((j / max(k - 1, 1)) * (len(ramp) - 1)))]
+            if r in tint:
+                c = c * (1 - tint[r][1]) + np.array(tint[r][0]) * tint[r][1]
+            cols.append(c)
+    return np.array(cols, dtype=np.uint8)
+
+
+def category_palette(n):
+    return np.array([[int(255 * c) for c in colorsys.hsv_to_rgb((i * 0.618034) % 1.0, 0.55, 0.9)]
+                     for i in range(max(n, 1))], dtype=np.uint8)
+
+
+def _region_palette(cols):
+    """38 zone colours from per-region (dry, wet) colour pairs."""
+    out = []
+    for r in REGIONS:
+        dry, wet = (np.array(c, dtype=np.float64) for c in cols[r])
+        k = len(ZONES[r])
+        out += [dry + (wet - dry) * (j / max(k - 1, 1)) for j in range(k)]
+    return np.array(out, dtype=np.uint8)
+
+
+STYLES = {  # alien palettes (config [render] style); None = Earth-like
+    "tidal-lock": {
+        "zones": {"polar": ((34, 38, 46), (34, 38, 46)), "subpolar": ((70, 72, 78), (96, 104, 112)),
+                  "boreal": ((120, 96, 64), (150, 110, 60)), "cool temperate": ((170, 120, 60), (190, 140, 70)),
+                  "warm temperate": ((110, 70, 44), (130, 86, 50)), "subtropical": ((46, 36, 34), (60, 44, 38)),
+                  "tropical": ((22, 20, 22), (30, 26, 28))},
+        "ocean": [[4, 12, 16], [10, 34, 40], [40, 80, 84]], "lake": (60, 96, 104), "river": (70, 110, 118),
+        "ground": {6: (92, 86, 80)}, "sheet": (200, 222, 240), "sea_ice": (170, 200, 226), "sea_ice_seasonal": (120, 150, 170)},
+    "salt-mirror": {
+        "zones": {"polar": ((236, 240, 244), (236, 240, 244)), "subpolar": ((228, 228, 234), (214, 214, 228)),
+                  "boreal": ((238, 236, 228), (208, 208, 224)), "cool temperate": ((242, 238, 226), (200, 204, 222)),
+                  "warm temperate": ((240, 230, 204), (196, 204, 220)), "subtropical": ((238, 224, 186), (192, 206, 218)),
+                  "tropical": ((234, 214, 166), (186, 206, 216))},
+        "ocean": [[70, 120, 130], [150, 195, 200], [215, 235, 235]], "lake": (150, 196, 204), "river": (150, 196, 204),
+        "ground": {6: (252, 250, 244)}, "sheet": (248, 250, 252), "sea_ice": (236, 244, 246), "sea_ice_seasonal": (220, 236, 238)},
+}
+
+
+def style_palette(style):
+    return None if style is None else STYLES[style]
+
+
+def relief_rgb(z, hs, zone, ground, ice, lake, land=None, vary=None, style=None, low_memory=False):
+    """vary: optional (brightness factor, tint) per pixel for open land (not lakes or ice): tint > 0 drier/yellower,
+    < 0 lusher (deeper green). style: a STYLES key (alien palette) or None. low_memory: the same pixels, made
+    CHUNK rows at a time (no full-size float64 scratch)."""
+    if low_memory:
+        out = np.empty(np.shape(z) + (3,), np.uint8)
+        for y0 in range(0, np.shape(z)[0], CHUNK):
+            s = slice(y0, y0 + CHUNK)
+            out[s] = relief_rgb(z[s], hs[s], zone[s], ground[s], ice[s], lake[s], None if land is None else land[s],
+                                None if vary is None else tuple(np.asarray(v)[s] for v in vary), style)
+        return out
+    land = z > 0 if land is None else land
+    st = style_palette(style)
+    pal = holdridge_palette() if st is None else _region_palette(st["zones"])
+    rgb = pal[np.clip(zone, 0, 37)].astype(np.float64)
+    ocean = _ramp(z, -6500.0 if st is None else -1500.0, 0.0,
+                  [[11, 43, 90], [30, 90, 150], [143, 198, 224]] if st is None else st["ocean"]).astype(np.float64)
+    rgb = np.where(land[..., None], rgb, ocean)
+    grounds = ((1, (111, 143, 106)), (2, (120, 128, 100)), (5, (63, 111, 74)), (6, (239, 233, 220))) if st is None \
+        else tuple(st["ground"].items())
+    for code, col in grounds:
+        rgb = np.where((land & (ground == code))[..., None], col, rgb)
+    if vary is not None:
+        bright, tint = (np.asarray(v, dtype=np.float64)[..., None] for v in vary)
+        shift = np.where(tint > 0, tint * np.array([1.0, 0.6, -0.8]), -tint * np.array([-0.9, -0.2, -0.6]))
+        if st is not None:
+            shift = np.abs(tint) * np.array([0.4, 0.4, 0.4]) * np.sign(tint)
+        rgb = np.where(land[..., None], rgb * bright + shift, rgb)
+    rgb = np.where((lake & land)[..., None], (79, 143, 192) if st is None else st["lake"], rgb)
+    rgb = np.where(np.isin(ice, [1, 2])[..., None], (244, 248, 251) if st is None else st["sheet"], rgb)
+    rgb = np.where((ice == 4)[..., None], (225, 235, 242) if st is None else st["sea_ice"], rgb)
+    rgb = np.where((ice == 3)[..., None], 0.5 * rgb + 0.5 * np.array([220, 232, 240] if st is None else st["sea_ice_seasonal"]), rgb)
+    shade = np.where(land, 0.55 + 0.45 * hs, 0.85 + 0.15 * hs)
+    return np.clip(rgb * shade[..., None], 0, 255).astype(np.uint8)
+
+
+RIVER_RGB = (58, 112, 176)
+
+
+def draw_rivers(rgb, g, recv, river, strahler, min_order=2, colour=RIVER_RGB):
+    """Draw river segments (cell centre → receiver) of order ≥ min_order; width grows with order."""
+    H, W = rgb.shape[:2]
+    im = Image.fromarray(rgb)
+    draw = ImageDraw.Draw(im)
+    x = (np.asarray(g.lon) + 180.0) / 360.0 * W - 0.5
+    y = (90.0 - np.asarray(g.lat)) / 180.0 * H - 0.5
+    cells = np.flatnonzero(river & (recv != np.arange(g.n)) & (strahler >= min_order))
+    for i in cells[np.argsort(strahler[cells])]:
+        j = recv[i]
+        if abs(x[i] - x[j]) > W / 2:
+            continue
+        width = max(1, int(round(int(strahler[i]) * W / 8192)))
+        draw.line([(x[i], y[i]), (x[j], y[j])], fill=tuple(colour), width=width)
+    return np.array(im)
+
+
+def draw_currents(rgb, g, current, ocean, per_row=60, colour=(255, 255, 255)):
+    """Arrows along the surface current on a lattice ≈ `per_row` across the map; length ∝ speed (1 m/s ≈ one
+    lattice step), sea only; currents under 2 cm/s get none."""
+    H, W = rgb.shape[:2]
+    step = W / per_row
+    current = np.asarray(current, dtype=np.float64)
+    ocean = np.asarray(ocean, bool)
+    e, n = east_north(g.xyz)
+    tree = cKDTree(g.xyz)
+    ys, xs = np.meshgrid(np.arange(step / 2, H, step), np.arange(step / 2, W, step), indexing="ij")
+    lat, lon = 90.0 - (ys.ravel() + 0.5) / H * 180.0, (xs.ravel() + 0.5) / W * 360.0 - 180.0
+    _, cell = tree.query(latlon_to_xyz(lat, lon))
+    ue, un = np.sum(current[cell] * e[cell], axis=1), np.sum(current[cell] * n[cell], axis=1)
+    sp = np.hypot(ue, un)
+    im = Image.fromarray(rgb)
+    draw = ImageDraw.Draw(im)
+    for x, y, a, b, s, ok in zip(xs.ravel(), ys.ravel(), ue, un, sp, ocean[cell]):
+        if not ok or s < 0.02:
+            continue
+        L = min(s, 1.5) * 0.9 * step
+        dx, dy = a / s * L, -b / s * L
+        x0, y0, x1, y1 = x - dx / 2, y - dy / 2, x + dx / 2, y + dy / 2
+        draw.line([(x0, y0), (x1, y1)], fill=tuple(colour), width=1)
+        for ang in (2.6, -2.6):                                        # arrowhead: two barbs ±150°
+            c, sn = np.cos(ang), np.sin(ang)
+            draw.line([(x1, y1), (x1 + 0.35 * (c * dx - sn * dy), y1 + 0.35 * (sn * dx + c * dy))],
+                      fill=tuple(colour), width=1)
+    return np.array(im)
+
+
+RAMPS = {  # key: (unit, lo, hi, colour stops) — shared by previews, viewer textures and legends
+    "elevation": ("m", -6000.0, 6000.0, [[8, 30, 70], [30, 90, 150], [140, 200, 225], [60, 120, 60],
+                                          [150, 160, 90], [140, 110, 80], [235, 235, 235]]),
+    "T_mean": ("degC", -40.0, 40.0, [[40, 60, 160], [240, 240, 240], [180, 30, 30]]),
+    "T_range": ("degC", 0.0, 60.0, [[68, 1, 84], [33, 145, 140], [253, 231, 37]]),
+    "P": ("mm/yr", 0.0, 4000.0, [[150, 110, 60], [230, 220, 150], [60, 150, 70], [30, 70, 170]]),
+    "po2": ("bar", 0.08, 0.32, [[60, 40, 90], [240, 240, 240], [200, 90, 20]]),
+    "gravity": ("g", 0.3, 1.8, [[20, 120, 180], [240, 240, 240], [120, 40, 40]]),
+    "pressure": ("bar", 0.3, 2.5, [[40, 60, 120], [240, 240, 240], [150, 60, 30]]),
+    "o2_fraction": ("fraction", 0.10, 0.40, [[60, 40, 90], [240, 240, 240], [200, 90, 20]]),
+    "fire": ("x", 0.0, 1.5, [[40, 80, 160], [240, 240, 240], [220, 120, 20]]),
+    "vent_potential": ("0-1", 0.0, 1.0, [[10, 20, 50], [60, 90, 160], [250, 190, 60], [255, 250, 220]]),
+    "bottom_temp": ("degC", -2.0, 30.0, [[30, 40, 120], [80, 160, 200], [240, 200, 120]]),
+    "sediment": ("m", 0.0, 3000.0, [[40, 30, 60], [140, 110, 80], [240, 220, 170]]),
+    "plant_height": ("x", 0.5, 3.5, [[150, 120, 60], [240, 240, 240], [30, 110, 50]]),
+    "sst": ("degC", -2.0, 32.0, [[30, 40, 120], [60, 150, 200], [240, 240, 200], [220, 90, 40]]),
+    "productivity": ("0-1", 0.0, 1.0, [[10, 20, 60], [20, 110, 120], [120, 200, 90], [240, 240, 120]]),
+    "current_speed": ("m/s", 0.0, 1.0, [[10, 20, 50], [40, 90, 170], [120, 200, 230], [255, 255, 255]]),
+    "upwelling": ("m/yr", -200.0, 200.0, [[40, 60, 160], [240, 240, 240], [30, 140, 80]]),
+}
+VIEWER_LAYERS = [  # id, name, source (continuous raster name | categorical name | "relief")
+    ("relief", "Relief", "relief"), ("biomes", "Biomes (Holdridge)", "holdridge"),
+    ("elevation", "Elevation", "elevation"), ("temperature", "Mean temperature", "T_mean"),
+    ("rainfall", "Rainfall", "P_ann"), ("seasonality", "Seasonality", "seasonality"),
+    ("landform", "Landform", "landform"), ("ground", "Ground", "ground"), ("ice", "Ice", "ice"), ("deposits", "Mineral deposits", "deposits"),
+    ("plates", "Plates", "plates"), ("o2", "O₂ partial pressure", "po2"), ("gravity", "Gravity", "gravity"),
+    ("pressure", "Air pressure", "pressure"), ("fire", "Fire reactivity", "fire"),
+    ("seabed", "Sea-floor type", "seabed_type"), ("minerals", "Sea-floor minerals", "seabed_mineral"),
+    ("bottom_temp", "Bottom temperature", "bottom_temp"), ("sediment", "Sediment", "sediment"),
+    ("vent_potential", "Vent potential", "vent_potential"),
+    ("currents", "Ocean currents", "current_speed"), ("sst", "Sea-surface temperature", "sst"),
+    ("productivity", "Sea productivity", "productivity"),
+]
+
+
+def _ramp_key(name):
+    return "P" if name.startswith("P_") else name
+
+
+def _colorize(name, v):
+    _, lo, hi, stops = RAMPS[_ramp_key(name)]
+    return _ramp(v, lo, hi, stops)
+
+
+def _preview(img: Image.Image, width: int, nearest: bool) -> Image.Image:
+    return img.resize((width, width // 2), Image.NEAREST if nearest else Image.LANCZOS)
+
+
+def contact_sheet(tiles, tile_w):
+    cols = 3
+    th = tile_w // 2 + 14
+    rows = (len(tiles) + cols - 1) // cols
+    sheet = Image.new("RGB", (cols * tile_w, rows * th), (20, 20, 24))
+    draw = ImageDraw.Draw(sheet)
+    for k, (name, im) in enumerate(tiles):
+        x, y = (k % cols) * tile_w, (k // cols) * th
+        sheet.paste(im.resize((tile_w, tile_w // 2)), (x, y + 14))
+        draw.text((x + 4, y + 1), name, fill=(235, 235, 235))
+    return sheet
+
+
+class _Writer:
+    """Saves files on a background thread (PNG/zlib encoding releases the GIL) while the next layer is computed;
+    at most `depth` waiting, so memory stays bounded. The files are the same bytes as saved in line."""
+    def __init__(self, threads: int = 2, depth: int = 4):
+        from concurrent.futures import ThreadPoolExecutor
+        self.ex, self.pending, self.depth = ThreadPoolExecutor(threads), [], depth
+
+    def __call__(self, fn, *args, **kw):
+        while len(self.pending) >= self.depth:
+            self.pending.pop(0).result()
+        self.pending.append(self.ex.submit(fn, *args, **kw))
+
+    def close(self):
+        try:
+            for f in self.pending:
+                f.result()
+        finally:
+            self.ex.shutdown(wait=True, cancel_futures=True)
+
+
+def run(ctx) -> dict:
+    writer = _Writer()
+    try:
+        return _run(ctx, writer)
+    finally:
+        writer.close()
+
+
+def _run(ctx, save) -> dict:
+    g, cfg, d = ctx.grid, ctx.cfg, ctx.data
+    W = int(cfg["build"]["raster_width"])
+    if ctx.res < int(cfg["build"]["res_final"]):
+        W = int(cfg["build"].get("dev_raster_width", W))
+    H = W // 2
+    pw = int(cfg["build"]["preview_width"])
+    out = ctx.out_dir or ctx.root / "out" / f"r{ctx.res}"
+    rdir = out / "raster"
+    pdir = ctx.preview_dir or ctx.root / "previews" / f"r{ctx.res}"
+    rdir.mkdir(parents=True, exist_ok=True)
+    pdir.mkdir(parents=True, exist_ok=True)
+    idx, wts = pixel_neighbours(g, W, H)
+    legends = {**LEGENDS, "age_class": AGE_NAMES, "ice": ICE_NAMES, "plates": [p["id"] for p in ctx.tect["plate"]],
+               "seabed_type": SEABED_NAMES, "seabed_mineral": MINERAL_NAMES}
+    meta = {"width": W, "height": H, "projection": "equirectangular", "radius_km": g.radius_km,
+            "units": cfg.get("units", {}), "continuous": {}, "categorical": {}}
+    tiles = []
+    cats = {name: sample_cat(np.asarray(d[key]), idx) for name, key in CATEGORICAL.items()}
+
+    z = sample_cont(np.asarray(d["z_surface_m"], dtype=np.float64), idx, wts)
+    amp = DETAIL_M[np.clip(cats["landform"], 0, len(DETAIL_M) - 1)]
+    for y0 in range(0, H, CHUNK):
+        rows = np.arange(y0, min(H, y0 + CHUNK))
+        z[rows] += amp[rows] * fbm(_rows_xyz(rows, W, H), ctx.seed + 61, 5, 64.0).reshape(len(rows), W)
+    lake = sample_cat(np.asarray(d["lake"]), idx)
+    ocean_cells = np.asarray(d["ocean"]) if "ocean" in d else np.asarray(d["z_surface_m"]) <= 0
+    ocean_k = ocean_cells[idx]
+    coast = ocean_k.any(axis=-1) & ~ocean_k.all(axis=-1)
+    for y0 in range(0, H, CHUNK):                      # extra fine detail where the coastline runs
+        rows = np.arange(y0, min(H, y0 + CHUNK))
+        cz = COAST_DETAIL_M * fbm(_rows_xyz(rows, W, H), ctx.seed + 67, 5, 256.0).reshape(len(rows), W)
+        z[rows] += np.where(coast[rows], cz, 0.0)
+    land = pixel_land(ocean_k, z)
+    if "lake_level_m" in d:   # the water surface for drawing: lake cells at their level (beds stay in `elevation`)
+        lev = np.asarray(d["lake_level_m"], dtype=np.float64)
+        dz = sample_cont(np.where(np.isfinite(lev), lev - np.asarray(d["z_surface_m"], dtype=np.float64), 0.0), idx, wts)
+        save(Image.fromarray(encode(z + dz, *CONTINUOUS["elevation"][1:3])).save, rdir / "surface.png")
+        meta["continuous"]["surface"] = {"file": "surface.png", "scale": CONTINUOUS["elevation"][1],
+                                         "offset": CONTINUOUS["elevation"][2], "unit": "m"}
+    hs = hillshade(z, g.radius_km, low_memory=ctx.low_memory)
+    style = cfg.get("render", {}).get("style")
+    rgb = relief_rgb(z, hs, cats["holdridge"], cats["ground"], cats["ice"], lake, land, style=style,
+                     low_memory=ctx.low_memory)
+    rgb = draw_rivers(rgb, g, np.asarray(d["recv"]), np.asarray(d["river"]), np.asarray(d["strahler"]),
+                      colour=RIVER_RGB if style is None else STYLES[style]["river"])
+    relief = Image.fromarray(rgb)
+    projections.write_all(rgb, pdir, g, ocean_cells, pw, globe_size=max(pw // 2, 64))
+    save(relief.copy().save, rdir / "relief.png")
+    tiles.append(("relief", _preview(relief, pw, False)))
+
+    vdir = out / "viewer"
+    by_src = {}                                        # viewer layers are written as soon as their colours exist
+    for vid, vname, src in VIEWER_LAYERS:
+        by_src.setdefault(src, []).append((vid, vname))
+    ventries = {}
+
+    def emit(src, rgb_, legend):
+        for vid, vname in by_src.get(src, []):
+            ventries[vid] = viewer_export.write_layer(vdir, {"id": vid, "name": vname, "rgb": rgb_, "legend": legend})
+
+    emit("relief", rgb, None)
+    for name, (key, scale, offset, unit) in CONTINUOUS.items():
+        v = z if name == "elevation" else sample_cont(np.asarray(d[key], dtype=np.float64), idx, wts)
+        save(Image.fromarray(encode(v, scale, offset)).save, rdir / f"{name}.png")
+        meta["continuous"][name] = {"file": f"{name}.png", "scale": scale, "offset": offset, "unit": unit}
+        crgb = _colorize(name, v)
+        if name == "current_speed" and "current" in d:
+            crgb = draw_currents(crgb, g, d["current"], ocean_cells)
+        if name != "elevation":
+            tiles.append((name, _preview(Image.fromarray(crgb), pw, False)))
+        if name in by_src:
+            unit_, lo, hi, stops = RAMPS[_ramp_key(name)]
+            emit(name, crgb, viewer_export.continuous_legend(unit_, lo, hi, stops))
+        del crgb
+    for name, v in cats.items():
+        save(Image.fromarray(v.astype(np.uint8), "L").save, rdir / f"{name}.png")
+        meta["categorical"][name] = {"file": f"{name}.png", "legend": legends[name]}
+        pal = holdridge_palette() if name == "holdridge" else category_palette(len(legends[name]))
+        crgb = pal[np.clip(v, 0, len(pal) - 1)]
+        tiles.append((name, _preview(Image.fromarray(crgb), pw, True)))
+        emit(name, crgb, viewer_export.categorical_legend(legends[name], pal))
+        del crgb
+    (out / "fields.json").write_text(json.dumps(meta, indent=1))
+
+    per_cell = {k: v for k, v in d.items() if isinstance(v, np.ndarray) and v.shape[:1] == (g.n,)}
+    save(lambda: np.savez_compressed(out / "cells.npz", **per_cell))
+    (out / "cells_meta.json").write_text(json.dumps(
+        {"res": ctx.res, "radius_km": g.radius_km, "n_cells": g.n, "units": cfg.get("units", {}),
+         "planet": {**cfg["planet"], **({"sun_lock": cfg["climate"].get("lock_at", [0.0, 0.0])}
+                                        if cfg.get("climate", {}).get("lock") else {})}, "name": cfg.get("render", {}).get("name", "World"), "style": style,
+         "legends": legends, "fields": sorted(per_cell),
+         "plateaus": [PL.features(p, ctx.seed, g.radius_km) for p in ctx.tect.get("plateau", [])]}, indent=1))
+
+    for name, im in tiles:
+        im.save(pdir / f"{name}.png")
+    contact_sheet(tiles, max(pw // 3, 64)).save(pdir / "contact_sheet.png")
+    viewer_export.write_index(vdir, [ventries[vid] for vid, _, _ in VIEWER_LAYERS])
+    geo.write_all(ctx, out, land, lake)
+    return {}
diff --git a/mapgen/seabed.py b/mapgen/seabed.py
new file mode 100644
index 0000000..7715c44
--- /dev/null
+++ b/mapgen/seabed.py
@@ -0,0 +1,95 @@
+"""Stage `seabed`: sea-floor data for every ocean cell — vent potential,
+sea-floor type, minerals, bottom temperature, sediment. Land cells get 0 / "—" (the fields describe the sea floor)."""
+from __future__ import annotations
+
+import numpy as np
+from scipy.spatial import cKDTree
+
+from . import plateaus as PL
+from .config import params
+from .elevation import hotspot_track
+from .graph import distance_to
+from .noise import fbm
+from .sphere import latlon_to_xyz
+
+(SB_NONE, SB_TERRIGENOUS, SB_CARBONATE, SB_CLAY, SB_BASALT, SB_VOLCANIC, SB_CONTINENTAL, SB_VENTS,
+ SB_TRENCH) = range(9)
+SEABED_NAMES = ["—", "terrigenous sediment", "carbonate ooze", "pelagic clay", "young basalt", "volcanic",
+                "continental / plateau", "vent field", "trench"]
+MI_NONE, MI_SULFIDES, MI_NODULES, MI_COBALT, MI_PHOSPHORITE = range(5)
+MINERAL_NAMES = ["none", "polymetallic sulfides", "manganese nodules", "cobalt crusts", "phosphorite"]
+DEFAULTS = {"ridge_km": 100.0, "back_arc_km": 300.0, "back_arc_width_km": 150.0, "back_arc": 0.6,
+            "hotspot_km": 150.0, "hotspot_spacing_km": 150.0, "terrigenous_km": 300.0, "carbonate_depth_m": 4500.0,
+            "carbonate_t_c": 10.0, "young_myr": 10.0, "trench_km": 80.0, "trench_depth_m": 5000.0,
+            "vent_field": 0.8, "sulfides": 0.7, "volcanic": 0.3, "thermocline_m": 800.0}
+
+
+def _cells_near(tree, radius_km, p, reach_km):
+    return np.asarray(tree.query_ball_point(p, 2.0 * np.sin(min(np.pi, reach_km / radius_km) / 2.0)), dtype=np.int64)
+
+
+def _bump(v, g, tree, p, r_km, amp):
+    i = _cells_near(tree, g.radius_km, p, 3.0 * r_km)
+    if len(i):
+        d = g.radius_km * np.arccos(np.clip(g.xyz[i] @ p, -1.0, 1.0))
+        v[i] = np.maximum(v[i], amp * np.exp(-(d / r_km) ** 2))
+
+
+def volcanic_potential(g, tect: dict, vel, seed: int, P: dict):
+    """0–1 per cell from hotspot chains (strongest at the active end) and plateau volcanic fields."""
+    tree, v = cKDTree(g.xyz), np.zeros(g.n)
+    for h in tect.get("hotspot", []):
+        for p, s in hotspot_track(g, h, vel, P["hotspot_spacing_km"]):
+            _bump(v, g, tree, p, P["hotspot_km"], 1.0 - s / h["length_km"])
+    for pl in tect.get("plateau", []):
+        f, vent = PL.features(pl, seed, g.radius_km), float(pl.get("vent", 1.0))
+        _bump(v, g, tree, latlon_to_xyz(*f["hotspot"]), 200.0, vent)
+        for c in f["cones"]:
+            _bump(v, g, tree, latlon_to_xyz(c["lat"], c["lon"]), c["radius_km"] + 30.0, 0.8 * vent)
+        for c in f["calderas"]:
+            _bump(v, g, tree, latlon_to_xyz(c["lat"], c["lon"]), c["radius_km"] + 30.0, vent)
+    return v
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "seabed", DEFAULTS)
+    z, ocean, cont, age, d_div, d_over, d_sub, t_mean, vel = ctx.need(
+        "elevation_eroded_m", "ocean", "continental", "ocean_age_myr", "d_div_km", "d_over_km", "d_sub_km", "T_mean",
+        "vel")
+    z, ocean, cont = np.asarray(z, dtype=np.float64), np.asarray(ocean, bool), np.asarray(cont, bool)
+    age, t_mean = np.asarray(age, dtype=np.float64), np.asarray(t_mean, dtype=np.float64)
+    plateau_id = np.asarray(ctx.data.get("plateau_id", np.full(g.n, -1)))
+    depth = np.maximum(-z, 0.0)
+    volc = volcanic_potential(g, ctx.tect, np.asarray(vel), ctx.seed, P)
+    ridge = np.exp(-(np.asarray(d_div) / P["ridge_km"]) ** 2)
+    back_arc = P["back_arc"] * np.exp(-((np.asarray(d_over) - P["back_arc_km"]) / P["back_arc_width_km"]) ** 2)
+    cluster = np.clip(0.5 + 1.5 * fbm(g.xyz, ctx.seed + 301, 4, 30.0), 0.0, 1.0)    # vents come in fields
+    vent = np.where(ocean, np.maximum.reduce([ridge, back_arc, volc]) * (0.5 + 0.5 * cluster), 0.0)
+    d_land = distance_to(g, ~ocean) if (~ocean).any() else np.full(g.n, np.inf)
+
+    t = np.full(g.n, SB_CLAY, np.int8)
+    t[(depth < P["carbonate_depth_m"]) & (t_mean > P["carbonate_t_c"])] = SB_CARBONATE
+    t[cont] = SB_CONTINENTAL
+    t[d_land < P["terrigenous_km"]] = SB_TERRIGENOUS
+    t[~cont & (age < P["young_myr"])] = SB_BASALT
+    t[volc > P["volcanic"]] = SB_VOLCANIC
+    t[(np.asarray(d_sub) < P["trench_km"]) & (depth > P["trench_depth_m"])] = SB_TRENCH
+    t[vent > P["vent_field"]] = SB_VENTS
+    t[~ocean] = SB_NONE
+
+    m = np.zeros(g.n, np.int8)
+    m[(t == SB_CLAY) & (age > 30.0) & (d_land > 1000.0) & (depth > 4000.0)] = MI_NODULES
+    m[((t == SB_VOLCANIC) | (plateau_id >= 0)) & (depth >= 800.0) & (depth <= 2500.0)] = MI_COBALT
+    m[cont & (depth < 500.0)] = MI_PHOSPHORITE
+    m[vent > P["sulfides"]] = MI_SULFIDES
+    m[~ocean] = MI_NONE
+
+    deep = 1.0 + 3.0 * np.cos(np.radians(g.lat)) ** 2                          # ≈ 1 °C polar … 4 °C tropical
+    w = np.clip((P["thermocline_m"] - depth) / P["thermocline_m"], 0.0, 1.0)   # shallow floors: towards the surface
+    bottom = np.where(ocean, deep + w * (np.maximum(t_mean, -1.8) - deep), 0.0)
+    terr = 2500.0 * np.exp(-d_land / 250.0)                                     # aprons off the continents
+    sed = np.where(cont, 300.0 + terr, np.minimum(5.0 * age, 800.0) + terr)     # pelagic rain ≈ 5 m per Myr
+    sed = np.where(ocean, sed, 0.0)
+    return {"vent_potential": vent.astype(np.float32), "seabed_type": t, "seabed_mineral": m,
+            "bottom_temp_c": bottom.astype(np.float32), "sediment_m": sed.astype(np.float32)}
diff --git a/mapgen/sketch.py b/mapgen/sketch.py
new file mode 100644
index 0000000..a0249bf
--- /dev/null
+++ b/mapgen/sketch.py
@@ -0,0 +1,121 @@
+"""Stage `sketch`: sample the sketch + user masks onto cells."""
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+
+from . import zones as ZN
+from .config import params
+from .noise import fbm
+from .pipeline import StageError
+from .sphere import latlon_to_xyz, xyz_to_latlon
+
+W, H = 2000, 1000
+SKETCH = ("land", "mountains", "desert", "rainforest", "trench")
+MASKS = ("land_hint", "mountain_hint", "o2_zones", "gravity_zones", "lock")
+DEFAULTS = {"warp_km": 1000.0, "warp_freq": 1.2, "detail_warp_km": 200.0, "detail_warp_freq": 5.0, "moves": []}
+
+
+def load_png01(path: Path) -> np.ndarray:
+    return np.array(Image.open(path).convert("L"), dtype=np.float64) / 255.0
+
+
+def load_mask(path: Path) -> np.ndarray:
+    """Grayscale override mask → [−1, 1]; mid-grey (128 / 32768), transparency and missing detail = neutral 0."""
+    try:
+        im = Image.open(path)
+        im.load()
+    except (OSError, ValueError) as e:
+        raise StageError(f"cannot read mask {path.name}: {e}") from e
+    if im.mode in ("I", "I;16", "I;16B", "I;16L"):
+        v = (np.array(im).astype(np.float64) - 32768.0) / 32768.0
+        dead = 1.0 / 32768.0
+    else:
+        la = np.array(im.convert("LA")).astype(np.float64)
+        v = np.where(la[..., 1] > 0, (la[..., 0] - 127.5) / 127.5, 0.0)
+        dead = 1.0 / 255.0
+    return np.clip(np.where(np.abs(v) <= dead, 0.0, v), -1.0, 1.0)
+
+
+def sample_equirect(img: np.ndarray, lat, lon) -> np.ndarray:
+    h, w = img.shape
+    x = (np.asarray(lon, dtype=np.float64) + 180.0) / 360.0 * w - 0.5
+    y = (90.0 - np.asarray(lat, dtype=np.float64)) / 180.0 * h - 0.5
+    x0 = np.floor(x).astype(np.int64)
+    y0 = np.floor(y).astype(np.int64)
+    fx, fy = x - x0, y - y0
+    xa, xb = x0 % w, (x0 + 1) % w
+    ya, yb = np.clip(y0, 0, h - 1), np.clip(y0 + 1, 0, h - 1)
+    top = img[ya, xa] * (1 - fx) + img[ya, xb] * fx
+    bot = img[yb, xa] * (1 - fx) + img[yb, xb] * fx
+    return top * (1 - fy) + bot * fy
+
+
+def rotation(a, b):
+    """3×3 rotation carrying [lat, lon] a onto b along the great circle (Rodrigues)."""
+    pa, pb = latlon_to_xyz(*np.array(a, dtype=np.float64)), latlon_to_xyz(*np.array(b, dtype=np.float64))
+    k = np.cross(pa, pb)
+    s, c = np.linalg.norm(k), float(np.dot(pa, pb))
+    if s < 1e-12:
+        return np.eye(3)
+    k /= s
+    K = np.array([[0, -k[2], k[1]], [k[2], 0, -k[0]], [-k[1], k[0], 0]])
+    return np.eye(3) + s * K + (1 - c) * K @ K
+
+
+def continent_mask(land, lat, lon):
+    """Pixels of the sketch landmass containing (lat, lon), 2 px coastal fringe included; lon wraps."""
+    from scipy import ndimage
+    lab, _ = ndimage.label(land)
+    for a, b in zip(lab[:, 0], lab[:, -1]):
+        if a and b and a != b:
+            lab[lab == b] = a
+    h, w = land.shape
+    y, x = min(max(int((90.0 - lat) / 180.0 * h), 0), h - 1), int((lon + 180.0) / 360.0 * w) % w
+    if not lab[y, x]:
+        raise StageError(f"sketch move: no sketch land at [{lat}, {lon}]")
+    return ndimage.binary_dilation(lab == lab[y, x], iterations=2)
+
+
+def warped_latlon(g, P, seed):
+    """Sample positions displaced by a two-scale domain warp (RMS ≈ warp_km, detail_warp_km)."""
+    p = g.xyz.copy()
+    for amp, freq, s in ((P["warp_km"], P["warp_freq"], 201), (P["detail_warp_km"], P["detail_warp_freq"], 211)):
+        if amp <= 0:
+            continue
+        disp = np.stack([fbm(g.xyz, seed + s + k, 4, freq) for k in range(3)], axis=1)
+        p = p + (amp / g.radius_km) * disp / max(float(disp.std()), 1e-12)
+    return xyz_to_latlon(p)
+
+
+def run(ctx) -> dict:
+    g = ctx.grid
+    P = params(ctx.cfg, "sketch", DEFAULTS)
+    sk = ctx.root / "sketch"
+    missing = [n for n in SKETCH if not (sk / f"{n}.png").exists()]
+    if missing:
+        raise StageError(f"sketch files missing {missing}: draw them or run `mapgen.py new-world` first")
+    lat, lon = warped_latlon(g, P, ctx.seed)     # the drawing is a loose guide: warp it (user masks are not)
+    imgs = {n: load_png01(sk / f"{n}.png") for n in SKETCH}
+    moved = []                                   # [[sketch.moves]]: rotate whole landmasses across the sphere
+    for mv in P["moves"]:
+        comp = continent_mask(imgs["land"] > 0.5, *mv["at"])
+        moved.append((rotation(mv["at"], mv["to"]), {n: np.where(comp, im, 0.0) for n, im in imgs.items()}))
+        imgs = {n: np.where(comp, 0.0, im) for n, im in imgs.items()}
+    p = latlon_to_xyz(lat, lon)
+    out = {}
+    for n in SKETCH:
+        v = sample_equirect(imgs[n], lat, lon)
+        for rot, parts in moved:
+            qlat, qlon = xyz_to_latlon(p @ rot)          # rot.T applied to row vectors
+            v = np.maximum(v, sample_equirect(parts[n], qlat, qlon))
+        out[f"sk_{n}"] = v
+    for m in MASKS:
+        p = ctx.root / "masks" / f"{m}.png"
+        out[f"m_{m}"] = sample_equirect(load_mask(p), g.lat, g.lon) if p.exists() else np.zeros(g.n)
+    zo2, zg = ZN.masks(g.xyz, ctx.tect.get("zone", []), ctx.seed, g.radius_km)   # config zones on top of the masks
+    out["m_o2_zones"] = np.clip(out["m_o2_zones"] + zo2, -1.0, 1.0)
+    out["m_gravity_zones"] = np.clip(out["m_gravity_zones"] + zg, -1.0, 1.0)
+    return out
diff --git a/mapgen/sphere.py b/mapgen/sphere.py
new file mode 100644
index 0000000..394b9ef
--- /dev/null
+++ b/mapgen/sphere.py
@@ -0,0 +1,68 @@
+"""Unit-sphere geometry and plate kinematics. Vectors are (..., 3); lat/lon in degrees."""
+from __future__ import annotations
+
+import numpy as np
+
+
+def latlon_to_xyz(lat, lon):
+    la, lo = np.radians(lat), np.radians(lon)
+    return np.stack([np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)], axis=-1)
+
+
+def xyz_to_latlon(p):
+    p = p / np.linalg.norm(p, axis=-1, keepdims=True)
+    return np.degrees(np.arcsin(np.clip(p[..., 2], -1, 1))), np.degrees(np.arctan2(p[..., 1], p[..., 0]))
+
+
+def east_north(p):
+    """Local unit east/north vectors. At the poles east is taken as +y (finite, arbitrary)."""
+    p = np.asarray(p, dtype=np.float64)
+    e = np.stack([-p[..., 1], p[..., 0], np.zeros_like(p[..., 0])], axis=-1)  # z × p
+    ne = np.linalg.norm(e, axis=-1, keepdims=True)
+    polar = ne[..., 0] < 1e-12
+    e = np.where(polar[..., None], np.array([0.0, 1.0, 0.0]), e / np.maximum(ne, 1e-300))
+    n = np.cross(p, e)
+    return e, n
+
+
+def tangent_dir(p, q):
+    """Unit tangent at p pointing along the great circle toward q."""
+    t = q - np.sum(p * q, axis=-1, keepdims=True) * p
+    return t / np.maximum(np.linalg.norm(t, axis=-1, keepdims=True), 1e-15)
+
+
+def gc_dist_km(p, q, radius_km):
+    return radius_km * np.arccos(np.clip(np.sum(p * q, axis=-1), -1.0, 1.0))
+
+
+def motion_to_omega(lat, lon, azimuth_deg, speed_cm_yr, radius_km):
+    """Angular velocity (rad/yr) that moves the point (lat, lon) along azimuth at speed."""
+    p = latlon_to_xyz(np.array([lat]), np.array([lon]))
+    e, n = east_north(p)
+    az = np.radians(azimuth_deg)
+    v = (np.sin(az) * e[0] + np.cos(az) * n[0]) * (speed_cm_yr / 100.0)  # m/yr
+    return np.cross(p[0], v) / (radius_km * 1000.0)
+
+
+def velocity(p, omega, radius_km):
+    """Surface velocity (m/yr) of points p on plates with angular velocity omega (broadcast)."""
+    return np.cross(omega, p) * (radius_km * 1000.0)
+
+
+def rotate_about(p, v, ang):
+    """Rotate tangent vectors v about normals p by ang radians (counter-clockwise seen from outside)."""
+    ang = np.asarray(ang, dtype=np.float64)[..., None]
+    return v * np.cos(ang) + np.cross(p, v) * np.sin(ang)
+
+
+def azimuth_deg(center, p):
+    """Bearing (deg clockwise from north, [0, 360)) from center (3,) to points p (..., 3)."""
+    e, n = east_north(center[None])
+    t = tangent_dir(np.broadcast_to(center, p.shape), p)
+    return np.degrees(np.arctan2(t @ e[0], t @ n[0])) % 360.0
+
+
+def great_circle_point(p0, t0, dist_km, radius_km):
+    """Point reached from p0 moving dist_km along unit tangent t0."""
+    a = dist_km / radius_km
+    return np.cos(a) * p0 + np.sin(a) * t0
diff --git a/mapgen/testdata/world.toml b/mapgen/testdata/world.toml
new file mode 100644
index 0000000..e650e39
--- /dev/null
+++ b/mapgen/testdata/world.toml
@@ -0,0 +1,31 @@
+# Test fixture world (the unit tests' planet); not an example — see example/config/world.toml.
+[planet]
+radius_km = 12742.0
+gravity_g = 1.05
+day_hours = 31.149
+year_days = 216
+tilt_deg = 20.0
+sea_level_pressure_bar = 1.0
+scale_height_km = 8.0
+o2_fraction = 0.21
+
+[units]
+span_m = 0.331013
+moment_s = 2.403473
+league_km = 15.4437
+
+[build]
+seed = 1296
+res_dev = 4
+res_final = 5
+raster_width = 8192
+dev_raster_width = 2048
+preview_width = 2048
+land_fraction = 0.19
+
+[climate]
+hadley_edge_deg = 20.0
+ferrel_edge_deg = 55.0
+
+[sketch]
+moves = []
diff --git a/mapgen/testing.py b/mapgen/testing.py
new file mode 100644
index 0000000..b70e5fe
--- /dev/null
+++ b/mapgen/testing.py
@@ -0,0 +1,59 @@
+"""Test fixtures shared with tools built on mapgen (e.g. a viewer's tests): a tiny two-continent world at H3 res 2."""
+from __future__ import annotations
+
+import atexit
+import functools
+import re
+import shutil
+import tempfile
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+
+FIXTURE_TOML = Path(__file__).resolve().parent / "testdata" / "world.toml"
+
+TECT = """
+[[plate]]
+id = "a"
+seed = [10.0, -30.0]
+kind = "continental"
+motion = [90.0, 4.0]
+[[plate]]
+id = "b"
+seed = [10.0, 30.0]
+kind = "continental"
+motion = [270.0, 4.0]
+[[plate]]
+id = "c"
+seed = [0.0, 150.0]
+kind = "oceanic"
+motion = [0.0, 3.0]
+"""
+
+
+def small_world(tmp: Path) -> None:
+    (tmp / "config").mkdir()
+    w = FIXTURE_TOML.read_text()
+    w = w.replace("raster_width = 8192", "raster_width = 256").replace("preview_width = 2048", "preview_width = 128")
+    w = w.replace("dev_raster_width = 2048", "dev_raster_width = 256")
+    w = re.sub(r"moves = \[.*?\n\]", "moves = []", w, flags=re.S)
+    (tmp / "config" / "world.toml").write_text(w)
+    (tmp / "config" / "tectonics.toml").write_text(TECT)
+    (tmp / "sketch").mkdir()
+    lat = 90 - (np.arange(100) + 0.5) * 1.8
+    lon = (np.arange(200) + 0.5) * 1.8 - 180
+    LA, LO = np.meshgrid(lat, lon, indexing="ij")
+    land = ((np.abs(LA - 10) < 35) & (np.abs(LO) < 70)).astype(np.uint8) * 255
+    for n in ("land", "mountains", "desert", "rainforest", "trench"):
+        Image.fromarray(land if n == "land" else np.zeros_like(land), "L").save(tmp / "sketch" / f"{n}.png")
+
+
+@functools.lru_cache(maxsize=1)
+def built_world() -> Path:
+    from mapgen import pipeline as P
+    tmp = Path(tempfile.mkdtemp(prefix="worldgen-fixture-"))
+    atexit.register(shutil.rmtree, tmp, True)          # 5 MB per build: never leave them in /tmp
+    small_world(tmp)
+    P.build(tmp, 2, log=lambda m: None)
+    return tmp
diff --git a/mapgen/viewer_export.py b/mapgen/viewer_export.py
new file mode 100644
index 0000000..e6f51b4
--- /dev/null
+++ b/mapgen/viewer_export.py
@@ -0,0 +1,39 @@
+"""Viewer textures: colourised equirectangular JPEG layers + layers.json."""
+from __future__ import annotations
+
+import json
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+
+
+def hex_color(c) -> str:
+    return "#%02x%02x%02x" % tuple(int(v) for v in list(c)[:3])
+
+
+def categorical_legend(names, palette) -> dict:
+    return {"type": "categorical",
+            "items": [{"name": n, "color": hex_color(palette[i % len(palette)])} for i, n in enumerate(names)]}
+
+
+def continuous_legend(unit, lo, hi, stops) -> dict:
+    vals = np.linspace(lo, hi, len(stops))
+    return {"type": "continuous", "unit": unit,
+            "stops": [[round(float(v), 4), hex_color(s)] for v, s in zip(vals, stops)]}
+
+
+def write_layer(vdir: Path, L: dict) -> dict:
+    """Write one layer's JPEG; returns its layers.json entry."""
+    vdir.mkdir(parents=True, exist_ok=True)
+    Image.fromarray(np.ascontiguousarray(L["rgb"], dtype=np.uint8)).save(vdir / f"{L['id']}.jpg", quality=90)
+    return {"id": L["id"], "name": L["name"], "file": f"{L['id']}.jpg", "legend": L["legend"]}
+
+
+def write_index(vdir: Path, entries: list[dict]) -> None:
+    vdir.mkdir(parents=True, exist_ok=True)
+    (vdir / "layers.json").write_text(json.dumps(entries, indent=1))
+
+
+def write_all(vdir: Path, layers: list[dict]) -> None:
+    write_index(vdir, [write_layer(vdir, L) for L in layers])
diff --git a/mapgen/zones.py b/mapgen/zones.py
new file mode 100644
index 0000000..46386f7
--- /dev/null
+++ b/mapgen/zones.py
@@ -0,0 +1,38 @@
+"""Config O₂ / low-gravity zones: smooth, noise-warped discs in mask units
+(O₂ ×(1 + 0.5 v), gravity ×(1 + 0.7 v)), added to the painted masks in the `sketch` stage. Very gradual: at most
+≈ 1.9 · |effect| / radius per km (≤ 0.23 O₂ points and ≤ 0.014 g per 30 km for the configured zones)."""
+from __future__ import annotations
+
+import numpy as np
+
+from .noise import fbm, name_seed
+from .sphere import gc_dist_km, latlon_to_xyz
+
+WARP = 0.2           # outline wobble: distances stretched or shrunk by up to 20 % (low-frequency noise)
+WARP_FREQ = 1.0
+
+
+def smootherstep(t):
+    t = np.clip(t, 0.0, 1.0)
+    return t * t * t * (t * (6.0 * t - 15.0) + 10.0)
+
+
+def contribution(xyz, zone: dict, seed: int, radius_km: float):
+    """v · (1 − smootherstep(d′ / r)), d′ = d · (1 + WARP · warp): v at the centre, 0 beyond r / (1 − WARP)."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    d = gc_dist_km(xyz, latlon_to_xyz(*zone["center"]), radius_km)
+    r = float(zone["radius_km"])
+    out = np.zeros(len(xyz))
+    near = d < r / (1.0 - WARP)
+    if near.any():
+        w = np.clip(2.0 * fbm(xyz[near], name_seed(seed, zone["name"]), 3, WARP_FREQ), -1.0, 1.0)
+        out[near] = zone["v"] * (1.0 - smootherstep(d[near] * (1.0 + WARP * w) / r))
+    return out
+
+
+def masks(xyz, zones: list, seed: int, radius_km: float):
+    """(m_o2, m_gravity): the config zones summed (the caller adds the painted masks and clips to [−1, 1])."""
+    m = {"o2": np.zeros(len(xyz)), "gravity": np.zeros(len(xyz))}
+    for z in zones:
+        m[z["field"]] += contribution(xyz, z, seed, radius_km)
+    return m["o2"], m["gravity"]