"""Deep-zoom tiles for the map viewer: procedural detail below the build's resolution.

Deterministic, seeded and continuous across tiles and zoom levels, but invented rather than data. Lives outside
mapgen/ so changing it never invalidates the build cache. Zoom z has 2^(z+1) × 2^z tiles of 256 px (x east from
−180°, y south from +90°); pixel size = 360° / (2^(z+1)·256).
"""
from __future__ import annotations

import hashlib
import io
import json
import multiprocessing
import os
import queue
import shutil
import signal
import tempfile
import threading
from collections import OrderedDict
from pathlib import Path

import numpy as np
from PIL import Image
from scipy.ndimage import binary_dilation, map_coordinates
from scipy.spatial import cKDTree

import worldgen_path  # noqa: F401  (mapgen on sys.path)
from mapgen import render as RN
from mapgen.noise import value_noise
from mapgen.sphere import east_north
from mapgen.sphere import latlon_to_xyz

import refine
import rivers as RV
import servecache

TILE = 256
MAX_Z = 16
VERSION = 6                                    # 2: river valleys; 3: per-valley surfaces; 4: 129² heights;
                                               # 5: meshes with true heights + water, era sharing; 6: erosion gullies
DRAW = 5       # how tiles draw a build, apart from VERSION (which refine's area keys hash): bump → tiles re-render,
               # refined areas stay; 1: river banks never below their water; 2: micro-relief + land colour variation;
               # 3: stronger (visible from the ground view), lush patches green not teal; 4: drawn from the water
               # surface raster, world lakes flat at their level (no pits in 3D); 5: land colour = terrain shading +
               # vegetation mosaic (author 2026-10-03)
TOP_KM, MIN_KM, GAIN = 8.0, 0.005, 0.55        # detail octaves: 8 km wavelength down to 5 m
MICRO_KM, MICRO_GAIN = 2.0, 0.75               # drawn-only micro-relief (not in relief_at): 2 km down to 5 m, rougher
MICRO_M = (10.0, 0.05)                         # its 2-km amplitude (m) = a + b × landform amplitude: swells on plains
TINT_KM, TINT = 4.0, (0.18, 22.0)              # land colour variation: 4 km down to the pixel; ± brightness, ± tint
WARP_KM, WARP_TOP_KM = 12.0, 48.0              # ≈ ⅓ of a res-5 cell: organic borders between cells
LAND_BORDER = 8                                # px around a relief tile for its land colour (slope, water)
CREST_M = 40.0                                 # crests/hollows: height above the smooth base surface, in units of
                                               # max(this, ½ × the landform's detail amplitude)
STEEP = 0.35                                   # slope (m/m) ≈ full "steep"
RIPARIAN_KM = 0.4                              # woods and green along water fade over this distance
VEG_KM, VEG_WARP = 16.0, 0.8                   # vegetation patches: noise from 16 km down, domain-warped (Quilez)
VEG_SD = 0.455                                 # the mosaic noise is ≈ normal with this SD at every zoom (measured on
                                               # r4, 2026-10-03): woods where it is below SD·Φ⁻¹(cover)
VEG_COVER = (("polar desert", 0.0), ("dry tundra", 0.02), ("moist tundra", 0.06), ("wet tundra", 0.1),
             ("rain tundra", 0.12), ("desert scrub", 0.08), ("desert", 0.02), ("dry scrub", 0.12),
             ("thorn steppe", 0.12), ("steppe", 0.12), ("thorn woodland", 0.3), ("very dry forest", 0.4),
             ("dry forest", 0.5), ("moist forest", 0.65), ("wet forest", 0.75), ("rain forest", 0.82))
GROUND_COVER = {"salt flat": 0.0, "mangrove": 0.9}  # wooded share of open land by Holdridge zone (first suffix match)
GULLY_KM, GULLY_GAIN, GULLY_AMP = 4.0, 0.5, 0.55  # erosion gullies: waves 4 km down, × detail amplitude, halving
GULLY_SLOPE = 0.15                             # slope (m/m) where gullies are ≈ ¾ on (tanh)
GULLY_SHARE = 0.8                              # on full slopes, this share of the finer fractal detail gives way
PH_JIT = 0.8                                   # gully lattice points jitter over this share of a cell
GULLY_STEP_KM = 0.5                            # finite-difference step for the downhill direction
RIDGE_MEAN = 0.394                             # E[1 − 2|n|] of value noise: keeps ridged detail unbiased
DEEP = {"ocean": (150, 0.0), "plain": (25, 0.0), "hills": (120, 0.0), "mountains": (450, 1.0), "plateau": (60, 0.0),
        "rift valley": (150, 0.3), "escarpment": (250, 0.5), "volcanic arc": (350, 1.0), "volcanic massif": (350, 1.0),
        "basalt plateau": (60, 0.0), "dunes": (30, 0.0), "badlands": (120, 0.3)}   # 8-km amplitude (m), ridged share
CONT = {"temperature": "T_mean", "rainfall": "P_ann", "o2": "po2", "gravity": "gravity", "pressure": "pressure",
        "fire": "fire", "bottom_temp": "bottom_temp", "sediment": "sediment", "vent_potential": "vent_potential",
        "currents": "current_speed", "sst": "sst", "productivity": "productivity"}
CAT = {"biomes": ("holdridge", "holdridge"), "seasonality": ("seasonality", "seasonality"), "landform": ("landform", "landform"),
       "ground": ("ground", "ground"), "ice": ("ice", "ice"), "plates": ("plate", "plates"),
       "seabed": ("seabed_type", "seabed_type"), "minerals": ("seabed_mineral", "seabed_mineral"),
       "deposits": ("deposit_main", "deposits")}
LAYERS = ("relief", "seafloor", "elevation", *CONT, *CAT)
MEM_TILES = 512
REGION_WARP = 0.4                              # refined areas: the organic warp at ≈ their 5 km cell scale:
REGION_FREQ = 7.0                              # 0.4 × the displacement, waves 7× shorter
MESH_N = 129                                   # heights per tile side for 3D patches (128 segments, corners included)
MESH_BYTES = MESH_N * MESH_N * 5               # float32 heights + uint8 water
EXT = {"height": "png", "mesh": "bin"}         # every other layer is a JPEG image
CAT_KEYS = ("holdridge", "ground", "ice", "landform", "seasonality", "plate", "lake", "seabed_type", "seabed_mineral", "deposit_main")
SEABED_RGB = np.array([[0, 0, 0], [150, 120, 80], [225, 225, 205], [150, 110, 95], [70, 70, 80], [95, 60, 70],
                       [120, 130, 110], [230, 110, 40], [40, 40, 60]], dtype=np.float64)
# by mapgen.seabed code: none, terrigenous, carbonate, clay, basalt, volcanic, continental, vents, trench
FLOOR_RAMP = [[5, 15, 45], [20, 60, 120], [90, 150, 190], [170, 215, 235]]   # −8,000 m … 0 m


def tile_ok(z: int, x: int, y: int) -> bool:
    return 0 <= z <= MAX_Z and 0 <= x < 2 ** (z + 1) and 0 <= y < 2 ** z


def seafloor_rgb(z, hs, seabed, land):
    """The sea floor as if the water were gone: depth ramp tinted 35 % by sea-floor type, hillshaded like land;
    land as the relief layer's colours would shade it (here: a plain earth tone, hillshaded)."""
    z = np.asarray(z, dtype=np.float64)
    depth = RN._ramp(z, -8000.0, 0.0, FLOOR_RAMP).astype(np.float64)
    tint = SEABED_RGB[np.clip(np.asarray(seabed), 0, len(SEABED_RGB) - 1)]
    sea = 0.65 * depth + 0.35 * tint
    earth = RN._ramp(z, 0.0, 4000.0, [[110, 125, 80], [150, 135, 100], [225, 225, 225]]).astype(np.float64)
    rgb = np.where(np.asarray(land)[..., None], earth, sea)
    return np.clip(rgb * (0.45 + 0.55 * np.asarray(hs))[..., None], 0, 255).astype(np.uint8)


def _ll(p):
    """(lat, lon) in degrees of unit vectors (N, 3)."""
    return np.degrees(np.arcsin(np.clip(p[:, 2], -1.0, 1.0))), np.degrees(np.arctan2(p[:, 1], p[:, 0]))


def _hash32(x, y, z, seed: int) -> np.ndarray:
    """Integer lattice coordinates → 30 well-mixed bits (int64), deterministic."""
    h = (x * 73856093) ^ (y * 19349663) ^ (z * 83492791) ^ (seed * 2654435761)
    h = h & 0xFFFFFFFF
    h = ((h ^ (h >> 15)) * 0x2C1B3C6D) & 0xFFFFFFFF
    h = ((h ^ (h >> 12)) * 0x297A2D39) & 0xFFFFFFFF
    return (h ^ (h >> 15)) & 0x3FFFFFFF


