worldmap-viewer: initial public history

1 files changed, 323 insertions, 0 deletions

diff --git a/seafloor.py b/seafloor.py
new file mode 100644
index 0000000..86791e5
--- /dev/null
+++ b/seafloor.py
@@ -0,0 +1,323 @@
+"""Sea-floor detail for refined areas:
+turbidity canyons and their fans, sediment ponds, fine volcanic fields on the sunken plateaus, vents (data) and the
+refined sea-floor fields. Only sea cells inside the area change, and never above SEA_TOP_M; land and the halo keep
+their heights. Deterministic: the same cells and world values give the same floor. Generated terrain (`idea`)."""
+from __future__ import annotations
+
+import json
+from pathlib import Path
+
+import numpy as np
+from scipy.spatial import cKDTree
+
+import h3par
+import worldgen_path  # noqa: F401  (mapgen on sys.path)
+from mapgen import plateaus as PL, seabed as SB
+from mapgen.graph import accumulate, components, priority_flood, receiver_levels, smooth_km, steepest_receivers
+from mapgen.noise import name_seed
+from mapgen.sphere import east_north, gc_dist_km, great_circle_point, latlon_to_xyz
+from mapgen.zones import smootherstep
+
+SEA_TOP_M = -5.0               # the pass never lifts a sea cell above this (no new land)
+WALL_M = 1.0e7                 # land in the sea-floor routing
+CANYON = {"k": 0.1, "m": 0.5, "n": 1.0, "min_area_km2": 500.0, "min_slope": 2.0, "cap_m": 1500.0}   # slope: m/km
+FAN = {"flat_slope": 1.0, "max_m": 100.0, "km": 15.0}
+POND = {"min_depth_m": 50.0, "fill": 0.3}
+VENT_TYPES = ["black smoker", "white smoker", "diffuse", "cold seep"]
+VENT_MINERALS = ["copper-iron sulfides", "zinc sulfides and barite", "iron-manganese oxides", "methane carbonates"]
+VENT = {"rate": 0.03, "seep_rate": 0.002, "min_depth_m": 200.0, "ridge_km": 60.0, "reach_km": 3.0}
+TEMP_C = (np.array([300.0, 100.0, 10.0, 0.0]), np.array([100.0, 200.0, 90.0, 5.0]))   # per type: from, span
+
+
+def canyons(g, z, sea, halo, P=CANYON, F=FAN):
+    """(heights, canyon depth, fan thickness). Turbidity flows follow the steepest descent over the sea floor (land is
+    a wall); where their contributing sea area and slope are high they cut (stream-power form, like rivers, capped
+    at cap_m); the material they carry settles as a fan where they first reach flat floor (slope < flat_slope),
+    spread over ≈ km and at most max_m thick. Canyons only lower, fans only raise; the halo stays."""
+    z = np.asarray(z, dtype=np.float64)
+    sea, halo = np.asarray(sea, bool), np.asarray(halo, bool)
+    active = sea & ~halo
+    recv, slope, _ = steepest_receivers(g, np.where(sea, z, WALL_M))
+    levels = receiver_levels(recv)
+    area = accumulate(recv, levels, np.where(sea, g.area_km2, 0.0))
+    cut = np.where(active & (area >= P["min_area_km2"]) & (slope >= P["min_slope"]),
+                   np.minimum(P["cap_m"], P["k"] * area ** P["m"] * slope ** P["n"]), 0.0)
+    vol = accumulate(recv, levels, cut * g.area_km2 / 1000.0)                  # km³ carried down
+    flat = active & (slope < F["flat_slope"])
+    ar = np.arange(g.n)
+    src = active & ~flat & (recv != ar) & (vol > 0)
+    entry = np.zeros(g.n, bool)
+    entry[recv[src]] = True
+    entry &= flat
+    dep = np.where(entry, vol * 1000.0 / g.area_km2, 0.0)                       # m over the entry cell
+    fan = np.where(active & (cut == 0), np.clip(smooth_km(g, dep, F["km"]), 0.0, F["max_m"]), 0.0)
+    zc = z - cut
+    zn = np.where(fan > 0, np.maximum(zc, np.minimum(zc + fan, SEA_TOP_M)), zc)
+    return zn, cut, zn - zc
+
+
+def ponds(g, z, sea, halo, sediment_m, P=POND):
+    """(heights, fill thickness): closed sea-floor hollows deeper than min_depth_m below their spill level fill
+    with sediment to a flat floor at min(spill level, lowest point + fill × the mean sediment there). Never above
+    the spill level, never lower."""
+    z = np.asarray(z, dtype=np.float64)
+    sea, halo = np.asarray(sea, bool), np.asarray(halo, bool)
+    out = z.copy()
+    if not (sea & ~halo).any():
+        return out, np.zeros(g.n)
+    sinks = sea & halo
+    if not sinks.any():                                       # sea closed inside the area: its deepest cell drains
+        sinks = np.zeros(g.n, bool)
+        s = np.flatnonzero(sea)
+        sinks[s[np.argmin(z[s])]] = True
+    zf = priority_flood(g, np.where(sea, z, WALL_M), sinks, eps=0.0)
+    hollow = sea & ~halo & (zf - z > 0.5)
+    if not hollow.any():
+        return out, np.zeros(g.n)
+    lab = components(g, hollow)
+    idx = np.flatnonzero(hollow)
+    L = lab[idx]
+    n = int(L.max()) + 1
+    low = np.full(n, np.inf)
+    np.minimum.at(low, L, z[idx])
+    spill = np.full(n, -np.inf)
+    np.maximum.at(spill, L, zf[idx])
+    cnt = np.bincount(L, minlength=n)
+    sup = np.bincount(L, weights=np.asarray(sediment_m, dtype=np.float64)[idx], minlength=n) / np.maximum(cnt, 1)
+    level = np.minimum(spill, low + P["fill"] * sup)
+    deep = (spill - low) >= P["min_depth_m"]
+    lv = np.where(deep[L], np.minimum(level[L], SEA_TOP_M), -np.inf)
+    out[idx] = np.maximum(z[idx], lv)
+    return out, out - z
+
+
+def _u(ids, salt: int) -> np.ndarray:
+    """Uniform [0, 1) per H3 cell id (splitmix64): a cell draws the same number in every area and every era."""
+    x = np.asarray(ids, dtype=np.uint64) + np.uint64((int(salt) * 0x9E3779B97F4A7C15) & 0xFFFFFFFFFFFFFFFF)
+    x = x + np.uint64(0x9E3779B97F4A7C15)
+    x = (x ^ (x >> np.uint64(30))) * np.uint64(0xBF58476D1CE4E5B9)
+    x = (x ^ (x >> np.uint64(27))) * np.uint64(0x94D049BB133111EB)
+    x = x ^ (x >> np.uint64(31))
+    return (x >> np.uint64(11)).astype(np.float64) / float(1 << 53)
+
+
+def _around(rng, lat, lon, reach_km, radius_km):
+    c = latlon_to_xyz(lat, lon)
+    e, n = east_north(c[None])
+    az = rng.uniform(0.0, 2.0 * np.pi)
+    return great_circle_point(c, np.cos(az) * n[0] + np.sin(az) * e[0], rng.uniform(0.0, reach_km), radius_km)
+
+
+def small_features(p: dict, feat: dict, seed: int, radius_km: float) -> list:
+    """A plateau's fine volcanic features, deterministic per name and independent of any area: 2–6 cones (2–15 km
+    wide, 0.2–2 km tall) around each world-scale cone, 1–2 pit calderas in each world caldera, and a fissure ridge
+    from the hotspot to each caldera (the hotspot's feeding track)."""
+    rng = np.random.default_rng(name_seed(seed, p["name"]) + 11)
+    vent = float(p.get("vent", 1.0))
+    out = []
+    for c in feat["cones"]:
+        for _ in range(int(rng.integers(2, 7))):
+            out.append({"kind": "cone", "at": _around(rng, c["lat"], c["lon"], c["radius_km"], radius_km),
+                        "radius_km": float(rng.uniform(1.0, 7.5)), "height_m": float(rng.uniform(200.0, 2000.0)) * vent})
+    for c in feat["calderas"]:
+        for _ in range(int(rng.integers(1, 3))):
+            out.append({"kind": "pit", "at": _around(rng, c["lat"], c["lon"], 0.5 * c["radius_km"], radius_km),
+                        "radius_km": float(rng.uniform(2.0, 6.0)), "height_m": 150.0 * vent, "floor_m": -250.0 * vent})
+    hot = latlon_to_xyz(*feat["hotspot"])
+    for c in feat["calderas"]:
+        out.append({"kind": "fissure", "at": hot, "to": latlon_to_xyz(c["lat"], c["lon"]), "radius_km": 1.5,
+                    "height_m": float(rng.uniform(100.0, 400.0)) * vent})
+    return out
+
+
+def _segment_km(xyz, a, b, radius_km, step_km=1.0):
+    n = max(2, int(gc_dist_km(a, b, radius_km) / step_km) + 1)
+    t = np.linspace(0.0, 1.0, n)[:, None]
+    pts = a * (1.0 - t) + b * t
+    pts /= np.linalg.norm(pts, axis=1, keepdims=True)
+    chord, _ = cKDTree(pts).query(xyz, workers=max(1, h3par.workers()))
+    return 2.0 * np.arcsin(np.clip(chord / 2.0, 0.0, 1.0)) * radius_km
+
+
+def volcanic(g, sea, halo, plateau_id, plateaus: list, seed: int):
+    """(height to add, nearness 0–1 to a feature) from the fine volcanic features of the plateaus under the area's
+    sea cells (vent = 0: none)."""
+    add, near = np.zeros(g.n), np.zeros(g.n)
+    for k, p in enumerate(plateaus):
+        on = np.flatnonzero((np.asarray(plateau_id) == k) & sea & ~halo)
+        if not len(on) or float(p.get("vent", 1.0)) <= 0:
+            continue
+        xyz = g.xyz[on]
+        for f in small_features(p, PL.features(p, seed, g.radius_km), seed, g.radius_km):
+            r = f["radius_km"]
+            if f["kind"] == "fissure":
+                d = _segment_km(xyz, f["at"], f["to"], g.radius_km)
+                h = f["height_m"] * np.exp(-(d / r) ** 2)
+            else:
+                d = gc_dist_km(xyz, f["at"], g.radius_km)
+                if f["kind"] == "cone":
+                    h = f["height_m"] * np.clip(1.0 - d / r, 0.0, 1.0) ** 1.5
+                else:
+                    h = f["height_m"] * np.exp(-((d - r) / (0.25 * r)) ** 2) + f["floor_m"] * (1.0 - smootherstep(d / r))
+            add[on] += h
+            near[on] = np.maximum(near[on], np.exp(-(d / (r + 3.0)) ** 2))
+    return add, near
+
+
+def clamp_plateaus(z, sea, halo, plateau_id, plateaus: list):
+    """Hidden plateaus keep every sea cell at least HIDDEN_MAX_M (−550 m) deep."""
+    z = np.asarray(z, dtype=np.float64).copy()
+    for k, p in enumerate(plateaus):
+        if not p.get("islands", False):
+            m = (np.asarray(plateau_id) == k) & sea & ~halo
+            z[m] = np.minimum(z[m], PL.HIDDEN_MAX_M)
+    return z
+
+
+def vents(g, z, sea, halo, potential, near, d_div_km, sediment_m, bottom_c, seed: int, V=VENT):
+    """(per-cell vent strength 0–1, vents). A sea cell deeper than min_depth_m holds a vent with probability
+    rate × potential × (0.5 + nearness to a volcanic feature); on a ridge axis (d_div_km < ridge_km) vents are black
+    smokers, near a feature black (60 %) or white smokers, elsewhere diffuse; cold seeps sit in thick sediment where
+    the potential is low. Temperature and flow class (0–2) per vent; the strength spreads ≈ reach_km around each."""
+    depth = -np.asarray(z, dtype=np.float64)
+    cand = sea & ~halo & (depth > V["min_depth_m"])
+    pot = np.clip(np.asarray(potential, dtype=np.float64), 0.0, 1.0)
+    near = np.asarray(near, dtype=np.float64)
+    u = [_u(g.ids, seed * 8 + s) for s in range(5)]
+    hot = cand & (u[0] < V["rate"] * pot * (0.5 + near))
+    seep = cand & ~hot & (pot < 0.2) & (np.asarray(sediment_m) > 1000.0) & (u[1] < V["seep_rate"])
+    i = np.flatnonzero(hot | seep)
+    ridge = np.asarray(d_div_km)[i] < V["ridge_km"]
+    close = near[i] > 0.3
+    t = np.where(seep[i], 3, np.where(ridge | (close & (u[2][i] < 0.6)), 0, np.where(close, 1, 2))).astype(np.int8)
+    lo, span = TEMP_C
+    temp = np.where(t == 3, np.asarray(bottom_c, dtype=np.float64)[i], 0.0) + lo[t] + span[t] * u[3][i]
+    flow = np.minimum((3.0 * u[4][i] * np.sqrt(pot[i])).astype(np.int8), 2).astype(np.int8)
+    strength = np.zeros(g.n)
+    if len(i):
+        reach = V["reach_km"] / g.radius_km
+        d, k = cKDTree(g.xyz[i]).query(g.xyz, distance_upper_bound=3.0 * reach, workers=max(1, h3par.workers()))
+        ok = np.isfinite(d) & sea
+        w = 0.5 + 0.25 * flow
+        strength[ok] = np.clip(w[k[ok]] * np.exp(-(d[ok] / reach) ** 2), 0.0, 1.0)
+    out = {"cell": g.ids[i].astype(np.uint64), "lat": g.lat[i].astype(np.float64), "lon": g.lon[i].astype(np.float64),
+           "type": t, "temp_c": temp.astype(np.float32), "flow": flow, "mineral": t.copy()}
+    return strength, out
+
+
+def fields(pv: dict, sea, halo, cut, fan, fill, volc, strength) -> dict:
+    """The refined sea-floor fields (0 / none on land): the world cell's values; sea over a world land cell (a refined
+    coast) starts as terrigenous sediment; canyons scour (their fans and pond fill add sediment) and read as
+    terrigenous; fine volcanic features read as volcanic floor; vent fields over the vents (sulfides, no sediment,
+    warmer water)."""
+    sea = np.asarray(sea, bool)
+    t = np.asarray(pv["seabed_type"]).astype(np.int8).copy()
+    m = np.asarray(pv["seabed_mineral"]).astype(np.int8).copy()
+    sed = np.asarray(pv["sediment_m"], dtype=np.float64).copy()
+    bt = np.asarray(pv["bottom_temp_c"], dtype=np.float64).copy()
+    orphan = sea & (t == SB.SB_NONE)
+    t[orphan] = SB.SB_TERRIGENOUS
+    sed[orphan] = 300.0
+    bt[orphan] = np.maximum(np.asarray(pv["T_mean"], dtype=np.float64)[orphan], -1.8)
+    cut, fan, fill = (np.asarray(a, dtype=np.float64) for a in (cut, fan, fill))
+    volc, strength = np.asarray(volc, dtype=np.float64), np.asarray(strength, dtype=np.float64)
+    sed = np.where(cut > 0, 0.3 * sed, sed) + fan + fill
+    t[(cut > 50.0) | (fan > 1.0)] = SB.SB_TERRIGENOUS
+    t[np.abs(volc) > 50.0] = SB.SB_VOLCANIC
+    vf = strength > 0.5
+    t[vf] = SB.SB_VENTS
+    m[vf] = SB.MI_SULFIDES
+    sed[vf] = 0.0
+    bt = bt + 20.0 * strength
+    land = ~sea
+    t[land], m[land] = SB.SB_NONE, SB.MI_NONE
+    z32 = lambda a: np.where(land, 0.0, a).astype(np.float32)
+    return {"seabed_type": t, "seabed_mineral": m, "bottom_temp_c": z32(bt), "sediment_m": z32(sed),
+            "vent": z32(strength), "canyon_m": z32(cut), "fan_m": z32(fan)}
+
+
+PV_KEYS = ("vent_potential", "seabed_type", "seabed_mineral", "bottom_temp_c", "sediment_m", "d_div_km", "T_mean")
+
+
+def run(g, z, sea, halo, parent, world_a, plateaus: list, seed: int):
+    """The whole pass: (heights, refined sea-floor fields incl. plateau_id, vents)."""
+    sea, halo = np.asarray(sea, bool), np.asarray(halo, bool)
+    pv = {k: np.asarray(world_a[k])[parent] for k in PV_KEYS}
+    plateau_id = (PL.cell_ids(g.xyz, plateaus, seed, g.radius_km) if plateaus
+                  else np.full(g.n, -1, np.int16))
+    z0 = np.asarray(z, dtype=np.float64)
+    zc, cut, fan = canyons(g, z0, sea, halo)
+    zp, fill = ponds(g, zc, sea, halo, pv["sediment_m"])
+    add, near = volcanic(g, sea, halo, plateau_id, plateaus, seed)
+    zv = np.where(sea & ~halo, np.minimum(zp + add, np.maximum(zp, SEA_TOP_M)), zp)
+    zn = clamp_plateaus(zv, sea, halo, plateau_id, plateaus)
+    strength, vt = vents(g, zn, sea, halo, pv["vent_potential"], near, pv["d_div_km"], pv["sediment_m"],
+                         pv["bottom_temp_c"], seed)
+    f = fields(pv, sea, halo, cut, fan, fill, zv - zp, strength)
+    f["plateau_id"] = plateau_id.astype(np.int16)
+    return zn, f, vt
+
+
+VENT_RGB = [(230, 40, 30), (245, 245, 245), (250, 160, 40), (40, 220, 230)]   # by VENT_TYPES
+
+
+def preview(area_dir: Path, path: Path, px: int = 1400) -> Path:
+    """A shaded depth image of a refined area (its box, north up) with its vents as dots — black smokers red, white
+    smokers white, diffuse orange, cold seeps cyan: a first look at a sea floor before the viewer shows it."""
+    from PIL import Image, ImageDraw
+    area_dir, path = Path(area_dir), Path(path)
+    with np.load(area_dir / "cells.npz") as a:
+        inner = ~a["halo"].astype(bool)
+        lat, lon = a["g_lat"][inner], a["g_lon"][inner]
+        z, sea = a["elevation_eroded_m"][inner].astype(np.float64), a["ocean"][inner].astype(bool)
+    c0 = float(np.degrees(np.arctan2(np.sin(np.radians(lon)).mean(), np.cos(np.radians(lon)).mean())))
+    x = (lon - c0 + 180.0) % 360.0 - 180.0
+    k = float(np.cos(np.radians(lat.mean())))
+    W = int(px)
+    H = max(2, int(round(W * (lat.max() - lat.min()) / max((x.max() - x.min()) * k, 1e-9))))
+    GX, GY = np.meshgrid(np.linspace(x.min(), x.max(), W), np.linspace(lat.max(), lat.min(), H))
+    pts = np.column_stack([x * k, lat])
+    tree = cKDTree(pts)
+    dist, j = tree.query(np.column_stack([GX.ravel() * k, GY.ravel()]))
+    spacing = float(np.median(tree.query(pts, k=2)[0][:, 1]))
+    Z, S = z[j].reshape(H, W), sea[j].reshape(H, W)
+    px_km = (lat.max() - lat.min()) * 111.2 / H
+    gy, gx = np.gradient(Z, px_km)
+    shade = np.clip(1.0 - 0.015 * (gx - gy), 0.45, 1.35)[..., None]
+    t = np.clip(-Z / 6000.0, 0.0, 1.0)[..., None]
+    rgb = np.where(S[..., None], (1.0 - t) * np.array([150.0, 210.0, 235.0]) + t * np.array([10.0, 30.0, 80.0]),
+                   np.array([125.0, 140.0, 95.0])) * shade
+    rgb[(dist > 2.0 * spacing).reshape(H, W)] = 40.0
+    im = Image.fromarray(np.clip(rgb, 0, 255).astype(np.uint8))
+    if (area_dir / "vents.npz").exists():
+        d = ImageDraw.Draw(im)
+        with np.load(area_dir / "vents.npz") as v:
+            vx = ((v["lon"] - c0 + 180.0) % 360.0 - 180.0 - x.min()) / max(x.max() - x.min(), 1e-9) * (W - 1)
+            vy = (lat.max() - v["lat"]) / max(lat.max() - lat.min(), 1e-9) * (H - 1)
+            for a_, b_, ty in zip(vx.tolist(), vy.tolist(), v["type"].tolist()):
+                d.ellipse([a_ - 3, b_ - 3, a_ + 3, b_ + 3], fill=VENT_RGB[int(ty)], outline=(0, 0, 0))
+    path.parent.mkdir(parents=True, exist_ok=True)
+    im.save(path, quality=90)
+    return path
+
+
+def preview_region(root: Path, res: int, region_id: str, regions_path: Path | None = None,
+                   regions_root: Path | None = None, out_dir: Path | None = None) -> Path:
+    """preview() of a region's refined area in every world that has it: <out_dir>/<region>-<era>.jpg."""
+    import refine as RF
+    root = Path(root)
+    regions = RF.load_regions(regions_path or root / "places" / "regions.json")
+    if not any(RF.region_key(r) == region_id for r in regions):
+        raise SystemExit(f"no region {region_id} in {regions_path or 'places/regions.json'}")
+    fine = res + RF.FINER
+    h = RF.region_areas(regions, RF.area_cells(regions, fine), fine)[region_id]["area"]
+    rroot = Path(regions_root or root / "out" / f"r{res}" / "regions")
+    out_dir = Path(out_dir or root / "previews" / f"r{res}" / "seafloor")
+    made = []
+    for s in sorted(p for p in rroot.glob("*-m*") if p.is_dir()):
+        if h and (s / h / "cells.npz").exists():
+            era = json.loads((s / "era.json").read_text())["era"] if (s / "era.json").exists() else s.name
+            made.append(preview(s / h, out_dir / f"{region_id}-{era}.jpg"))
+    if not made:
+        raise SystemExit(f"{region_id} is not refined at res {res}: run mapgen.py refine --res {res}")
+    return out_dir