def _make_phacelle_jit():
    """TileSource._phacelle's lattice blend compiled (numba optional; WORLDGEN_NO_JIT=1 turns it off): the same
    hash, weights and products per point. exp may differ from numpy's in the last bit (≪ 1e-9 m of height)."""
    if os.environ.get("WORLDGEN_NO_JIT"):
        return None
    try:
        import numba
    except ImportError:
        return None

    jit = PH_JIT

    @numba.njit(cache=True, error_model="numpy")
    def phacelle(q, a, seed):
        n = q.shape[0]
        cos_out, sin_out = np.empty(n), np.empty(n)
        two_pi = 2.0 * np.pi
        st = np.empty((3, 2))
        pw = np.empty((3, 4, 2))
        sd = np.int64(seed) * np.int64(2654435761)
        for k in range(n):
            ax, ay, az = a[k, 0], a[k, 1], a[k, 2]
            fx, fy, fz = np.floor(q[k, 0]), np.floor(q[k, 1]), np.floor(q[k, 2])
            ix, iy, iz = np.int64(fx), np.int64(fy), np.int64(fz)
            tx, ty, tz = q[k, 0] - fx, q[k, 1] - fy, q[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)
            ph0 = two_pi * ((fx * ax + fy * ay + fz * az) - 0.5 * jit * (ax + ay + az))
            e0r, e0i = np.cos(ph0), np.sin(ph0)
            for d in range(3):                              # step phasors: corner step, jitter step and its powers
                ang = two_pi * a[k, d]
                st[d, 0], st[d, 1] = np.cos(ang), np.sin(ang)
                ang = two_pi * jit / 3.0 * a[k, d]
                jr, ji = np.cos(ang), np.sin(ang)
                pw[d, 0, 0], pw[d, 0, 1] = 1.0, 0.0
                for m in range(1, 4):
                    pw[d, m, 0] = pw[d, m - 1, 0] * jr - pw[d, m - 1, 1] * ji
                    pw[d, m, 1] = pw[d, m - 1, 0] * ji + pw[d, m - 1, 1] * jr
            re, im = 0.0, 0.0
            for dx in range(2):
                wx = tx if dx else 1.0 - tx
                for dy in range(2):
                    wy = wx * (ty if dy else 1.0 - ty)
                    for dz in range(2):
                        w = wy * (tz if dz else 1.0 - tz)
                        h = ((ix + dx) * 73856093) ^ ((iy + dy) * 19349663) ^ ((iz + dz) * 83492791) ^ sd
                        h = h & 0xFFFFFFFF
                        h = ((h ^ (h >> 15)) * 0x2C1B3C6D) & 0xFFFFFFFF
                        h = ((h ^ (h >> 12)) * 0x297A2D39) & 0xFFFFFFFF
                        h = (h ^ (h >> 15)) & 0x3FFFFFFF
                        zr, zi = e0r, e0i
                        if dx:
                            zr, zi = zr * st[0, 0] - zi * st[0, 1], zr * st[0, 1] + zi * st[0, 0]
                        if dy:
                            zr, zi = zr * st[1, 0] - zi * st[1, 1], zr * st[1, 1] + zi * st[1, 0]
                        if dz:
                            zr, zi = zr * st[2, 0] - zi * st[2, 1], zr * st[2, 1] + zi * st[2, 0]
                        m0, m1, m2 = h & 3, (h >> 2) & 3, (h >> 4) & 3
                        zr, zi = zr * pw[0, m0, 0] - zi * pw[0, m0, 1], zr * pw[0, m0, 1] + zi * pw[0, m0, 0]
                        zr, zi = zr * pw[1, m1, 0] - zi * pw[1, m1, 1], zr * pw[1, m1, 1] + zi * pw[1, m1, 0]
                        zr, zi = zr * pw[2, m2, 0] - zi * pw[2, m2, 1], zr * pw[2, m2, 1] + zi * pw[2, m2, 0]
                        re += w * zr
                        im += w * zi
            A = two_pi * (q[k, 0] * ax + q[k, 1] * ay + q[k, 2] * az)
            ca, sa = np.cos(A), np.sin(A)
            r = max(np.hypot(re, im), 1e-12)
            re, im = re / r, im / r
            cos_out[k] = ca * re + sa * im
            sin_out[k] = sa * re - ca * im
        return cos_out, sin_out
    return phacelle


_phacelle_jit = _make_phacelle_jit()


def _resolved(wave_km, px_km) -> float:
    """0 → 1 as a pattern of this wavelength grows from 4 to 8 pixels: too small to read → drawn as its mean."""
    return float(np.clip(np.log2(wave_km / (4.0 * px_km)), 0.0, 1.0))


def veg_cover(zones) -> np.ndarray:
    """Wooded share of open land per Holdridge zone name (VEG_COVER: the first suffix that matches)."""
    return np.array([next((c for k, c in VEG_COVER if n.endswith(k)), 0.0) for n in zones], dtype=np.float64)


def coast_water(z, sea, a: float = 20.0):
    """Water where the detailed height is below a threshold that runs smoothly from −∞ (no sea nearby) through ≈0
    (mixed coast: the build's own shoreline, z = 0) to +∞ (open sea): fractal coasts, no hard cut where a sea cell
    leaves the neighbourhood, no inland seas in dry depressions."""
    s = np.clip(sea, 1e-6, 1.0 - 1e-6)
    return np.asarray(z) < a * (1.0 / (1.0 - s) - 1.0 / s)


def tile_lat_lon(z, x, y, idx):
    """Lat/lon (deg) of tile-local pixel-centre indices (negative or ≥ TILE reach into the neighbours)."""
    f = (np.asarray(idx, dtype=np.float64) + 0.5) / TILE
    span = 180.0 / 2 ** z
    return 90.0 - (y + f) * span, -180.0 + (x + f) * span


def _render_worker(conn, sources, parent: int) -> None:
    """A forked render process: renders what it is sent for the source it names; a job for an older region set is
    answered "stale" (the parent renders it itself while it forks fresh workers)."""
    signal.signal(signal.SIGINT, signal.SIG_IGN)               # Ctrl-C is for the server; when it dies, the pipe
    if os.getppid() != parent:                                  # ends (EOF below) and so does this worker
        os._exit(0)                                             # (not PR_SET_PDEATHSIG: that follows the forking
    for s in sources.values():                                  # thread, which may be a short-lived one)
        s.lock = threading.Lock()                               # forked while a server thread may have held one
    fd = conn.fileno()                                          # inherited: the server's listening socket, the other
    os.closerange(3, fd)                                        # workers' pipes… — holding them would keep a dead
    os.closerange(fd + 1, os.sysconf("SC_OPEN_MAX"))            # server's port and workers alive
    while True:
        try:
            name, layer, z, x, y, fp = conn.recv()
        except (EOFError, OSError):
            return
        try:
            src = sources[name]
            if src.regions.fingerprint != fp:
                conn.send(("stale", None))
                continue
            if layer == EXPORT_JOB:
                import export
                conn.send(("ok", export.tile_data(src, z, x, y)))
            else:
                conn.send(("ok", src.render(layer, z, x, y)))
        except Exception as e:                                  # noqa: BLE001 — reported to the parent
            conn.send(("err", f"{type(e).__name__}: {e}"))


RENDER_TIMEOUT_S = 300.0   # one job; a worker slower than this is taken as hung and stopped


def fingerprint(out: Path, seed: int) -> str:
    """A build's identity for its tiles and refined areas (a copy must keep the files' mtimes: rsync -a does)."""
    out = Path(out)
    h = hashlib.sha1(f"{VERSION}:{seed}".encode())
    for name in ("fields.json", "cells_meta.json"):
        h.update((out / name).read_bytes())
    for name in ("cells.npz", "raster/elevation.png"):
        st = (out / name).stat()
        h.update(f"{name}:{st.st_size}:{st.st_mtime_ns}".encode())
    return h.hexdigest()[:10]


class Busy(Exception):
    """Too many tiles already wait for a render (a public server's limit): ask again later."""


class RenderPool:
    """Worker processes forked from a warmed TileSource (the loaded world is shared copy-on-write), so tiles render
    in parallel on several cores. Each job names the regions fingerprint it expects."""
    def __init__(self, sources, n: int):
        sources = sources if isinstance(sources, dict) else {sources.name: sources}
        self.sources = sources
        import export  # noqa: F401 — loaded before forking: no import in a child (a parent thread may hold its lock)
        ctx = multiprocessing.get_context("fork")
        self.idle, self.procs, self.alive, self.lock = queue.Queue(), [], 0, threading.Lock()
        self.proc_of = {}
        for _ in range(n):
            a, b = ctx.Pipe()
            p = ctx.Process(target=_render_worker, args=(b, sources, os.getpid()), daemon=True)
            p.start()
            b.close()
            self.procs.append(p)
            self.proc_of[a] = p
            self.idle.put(a)
            self.alive += 1

    def render(self, name, layer, z, x, y, fp) -> bytes:
        while True:                                             # wait for a free worker, unless none are left
            if self.alive <= 0:
                raise RuntimeError("no render workers left")
            try:
                conn = self.idle.get(timeout=2.0)
                break
            except queue.Empty:
                continue
        try:
            conn.send((name, layer, z, x, y, fp))
            if not conn.poll(RENDER_TIMEOUT_S):
                self.proc_of[conn].kill()
                with self.lock:
                    self.alive -= 1
                raise RuntimeError(f"a render worker timed out ({layer} {z}/{x}/{y})")
            kind, data = conn.recv()
        except (EOFError, OSError):
            with self.lock:
                self.alive -= 1                                 # a dead worker is not handed out again
            raise RuntimeError("a render worker died")
        self.idle.put(conn)
        if kind != "ok":
            raise RuntimeError(data or kind)
        return data

    def close(self) -> None:
        for p in self.procs:
            p.terminate()
        for p in self.procs:
            p.join(timeout=5)
        self.alive = 0


REGION_VERSION = 4   # how tiles draw refined areas (bump: their saved tiles re-render); 2: lakes end at their bed;
                     # 3: judged on the cells' smooth ground (no dry pits); 4: that ground against the lowest bed of
                     # the lake's cells it mixes (no dry pits where a lake deepens), shores follow the ground
EXPORT_JOB = "__export__"   # a render-pool job: a tile's export layers (export.tile_data) instead of an image
CARRY_KM = 100.0   # a region change keeps saved tiles farther than this from the areas that changed
MANIFEST = "areas.npz"   # in a region set's tile folder: its areas (name, key, coarse points) for a later carry-over
MANIFEST_KM = 25.0       # the coarse points' grid


class TileSource:
    def __init__(self, world, seed: int, cache_dir: Path | None = None, regions_dir: Path | None = None,
                 with_regions: bool = True, name: str = "base", share=None):
        self.w, self.seed = world, int(seed)
        self.out = world.out
        self.R = float(world.cells_meta["radius_km"])
        self.enc = json.loads((self.out / "fields.json").read_text())["continuous"]
        self.fingerprint = self._fingerprint()
        self.regions_dir = Path(regions_dir or (self.out / "regions"))
        self.regions = (refine.RegionSet(self.regions_dir, self.fingerprint, world.legends, self.R) if with_regions
                        else refine.RegionSet.none())               # refined areas (a build needs none)
        self.base_tag = f"v{VERSION}.{DRAW}-{self.fingerprint}"
        self.tag = self.base_tag + self._region_tag(self.regions)   # new build or new refined
        # areas: never serve old tiles. On disk, tiles away from the areas stay in the world's own cache (base_tag).
        self.url = f"/tiles/{self.tag}/{{layer}}/{{z}}/{{x}}/{{y}}.jpg"
        self.own_cache = cache_dir is None
        self.cache_root = cache_dir or (self.out / "tiles")
        self.cache_dir = self.cache_root / self.tag
        self.lock = threading.Lock()
        self._rasters, self._tree, self._net, self.mem = {}, None, None, OrderedDict()
        self.name, self.share = name, share
        self.peers = {name: self}
        self._mask_tree = None                                  # an era: its influence mask (cells.npz era_mask)
        if share is not None:
            with np.load(self.out / "cells.npz") as a:
                m = a["era_mask"].astype(bool) if "era_mask" in a.files else np.ones(len(a["g_ids"]), bool)
                self._mask_tree = cKDTree(np.asarray(a["g_xyz"], dtype=np.float64)[m]) if m.any() else None
            self._spacing_km = float(np.sqrt(4.0 * np.pi * self.R ** 2 / len(self.w.ids)))
        self.pool = None                                        # RenderPool (start_workers), else render here
        self.cap = None                                         # tilecap.TileCap: a size cap on the saved tiles
        names = world.legends["landform"]
        self.amp = np.array([DEEP.get(n, (60, 0.0))[0] for n in names], dtype=np.float64)
        self.ridge = np.array([DEEP.get(n, (60, 0.0))[1] for n in names], dtype=np.float64)

    def _fingerprint(self) -> str:
        return fingerprint(self.out, self.seed)

    @staticmethod
    def _region_tag(rs) -> str:
        """The refined areas' part of the tile tag: their results and how tiles draw them (REGION_VERSION)."""
        return "" if rs.empty else hashlib.sha1(f"{rs.fingerprint}:{REGION_VERSION}".encode()).hexdigest()[:6]

    def reload_regions(self, cleanup: bool = True, refork: bool = True) -> None:
        """A region build finished: new refined areas, new tile URLs (old ones stop answering), fresh memory cache.
        Saved tiles of areas that did not change (and are far from the ones that did) move to the new set."""
        rs = refine.RegionSet(self.regions_dir, self.fingerprint, self.w.legends, self.R)
        tag = self.base_tag + self._region_tag(rs)
        with self.lock:
            old, old_rs = self.tag, self.regions
            self.regions, self.tag = rs, tag
            self.url = f"/tiles/{tag}/{{layer}}/{{z}}/{{x}}/{{y}}.jpg"
            self.cache_dir = self.cache_root / tag
            self.mem.clear()
        gone = [h for h in old_rs.area_names if h not in rs.area_names]
        new = [h for h in rs.area_names if old_rs.area_keys.get(h) != rs.area_keys[h]]   # new or rebuilt
        gone_trees = {h: old_rs.area_trees[h] for h in gone if h in old_rs.area_trees}   # all the carry-over needs
        old_spacing = getattr(old_rs, "spacing_km", 0.0)
        del old_rs                                              # the old set's memory goes before the fork
        rs.warm()                                               # before the fresh workers are forked
        if cleanup:
            self.prune_region_caches()                          # the previous set's (also when no area is left)
        if cleanup and refork:
            self.refork()
        if not cleanup or old in (tag, self.base_tag):
            return
        self._carry(self.cache_root / old, rs, tag, list(gone_trees.values()), old_spacing, new)
        self._save_manifest(rs, tag)

    def _carry(self, src_dir: Path, rs, tag: str, gone_trees: list, gone_spacing: float, new: list) -> None:
        """Move src_dir's saved tiles away from the gone and the new (or rebuilt) areas to region set `tag`; the
        rest of src_dir goes."""
        if not rs.empty and src_dir.is_dir():
            for f in src_dir.rglob("*"):
                if not f.is_file() or f.name.endswith(".part"):
                    continue
                rel = f.relative_to(src_dir)                    # layer/z/x/y.ext
                try:
                    z, x, y = int(rel.parts[1]), int(rel.parts[2]), int(rel.stem)
                except (IndexError, ValueError):
                    continue
                if (refine.touches_trees(gone_trees, gone_spacing, self.R, z, x, y, CARRY_KM)
                        or rs.touches(z, x, y, new, CARRY_KM)):
                    continue
                dst = self.cache_root / tag / rel
                dst.parent.mkdir(parents=True, exist_ok=True)
                os.replace(f, dst)
        shutil.rmtree(src_dir, ignore_errors=True)

    def _save_manifest(self, rs, tag: str) -> None:
        """Region set `tag`'s areas (name, key, points on a MANIFEST_KM grid) next to its tiles: a later start with
        other areas carries the tiles over (carry_stale) though these areas' results are gone by then."""
        f = self.cache_root / tag / MANIFEST
        if rs.empty or f.exists():
            return
        xyz, area, q = np.asarray(rs.arrays["g_xyz"]), np.asarray(rs.area), MANIFEST_KM / self.R
        pts = {f"p{k}": (np.unique(np.round(xyz[area == k] / q).astype(np.int32), axis=0) * q).astype(np.float32)
               for k in range(len(rs.area_names))}
        f.parent.mkdir(parents=True, exist_ok=True)
        tmp = f.with_name(f".{MANIFEST}-{os.getpid()}.npz")
        np.savez(tmp, names=np.array(rs.area_names), keys=np.array([rs.area_keys[h] for h in rs.area_names]),
                 version=REGION_VERSION, **pts)
        os.replace(tmp, f)

    def carry_stale(self) -> None:
        """At a start: saved tiles of this world's earlier region sets (regions rebuilt while the server was stopped)
        move to the current set where their areas did not change, as a reload does; then the current set's
        manifest. Sets without a manifest, or drawn by another REGION_VERSION, are left to warm()'s cleanup."""
        rs = self.regions
        if rs.empty:
            return
        for d in sorted(self.cache_root.glob(self.base_tag + "?*")):
            if d.name == self.tag or not (d / MANIFEST).is_file():
                continue
            try:
                with np.load(d / MANIFEST) as m:
                    if int(m["version"]) != REGION_VERSION:       # drawn the old way: left to warm()'s cleanup
                        continue
                    names = [str(n) for n in m["names"]]
                    keys = dict(zip(names, (str(k) for k in m["keys"])))
                    gone = [cKDTree(m[f"p{k}"].astype(np.float64)) for k, h in enumerate(names)
                            if rs.area_keys.get(h) != keys[h]]
            except (OSError, ValueError, KeyError):
                continue
            new = [h for h in rs.area_names if keys.get(h) != rs.area_keys[h]]
            self._carry(d, rs, self.tag, gone, MANIFEST_KM, new)
        self._save_manifest(rs, self.tag)

    def refork(self) -> None:
        """Fresh workers forked from the new state share its memory (loading it in each worker: six copies)."""
        if self.pool is None:
            return
        import gc
        gc.collect()
        refine.trim_memory()
        retired = self.pool
        fresh = RenderPool(self.peers, len(retired.procs))
        for s in self.peers.values():
            s.pool = fresh
        retired.close()

    def prune_region_caches(self) -> None:
        """Remove this world's region serve caches other than the current set's (all of them when it has no areas)."""
        keep = [self.regions.cache] if self.regions.cache is not None else []
        servecache.prune(refine.results_dir(self.regions_dir, self.fingerprint), keep)

    def meta(self) -> dict:
        return {"url": self.url, "mesh_url": f"/tiles/{self.tag}/mesh/{{z}}/{{x}}/{{y}}.bin", "min_z": 5, "max_z": MAX_Z,
                "rivers": True}

    # --- data -------------------------------------------------------------------------------------------------
    def raster(self, name: str) -> np.ndarray:
        with self.lock:
            if name not in self._rasters:
                e = self.enc[name]
                a = np.asarray(Image.open(self.out / "raster" / e["file"]), dtype=np.float32)
                self._rasters[name] = a * np.float32(e["scale"]) + np.float32(e["offset"])
            return self._rasters[name]

    def surface(self) -> np.ndarray:
        """The raster the tiles draw from: the water surface (lakes at their level) where the build has one."""
        return self.raster("surface" if "surface" in self.enc else "elevation")

    def _flat_lakes(self, z, lat, lon, e, nn):
        """The world's heights with its lakes flat at their level (refined areas blend their own in over them)."""
        a = self.w.arrays
        if "lake_level_m" not in a:
            return z
        lev = np.asarray(a["lake_level_m"], dtype=np.float64)[self._warped_cells(lat, lon, e, nn).reshape(np.shape(z))]
        return np.where(np.isfinite(lev), lev, z)

    def tree(self) -> cKDTree:
        with self.lock:
            if self._tree is None:
                self._tree = cKDTree(np.asarray(self.w.arrays["g_xyz"], dtype=np.float64))
            return self._tree

    @property
    def river_net(self) -> RV.RiverNet:
        with self.lock:
            if self._net is None:
                self._net = RV.RiverNet(self.w.arrays, self.w.legends, self.R)
            return self._net

    def warm(self) -> None:
        self.tree()
        self.river_net
        self.regions.warm()
        self.surface()
        if self.own_cache:                                    # tiles of older builds / models are dead weight; other
            for d in (self.out / "tiles").glob("v*"):         # region sets of this world only when we have our own
                if not d.is_dir() or d.name in (self.tag, self.base_tag):   # (a smoke server has none: it keeps them)
                    continue
                if not d.name.startswith(self.base_tag) or not self.regions.empty:
                    shutil.rmtree(d, ignore_errors=True)

    def bilinear(self, a, lat, lon):
        H, W = a.shape
        col = (np.asarray(lon) + 180.0) / 360.0 * W - 0.5
        row = np.clip((90.0 - np.asarray(lat)) / 180.0 * H - 0.5, 0, H - 1)
        c0 = np.floor(col).astype(np.int64)
        r0 = np.minimum(np.floor(row).astype(np.int64), H - 2)
        fc, fr = col - c0, row - r0
        ca, cb = c0 % W, (c0 + 1) % W
        top = a[r0, ca] * (1 - fc) + a[r0, cb] * fc
        bot = a[r0 + 1, ca] * (1 - fc) + a[r0 + 1, cb] * fc
        return top * (1 - fr) + bot * fr

    # --- procedural detail ------------------------------------------------------------------------------------
    def _octaves(self, xyz, px_km, top_km, min_km, gain, seed, ridge=None):
        """Σ gainᵏ·noiseₖ for wavelengths top_km/2ᵏ ≥ max(px_km, min_km); the finest ones fade in over 1–4 px."""
        out = np.zeros(len(xyz))
        a, lam, k = 1.0, top_km, 0
        while lam >= max(px_km, min_km):
            w = min(1.0, 0.5 * float(np.log2(lam / px_km)))    # fade in over 1–4 px: no pixel-scale grain
            v = value_noise(xyz * (self.R / lam), seed + 101 * k)
            if ridge is not None:
                v = (1.0 - ridge) * v + ridge * (1.0 - 2.0 * np.abs(v) - RIDGE_MEAN)
            out += w * a * v
            a *= gain
            lam /= 2.0
            k += 1
        return out

    def _detail(self, xyz, px_km, amp, ridge, lat, lon, base=None):
        """Procedural relief (m) at points: the fractal detail, and on slopes erosion gullies in its place — stripes
        running downhill, octave by octave, each finer one following the slopes the coarser ones made, so they
        branch (after Rune Skovbo Johansen's erosion filter and Phacelle noise, 2025). The downhill direction is that
        of the bedrock plus the two coarsest detail octaves, the same at every zoom. base(lat, lon): bedrock height
        (m), default the build's elevation raster. Each gully octave fades in as its waves grow from 4 to 8 px."""
        amp, ridge = np.asarray(amp, dtype=np.float64), np.asarray(ridge, dtype=np.float64)
        det = self._octaves(xyz, px_km, TOP_KM, MIN_KM, GAIN, self.seed + 7001, ridge)
        fz = _resolved(GULLY_KM, px_km)
        if fz == 0.0:
            return amp * det
        base = base or (lambda la, lo: self.bilinear(self.raster("elevation"), la, lo))
        ref = min(px_km, 0.25)                                       # zoom-independent coarse octaves

        def coarse(p):
            return self._octaves(p, ref, TOP_KM, TOP_KM / 2, GAIN, self.seed + 7001, ridge)

        e, n = east_north(xyz)
        h = GULLY_STEP_KM / self.R
        c0 = coarse(xyz)
        z0 = base(lat, lon) + amp * c0
        grad = []
        for t in (e, n):                                             # forward differences, m per m
            hi = xyz + h * t
            hi = hi / np.linalg.norm(hi, axis=1, keepdims=True)
            grad.append((base(*_ll(hi)) + amp * coarse(hi) - z0) / (GULLY_STEP_KM * 1000.0))
        ge, gn = grad
        on = fz * np.tanh(np.hypot(ge, gn) / GULLY_SLOPE)
        out = amp * (det + GULLY_SHARE * on * (c0 - det))   # finer fractal detail gives way on slopes
        lam, a, k = GULLY_KM, GULLY_AMP, 0
        while lam >= max(4.0 * px_km, MIN_KM):
            w = fz * _resolved(lam, px_km)
            sl = np.hypot(ge, gn)
            ue, un = ge / np.maximum(sl, 1e-12), gn / np.maximum(sl, 1e-12)
            across = -un[:, None] * e + ue[:, None] * n                 # tangent, across the slope
            c, sn, dphase = self._phacelle(xyz, across, lam, self.seed + 7801 + 31 * k)
            prof = 1.0 - 2.0 * (1.0 - np.abs(c)) ** 1.5               # sharp gully floors, rounded spurs
            m = np.tanh(sl / GULLY_SLOPE)
            height = w * a * amp * m
            out += height * prof
            dprof = -3.0 * np.sqrt(np.maximum(1.0 - np.abs(c), 0.0)) * np.sign(c) * sn * dphase
            ge = ge + height * dprof * (-un)                            # the slope this octave made (m/m) steers
            gn = gn + height * dprof * ue                               # the next
            lam /= 2.0
            a *= GULLY_GAIN
            k += 1
        return out

    def _phacelle(self, xyz, across, lam_km, seed):
        """Stripes of wavelength lam_km across `across` at points: the phases from the 8 jittered lattice points
        around each point, blended as unit phasors (Johansen's Phacelle noise, in 3-D on the sphere) so the
        stripes stay sharp and continuous. Returns (cos, sin) of the blended phase and its rate across (rad per m).
        Phase at a lattice point c: 2π (q − c)·across = 2π q·a − 2π (i·a + corner·a + jitter·a). The jitter takes
        4 steps per axis (PH_JIT wide), so every corner's phasor is a product of a few per-point ones: 7 sin/cos
        per point instead of 18."""
        q = xyz * (self.R / lam_km)
        if _phacelle_jit is not None and q.dtype == np.float64:
            c, sn = _phacelle_jit(np.ascontiguousarray(q), np.ascontiguousarray(across, dtype=np.float64), int(seed))
            return c, sn, 2.0 * np.pi / (lam_km * 1000.0)
        i = np.floor(q)
        t = q - i
        t = t * t * (3.0 - 2.0 * t)
        a = across
        ii = i.astype(np.int64)
        tp = 2.0 * np.pi
        e0 = np.exp(1j * tp * ((i * a).sum(1) - 0.5 * PH_JIT * a.sum(1)))   # lattice corner (0,0,0), jitter −½
        ex, ey, ez = (np.exp(1j * tp * a[:, d]) for d in range(3))           # one corner step along x, y, z
        jx, jy, jz = (np.exp(1j * tp * PH_JIT / 3.0 * a[:, d]) for d in range(3))   # one jitter step (of 3)
        jpow = [[np.ones(len(q), complex), j, j * j, j * j * j] for j in (jx, jy, jz)]
        acc = np.zeros(len(q), complex)
        for dx in (0, 1):
            wx = t[:, 0] if dx else 1.0 - t[:, 0]
            for dy in (0, 1):
                wy = wx * (t[:, 1] if dy else 1.0 - t[:, 1])
                for dz in (0, 1):
                    w = wy * (t[:, 2] if dz else 1.0 - t[:, 2])
                    h = _hash32(ii[:, 0] + dx, ii[:, 1] + dy, ii[:, 2] + dz, seed)
                    z = e0 * (ex if dx else 1.0) * (ey if dy else 1.0) * (ez if dz else 1.0)
                    z = z * np.choose(h & 3, jpow[0]) * np.choose((h >> 2) & 3, jpow[1]) * np.choose((h >> 4) & 3, jpow[2])
                    acc += w * z
        re, im = acc.real, acc.imag
        A = tp * (q * a).sum(1)
        ca, sa = np.cos(A), np.sin(A)
        r = np.maximum(np.hypot(re, im), 1e-12)
        re, im = re / r, im / r                                      # unit phasor of −(blended lattice phase)
        return ca * re + sa * im, sa * re - ca * im, 2.0 * np.pi / (lam_km * 1000.0)

    def _amp(self, xyz):
        """Landform amplitude, ridged share and sea fraction blended over the 4 nearest cells. The kernel falls to 0 at
        the 5th-nearest distance, so the blend stays continuous when the set of neighbours changes."""
        d, k = self.tree().query(xyz, k=5)
        w = np.clip(1.0 - d[:, :4] / np.maximum(d[:, 4:5], 1e-12), 0.0, None) ** 2
        s = w.sum(axis=1, keepdims=True)
        w = np.where(s > 0, w / np.maximum(s, 1e-300), 0.25)
        k4 = k[:, :4]
        lf = np.asarray(self.w.arrays["landform"])[k4]
        sea = np.asarray(self.w.arrays["ocean"])[k4].astype(np.float64)
        return (w * self.amp[lf]).sum(1), (w * self.ridge[lf]).sum(1), (w * sea).sum(1)

    def _micro(self, xyz, px_km, amp):
        """Micro-relief (m) the tiles draw on top of relief_at: low swells that keep plains from looking poured."""
        return (MICRO_M[0] + MICRO_M[1] * amp) * self._octaves(xyz, px_km, MICRO_KM, MIN_KM, MICRO_GAIN, self.seed + 7301)

    @staticmethod
    def _roughen(t, m, V, base):
        """Heights t + micro-relief m with valleys V cut in, never digging dry ground below a valley's water (its
        base − 1 m) that t alone kept above it. base: +inf / ≥ 1e8 where no valley reaches."""
        z0 = np.minimum(t, V)
        lo = np.where(base < 1e8, base - 1.0, -np.inf)
        return np.maximum(np.minimum(t + m, V), np.minimum(z0, lo))

    def _warp_offsets(self, xyz, px_km):
        """East/north warp noise (≈ ±1) for the organic cell borders; octaves down to px_km."""
        return (self._octaves(xyz, px_km, WARP_TOP_KM, 0.05, 0.5, self.seed + 9001),
                self._octaves(xyz, px_km, WARP_TOP_KM, 0.05, 0.5, self.seed + 9377))

    def _region_warp(self, xyz, px_km):   # the organic warp for refined areas (waves at their cell scale)
        return self._warp_offsets(np.asarray(xyz) * REGION_FREQ, 4 * px_km * REGION_FREQ)

    def _warped_cells(self, lat, lon, e, nn):
        return self.tree().query(latlon_to_xyz(*self._warp_ll(lat, lon, e, nn)).reshape(-1, 3))[1]

    def _refined(self, xyz, lat, lon, px_km, rough, warp, lattice=None):
        """Heights inside refined areas: the fine cells' valley-shoulder heights (compact 4-of-5 kernel) + procedural
        detail below 5 km; the refined rivers' valleys cut in with their rock/climate walls; lakes wherever the ground
        lies below a lake's surface. The organic warp (REGION_WARP of the world's) shapes rivers and lake shores.
        Computed only where an area shows (weight > 0). Returns (z, weight, channel, draw, lake, sea fraction):
        channels wherever the area shows (they meet the world's rivers in the blend band), lakes where it dominates."""
        rs = self.regions
        w = rs.weight(xyz)
        n = len(xyz)
        z, ch, dr, lake, sea = np.zeros(n), np.zeros(n, bool), np.zeros(n, bool), np.zeros(n, bool), np.zeros(n)
        on = np.where(w > 0)[0]
        if not len(on):
            return z, w, ch, dr, lake, sea
        q = xyz[on]
        idx5, d5 = rs.nearest(q, k=5)
        wt = np.clip(1.0 - d5[:, :4] / np.maximum(d5[:, 4:5], 1e-12), 0.0, None) ** 2
        s = wt.sum(axis=1, keepdims=True)
        wt = np.where(s > 0, wt / np.maximum(s, 1e-300), 0.25)
        idx = idx5[:, :4]
        zq = (rs.z_env[idx] * wt).sum(1)
        smooth = zq.copy()                                      # the cells' ground, before detail and valleys
        sea[on] = (rs.arrays["ocean"][idx].astype(np.float64) * wt).sum(1)
        amp, ridge, _ = self._amp(q)
        zq = zq + amp * 0.55 * self._octaves(q, px_km, TOP_KM / 2, MIN_KM, GAIN, self.seed + 7001, ridge)
        m = self._micro(q, px_km, amp)
        wla, wlo = self._warp_ll(np.asarray(lat)[on], np.asarray(lon)[on], warp[0][on], warp[1][on], scale=REGION_WARP)
        wxyz = latlon_to_xyz(wla, wlo).reshape(-1, 3)
        if rs.net is not None and rs.net.tree is not None:
            rg = 0.0 if rough is None else np.asarray(rough)[on] if np.ndim(rough) else rough
            if lattice is None:
                V, c, d, base = rs.net.valleys(wxyz, px_km, rg, with_base=True)
            else:   # tiles: valley surfaces on the aligned 4-px lattice (seams agree, fast), channels per pixel
                cxyz_w, up, near = lattice
                base, wall_up, rf = (np.full(len(cxyz_w), np.inf), np.zeros(len(cxyz_w)), np.zeros(len(cxyz_w)))
                if near.any():
                    base[near], wall_up[near], rf[near] = rs.net.parts(cxyz_w[near], px_km, 4 * px_km)
                base = up(np.where(np.isfinite(base), base, 1e9)).ravel()[on]
                V = base + np.maximum(0.0, up(wall_up).ravel()[on] + up(rf).ravel()[on] * rg)
                c, d, lev = rs.net.channel(wxyz, px_km)
                V = np.where(c, lev, V)
            ch[on], dr[on] = c, d
            zq = self._roughen(zq, m, V, base)
        else:
            zq = zq + m
        near = rs.nearest(wxyz)[0]
        level, floor = rs.lake_level[near], rs.lake_floor[near]
        kl = rs.lake_level[idx]                                          # a dry cell's point by a lake: the nearest
        first = np.argmax(np.isfinite(kl), axis=1)                       # lake cell of its kernel (its shore follows
        level = np.where(np.isfinite(level), level, kl[np.arange(len(kl)), first])   # the ground, not cell borders)
        same = np.abs(kl - level[:, None]) < 0.5            # smooth mixes the kernel's cells: their
        floor = np.fmin(floor, np.where(same, rs.lake_floor[idx], np.inf).min(1))   # beds count too
        # water below the lake's surface, unless the cells' own ground (not the fine detail, which would punch dry
        # pits) lies far below the lake's bed: then this is past a dam or a cliff, not under the lake
        lk = (np.isfinite(level) & (zq <= np.nan_to_num(level, nan=-1e9))
              & (smooth >= np.nan_to_num(floor, nan=1e9)))
        z[on] = np.where(lk, level, zq)
        lake[on] = lk & (w[on] > 0.5)
        return z, w, ch, dr, lake, sea

    def _z_sea(self, lat, lon, px_km):
        """Detailed height (m) with river valleys cut in, neighbourhood sea fraction, unit vectors and the river
        channel at points (1-D float64 arrays)."""
        xyz = latlon_to_xyz(np.clip(lat, -90, 90), lon).reshape(-1, 3)
        amp, ridge, sea = self._amp(xyz)
        D = self._detail(xyz, px_km, amp, ridge, lat, lon)
        wla, wlo = self._warp_ll(lat, lon, *self._warp_offsets(xyz, 4 * px_km))
        z, channel, _ = self._valleys(wla, wlo, self.bilinear(self.surface(), lat, lon) + D, px_km,
                                      D, self._micro(xyz, px_km, amp))
        z = self._flat_lakes(z, lat, lon, *self._warp_offsets(xyz, 4 * px_km))
        if not self.regions.empty:   # refined areas: their own heights, rivers, lakes and sea, blended in at the edge
            zr, w, chr_, _, lk, sear = self._refined(xyz, lat, lon, px_km, D, self._region_warp(xyz, px_km))
            z = (1 - w) * z + w * zr
            channel = (channel & (w < 1)) | chr_ | lk           # the world's rivers meet the refined ones in the band
            sea = np.where(w > 0.5, sear, sea)
        return z, sea, xyz, channel

    def _warp_ll(self, lat, lon, e, nn, scale=1.0):   # the organic warp as lat/lon (same as the cell borders)
        km_deg = self.R * np.pi / 180.0
        return (np.clip(lat + scale * WARP_KM * nn / km_deg, -90, 90),
                lon + scale * WARP_KM * e / (km_deg * np.maximum(np.cos(np.radians(lat)), 0.05)))

    def _valleys(self, lat, lon, z, px_km, rough=None, micro=0.0):
        """Cut river valleys into heights (+ micro-relief) at (already warped) points: min(z, the lowest
        valley surface reaching each point), fading in as a valley grows from 1 to 2 px wide. Returns (z, channel,
        drawn channel)."""
        xyz = latlon_to_xyz(np.clip(lat, -90, 90), lon).reshape(-1, 3)
        V, channel, draw, base = self.river_net.valleys(xyz, px_km, rough, with_base=True)
        return self._roughen(z, micro, V, base), channel, draw

    def surface_at(self, lat, lon, px_km: float = MIN_KM):
        """Detailed height and water at points, as the tiles draw them: sea by the coast rule, lakes from warped cells."""
        lat = np.atleast_1d(np.asarray(lat, dtype=np.float64))
        lon = np.atleast_1d(np.asarray(lon, dtype=np.float64))
        z, sea, xyz, channel = self._z_sea(lat, lon, px_km)
        idx = self._warped_cells(lat, lon, *self._warp_offsets(xyz, 4 * px_km))
        lake = np.asarray(self.w.arrays["lake"])[idx].astype(bool) if "lake" in self.w.arrays else np.zeros(len(lat), bool)
        if not self.regions.empty:                             # inside refined areas: their own lakes (in channel)
            lake &= self.regions.weight(xyz) <= 0.5
        return z, coast_water(z, sea) | lake | channel

    def z_at(self, lat, lon, px_km: float = MIN_KM) -> np.ndarray:
        lat = np.atleast_1d(np.asarray(lat, dtype=np.float64))
        lon = np.atleast_1d(np.asarray(lon, dtype=np.float64))
        return self._z_sea(lat, lon, px_km)[0]

    def mesh(self, z, x, y):
        """(heights m, water) at a 3D patch's MESH_N × MESH_N vertices: rows north → south, columns west → east, tile
        corners included, detail only down to the vertex spacing. True heights (the sea floor below 0 m); water = sea
        by the coast rule (the client draws its surface at 0 m); lakes flat at their level."""
        span = 180.0 / 2 ** z
        f = np.arange(MESH_N) / (MESH_N - 1)
        LA, LO = np.meshgrid(90.0 - (y + f) * span, -180.0 + (x + f) * span, indexing="ij")
        h, sea, _, _ = self._z_sea(LA.ravel(), LO.ravel(), span / (MESH_N - 1) * np.pi / 180.0 * self.R)
        return h.reshape(MESH_N, MESH_N), coast_water(h, sea).reshape(MESH_N, MESH_N)

    def relief_at(self, lat, lon, px_km: float) -> np.ndarray:
        """World surface + procedural relief down to px_km, without river valleys (regional refinement's start)."""
        lat = np.atleast_1d(np.asarray(lat, dtype=np.float64))
        lon = np.atleast_1d(np.asarray(lon, dtype=np.float64))
        xyz = latlon_to_xyz(np.clip(lat, -90, 90), lon).reshape(-1, 3)
        amp, ridge, _ = self._amp(xyz)
        return self.bilinear(self.raster("elevation"), lat, lon) + self._detail(xyz, px_km, amp, ridge, lat, lon)

    def fields(self, z, x, y, border=0, need=("z",), valleys=True) -> dict:
        n = TILE + 2 * border
        j = np.arange(n) - border
        la, _ = tile_lat_lon(z, x, y, j)
        _, lo = tile_lat_lon(z, x, y, j)
        LA, LO = np.meshgrid(la, lo, indexing="ij")
        px_km = 180.0 / 2 ** z / TILE * np.pi / 180.0 * self.R
        ext = 4 * max(0, -(-(border - 4) // 4))               # wide borders: the lattice reaches past them too
        C = np.arange(-4 - ext, TILE + 5 + ext, 4)           # coarse lattice every 4 px, aligned across tiles
        m = len(C)
        cla, _ = tile_lat_lon(z, x, y, C)
        _, clo = tile_lat_lon(z, x, y, C)
        CLA, CLO = np.meshgrid(cla, clo, indexing="ij")
        cxyz = latlon_to_xyz(np.clip(CLA, -90, 90), CLO).reshape(-1, 3)
        rc = (j + 4 + ext) / 4.0
        RR, CC = np.meshgrid(rc, rc, indexing="ij")

        def up(a):
            return map_coordinates(np.asarray(a, dtype=np.float64).reshape(m, m), [RR, CC], order=1, mode="nearest")

        out = {"lat": LA, "lon": LO, "px_km": px_km}
        if "z" in need or "water" in need:
            amp, ridge, sea = self._amp(cxyz)
            amp, ridge, out["sea"] = up(amp), up(ridge), up(sea)
            out["amp"] = amp
            xyz = latlon_to_xyz(np.clip(LA, -90, 90), LO).reshape(-1, 3)
            D = self._detail(xyz, px_km, amp.ravel(), ridge.ravel(), LA.ravel(), LO.ravel()).reshape(n, n)
            out["base"] = self.bilinear(self.surface(), LA, LO)
            out["z"] = out["base"] + D
            micro = self._micro(xyz, px_km, amp.ravel()).reshape(n, n)
            out["river"] = out["river_draw"] = np.zeros((n, n), bool)
            if valleys:   # valley surfaces on the aligned 4-px lattice (seams agree), channels per pixel
                e, nn = self._warp_offsets(cxyz, 4 * px_km)
                cla_w, clo_w = self._warp_ll(CLA.ravel(), CLO.ravel(), e, nn)
                base, wall_up, rf = self.river_net.parts(latlon_to_xyz(cla_w, clo_w).reshape(-1, 3), px_km, 4 * px_km)
                base = up(np.where(np.isfinite(base), base, 1e9))
                V = base + np.maximum(0.0, up(wall_up) + up(rf) * D)
                wla, wlo = self._warp_ll(LA, LO, up(e), up(nn))
                ch, dr, lev = self.river_net.channel(latlon_to_xyz(wla, wlo).reshape(-1, 3), px_km)
                ch, dr = ch.reshape(n, n), dr.reshape(n, n)
                V = np.where(ch, lev.reshape(n, n), V)
                out["z"], out["river"], out["river_draw"] = self._roughen(out["z"], micro, V, base), ch, dr
            else:
                out["z"] = out["z"] + micro
            e, nn = (up(v) for v in self._warp_offsets(cxyz, 4 * px_km))
            out["z"] = self._flat_lakes(out["z"], LA, LO, e, nn)
            if not self.regions.empty and self.regions.touches(z, x, y):   # refined areas: their own heights, rivers,
                e, nn = self._region_warp(cxyz, px_km)                       # lakes and sea, blended in at the edge
                cla_w, clo_w = self._warp_ll(CLA.ravel(), CLO.ravel(), e, nn, scale=REGION_WARP)
                near = binary_dilation((self.regions.weight(cxyz) > 0).reshape(m, m), iterations=1).ravel()
                zr, w, chr_, drr, lk, sear = self._refined(
                    xyz, LA.ravel(), LO.ravel(), px_km, D.ravel(), (up(e).ravel(), up(nn).ravel()),
                    lattice=(latlon_to_xyz(cla_w, clo_w).reshape(-1, 3), up, near))
                if w.any():
                    w, inside = w.reshape(n, n), (w > 0.5).reshape(n, n)
                    out["z"] = (1 - w) * out["z"] + w * zr.reshape(n, n)
                    world = w < 1                                   # the world's rivers meet the refined ones in the band
                    out["river"] = (out["river"] & world) | chr_.reshape(n, n)
                    out["river_draw"] = (out["river_draw"] & world) | drr.reshape(n, n)
                    out["sea"] = np.where(inside, sear.reshape(n, n), out["sea"])
                    out["_region"] = (inside, lk.reshape(n, n))
            out["water"] = coast_water(out["z"], out["sea"])
        if "cat" in need:
            e, nn = (up(v) for v in self._warp_offsets(cxyz, 4 * px_km))
            idx = self._warped_cells(LA, LO, e, nn).reshape(n, n)
            for k in CAT_KEYS:
                if k in self.w.arrays:
                    out[k] = np.asarray(self.w.arrays[k])[idx]
            rs = self.regions
            if rs.touches(z, x, y):   # refined areas: categories of the fine cells (organic borders at their own scale)
                xyz = latlon_to_xyz(np.clip(LA, -90, 90), LO).reshape(-1, 3)
                inside = (rs.weight(xyz) > 0.5).reshape(n, n)
                if inside.any():
                    re_, rn_ = (up(v) for v in self._region_warp(cxyz, px_km))
                    wla, wlo = self._warp_ll(LA, LO, re_, rn_, scale=REGION_WARP)
                    fi = rs.nearest(latlon_to_xyz(wla, wlo).reshape(-1, 3))[0].reshape(n, n)
                    for k in CAT_KEYS:
                        if k in out and k in rs.arrays:
                            out[k] = np.where(inside, np.asarray(rs.arrays[k])[fi], out[k])
                    if "_region" in out:   # lakes: wherever the ground lies below a lake's surface (organic shores)
                        out["lake"] = np.where(out["_region"][0], out["_region"][1], out["lake"])
        return out

    # --- images -----------------------------------------------------------------------------------------------
    def _hillshade(self, z, lat, px_km, zoom):
        exag = float(np.clip(2.5 - 0.15 * (zoom - 5), 1.2, 2.5))   # gentler than the base's 4×: detail, not grain
        dy = px_km * 1000.0
        dx = dy * np.maximum(np.cos(np.radians(lat)), 0.01)
        gy, gx = np.gradient(z)
        dzdx, dzdn = gx / dx * exag, -gy / dy * exag
        a, b = np.radians(315.0), np.radians(45.0)
        L = (np.sin(a) * np.cos(b), np.cos(a) * np.cos(b), np.sin(b))
        hs = np.clip((-dzdx * L[0] - dzdn * L[1] + L[2]) / np.sqrt(dzdx**2 + dzdn**2 + 1.0), 0.0, 1.0)
        return hs[1:-1, 1:-1]

    def _warped(self, xyz, px_km, top_km, gain, seed, s):
        """Domain-warped octaves (Inigo Quilez, "Domain warping", 2002): noise at p + s·q(p), q three octave sums;
        s in units of top_km. Swirled, flowing patch shapes instead of round blobs. The warp itself is smooth: its
        octaves stop at top_km / 8."""
        q = np.stack([self._octaves(xyz, px_km, top_km, top_km / 8, 0.5, seed + 11 + i) for i in range(3)], 1)
        p = xyz + s * q * (top_km / self.R)
        return self._octaves(p / np.linalg.norm(p, axis=1, keepdims=True), px_km, top_km, MIN_KM, gain, seed)

    def _cover(self) -> tuple[np.ndarray, dict]:
        if getattr(self, "_cover_cache", None) is None:
            g = list(self.w.legends.get("ground", ()))
            self._cover_cache = (veg_cover(self.w.legends.get("holdridge", ())),
                                 {g.index(k): c for k, c in GROUND_COVER.items() if k in g})
        return self._cover_cache

    def _veg(self, f, B) -> dict:
        """Land colour of a tile from its fields with a border of B px (B ≥ LAND_BORDER reaches everything it uses,
        so tiles meet without seams; f from fields(…, ("z", "water", "cat"))): terrain shading (steep darker and drier, crests lighter and drier, green along
        water) and a vegetation mosaic — woods vs open ground, wooded share by Holdridge zone, woods in hollows and
        along water, not on crests or steep ground. Returns wood (0–1), bright (×) and tint (+ drier, − lusher) for
        the tile's own pixels."""
        from scipy.ndimage import distance_transform_edt
        from scipy.special import ndtri
        px, zz = f["px_km"], np.asarray(f["z"], dtype=np.float64)
        R = min(B, LAND_BORDER)
        c = lambda a: a[B:-B, B:-B]
        gy, gx = np.gradient(zz, px * 1000.0)
        steep = c(np.tanh(np.hypot(gx, gy) / STEEP))
        ridge = c(np.tanh((zz - f["base"]) / np.maximum(CREST_M, 0.5 * f["amp"])))   # + crest, − hollow or valley
        ridge = ridge * _resolved(TOP_KM, px)
        wet = np.asarray(f["water"], bool) | (np.asarray(f["lake"]) > 0) | np.asarray(f["river_draw"], bool)
        d = np.minimum(distance_transform_edt(~wet), R) if wet.any() else np.full(zz.shape, float(R))
        near = c(np.where(d < R, np.exp(-d / min(RIPARIAN_KM / px, R / 3.0)), 0.0))
        cover_z, cover_g = self._cover()
        zone = c(np.asarray(f["holdridge"]))
        cover = cover_z[np.clip(zone, 0, len(cover_z) - 1)] if len(cover_z) else np.zeros(zone.shape)
        ground = c(np.asarray(f["ground"]))
        for code, v in cover_g.items():
            cover = np.where(ground == code, v, cover)
        lat, lon = c(f["lat"]), c(f["lon"])
        xyz = latlon_to_xyz(np.clip(lat, -90, 90), lon).reshape(-1, 3)
        m = self._warped(xyz, px, VEG_KM, 0.6, self.seed + 7601, VEG_WARP).reshape(lat.shape)
        m = m + 0.25 * ridge - 0.5 * near + 0.3 * steep
        edge = max(0.04, px / 8.0)                           # patch edges ≈ a pixel soft at any zoom
        wood = np.where(cover > 0, 1.0 / (1.0 + np.exp(np.clip((m - VEG_SD * ndtri(np.clip(cover, 1e-3, 1 - 1e-3))) / edge, -60, 60))), 0.0)
        wood = cover + _resolved(VEG_KM, px) * (wood - cover)   # far out: the zone's mean, not a few-pixel speckle
        n1 = self._octaves(xyz, px, TINT_KM, MIN_KM, 0.6, self.seed + 7501).reshape(lat.shape)
        n2 = self._octaves(xyz, px, TINT_KM, TINT_KM / 8, 0.6, self.seed + 7502).reshape(lat.shape)
        terrain_b = 1.0 - 0.10 * steep + 0.06 * ridge - 0.05 * near
        terrain_t = 0.6 * (14.0 * steep + 10.0 * ridge - 26.0 * near)
        wood_b, wood_t = 0.9 - 0.08 * cover, -6.0 - 12.0 * cover          # scrub grey-olive … dark forest green
        open_t = np.minimum(16.0, -wood_t * cover / np.maximum(1.0 - cover, 0.05))   # clearings: lighter, the mean kept
        bright = terrain_b * (1.0 + (wood_b - 1.0) * wood) + 0.3 * TINT[0] * n1
        tint = terrain_t + wood * wood_t + (1.0 - wood) * open_t + 0.3 * TINT[1] * n2
        return {"wood": wood, "bright": bright, "tint": tint}

    def rgb(self, layer, z, x, y) -> np.ndarray:
        if layer == "relief":
            B = LAND_BORDER
            f = self.fields(z, x, y, B, ("z", "water", "cat"))
            c = lambda a: a[B:-B, B:-B]
            o = lambda a: a[B - 1:1 - B, B - 1:1 - B]                 # the hillshade's one-pixel border
            hs = self._hillshade(o(f["z"]), o(f["lat"]), f["px_km"], z)
            v = self._veg(f, B)
            return RN.relief_rgb(c(f["z"]), hs, c(f["holdridge"]), c(f["ground"]), c(f["ice"]), c(f["lake"]) | c(f["river_draw"]), ~c(f["water"]),   # rivers: water
                                 vary=(v["bright"], v["tint"]), style=self.w.cells_meta.get("style"))
        if layer == "seafloor":
            f = self.fields(z, x, y, 1, ("z", "water", "cat"))
            c = lambda a: a[1:-1, 1:-1]
            hs = self._hillshade(f["z"], f["lat"], f["px_km"], z)
            seabed = c(f["seabed_type"]) if "seabed_type" in f else np.zeros(hs.shape, np.int64)
            return seafloor_rgb(c(f["z"]), hs, seabed, ~c(f["water"]))
        if layer == "elevation":
            return RN._colorize("elevation", self.fields(z, x, y, 0, ("z",))["z"])
        if layer in CONT:
            la, lo = tile_lat_lon(z, x, y, np.arange(TILE))
            LA, LO = np.meshgrid(la, lo, indexing="ij")
            return RN._colorize(CONT[layer], self.bilinear(self.raster(CONT[layer]), LA, LO))
        key, legend = CAT[layer]
        f = self.fields(z, x, y, 0, ("cat",))
        pal = RN.holdridge_palette() if key == "holdridge" else RN.category_palette(len(self.w.legends[legend]))
        return np.asarray(pal)[np.clip(f[key], 0, len(pal) - 1)].astype(np.uint8)

    def render(self, layer, z, x, y) -> bytes:
        if layer == "mesh":
            h, water = self.mesh(z, x, y)
            return h.astype("<f4").tobytes() + water.astype(np.uint8).tobytes()
        buf = io.BytesIO()
        if layer == "height":
            Image.fromarray(RN.encode(self.fields(z, x, y, 0, ("z",))["z"], 0.5, -12000.0)).save(buf, "PNG")
        else:
            Image.fromarray(np.ascontiguousarray(self.rgb(layer, z, x, y), dtype=np.uint8)).save(buf, "JPEG", quality=90)
        return buf.getvalue()

    def start_workers(self, n: int) -> None:
        """Fork n render processes for this source and its peers (call before other threads start): every world is
        loaded first, then shared copy-on-write."""
        for s in self.peers.values():
            s.tree()
            s.river_net
            s.raster("elevation")
            s.regions.warm()                                    # search trees in the parent: workers share them
        pool = RenderPool(self.peers, n) if n > 0 else None
        for s in self.peers.values():
            s.pool = pool

    def prerender(self, max_z: int = 11, layers=("relief", "mesh"), progress=None, threads: int = 4,
                  stop=lambda: False, max_tiles: int | None = None) -> int:
        """Render and save every tile of the refined areas from z5 to max_z (coarse first; saved tiles are skipped),
        leaving out the deeper zooms that would take the job count over max_tiles. progress(done, total).
        Returns the number rendered."""
        rs = self.regions
        if rs.empty:
            return 0
        lat, lon = np.asarray(rs.arrays["g_lat"]), np.asarray(rs.arrays["g_lon"])
        jobs = []
        for z in range(5, max_z + 1):
            span = 180.0 / 2 ** z
            xs = np.floor((lon + 180.0) / span).astype(np.int64)
            ys = np.floor((90.0 - lat) / span).astype(np.int64)
            cand = {(int(x + dx) % 2 ** (z + 1), int(y + dy)) for x, y in set(zip(xs.tolist(), ys.tolist()))
                    for dx in (-1, 0, 1) for dy in (-1, 0, 1)}
            level = [(layer, z, x, y) for x, y in sorted(cand) if 0 <= y < 2 ** z and rs.touches(z, x, y)
                     for layer in layers]
            if max_tiles is not None and jobs and len(jobs) + len(level) > max_tiles:
                break                                           # huge areas: deeper zooms render when viewed
            jobs += level
        done, rendered, lock, todo = 0, 0, threading.Lock(), iter(jobs)

        def run():                                              # daemon threads taking one job at a time: a stop
            nonlocal done, rendered                             # or the server's exit ends them at once
            while not stop() and self.regions is rs:
                with lock:
                    job = next(todo, None)
                if job is None:
                    return
                fresh = not self.cache_path(*job).exists()
                if fresh:
                    self.tile(*job, remember=False)
                with lock:
                    done += 1
                    rendered += fresh
                    if progress:
                        progress(done, len(jobs))

        ts = [threading.Thread(target=run, daemon=True) for _ in range(max(1, threads))]
        for t in ts:
            t.start()
        for t in ts:
            t.join()
        return rendered

    def export_tile(self, z, x, y) -> dict:
        """A tile's export layers (export.tile_data), on a render worker when there are some."""
        import export
        with self.lock:
            rs = self.regions
        if self.pool is not None:
            try:
                return self.pool.render(self.name, EXPORT_JOB, z, x, y, rs.fingerprint)
            except RuntimeError:
                pass
        return export.tile_data(self, z, x, y)

    def _own_areas(self) -> list:
        """Refined areas of this era whose content differs from the base's (new, or a different key)."""
        mine, theirs = self.regions.area_keys, self.share.regions.area_keys
        return [h for h, k in mine.items() if theirs.get(h) != k]

    def shared(self, z, x, y) -> bool:
        """An era's tile that is exactly the base's: away from the influence mask (+ CARRY_KM) and from refined
        areas the era changed; the base renders and saves it once for every era."""
        if self.share is None:
            return False
        if self._mask_tree is not None and refine.touches_trees([self._mask_tree], self._spacing_km, self.R, z, x, y,
                                                                CARRY_KM):
            return False
        own = self._own_areas()
        return not (own and self.regions.touches(z, x, y, own, CARRY_KM))

    def cache_path(self, layer, z, x, y) -> Path:
        if self.shared(z, x, y):
            return self.share.cache_path(layer, z, x, y)
        with self.lock:
            rs, tag = self.regions, self.tag
        return self._path(rs, tag, layer, z, x, y)

    def _path(self, rs, tag, layer, z, x, y) -> Path:
        return self.cache_root / (tag if rs.touches(z, x, y) else self.base_tag) / layer / str(z) / str(x) / \
            f"{y}.{EXT.get(layer, 'jpg')}"

    def _render_for(self, rs, layer, z, x, y) -> bytes:
        if self.pool is not None:
            try:
                return self.pool.render(self.name, layer, z, x, y, rs.fingerprint)
            except RuntimeError:
                pass                                            # a worker failed or lags behind: render here
        return self.render(layer, z, x, y)

    def tile(self, layer, z, x, y, remember: bool = True, gate=None) -> bytes:
        """gate: a semaphore of render slots (saved tiles never wait for one); none free → Busy."""
        if (layer not in LAYERS and layer not in EXT) or not tile_ok(z, x, y):
            raise KeyError(f"no tile {layer}/{z}/{x}/{y}")
        if self.shared(z, x, y):
            return self.share.tile(layer, z, x, y, remember, gate)
        data = b""
        for _ in range(3):                                      # a region set swapped mid-render: render again
            with self.lock:
                rs, tag = self.regions, self.tag
                key = (tag, layer, z, x, y)
                if key in self.mem:
                    self.mem.move_to_end(key)
                    return self.mem[key]
            f = self._path(rs, tag, layer, z, x, y)
            if f.exists() and (layer != "mesh" or f.stat().st_size == MESH_BYTES):   # an old mesh format: again
                data = f.read_bytes()
                if self.cap is not None:
                    self.cap.used(f)
            else:
                if gate is not None and not gate.acquire(blocking=False):
                    raise Busy(f"{layer} {z}/{x}/{y}")
                try:
                    data = self._render_for(rs, layer, z, x, y)
                finally:
                    if gate is not None:
                        gate.release()
                with self.lock:
                    swapped = self.tag != tag
                if swapped:
                    continue                                    # may mix both sets: never save it
                f.parent.mkdir(parents=True, exist_ok=True)
                fd, tmp = tempfile.mkstemp(dir=f.parent, suffix=".part")
                with os.fdopen(fd, "wb") as fh:
                    fh.write(data)
                os.replace(tmp, f)
                if self.cap is not None:
                    self.cap.saved(f, len(data))
            if remember:
                with self.lock:
                    if self.tag == tag:
                        self.mem[key] = data
                        while len(self.mem) > MEM_TILES:
                            self.mem.popitem(last=False)
            return data
        return data


def reload_all(sources: list) -> None:
    """A region build finished: every source reloads, then the shared pool forks once."""
    for s in sources:
        s.reload_regions(refork=False)
    if sources:
        sources[0].refork()


def join(sources: list) -> None:
    """Sources that share one render pool (the eras of a server): each knows the others by name."""
    peers = {s.name: s for s in sources}
    for s in sources:
        s.peers = peers