worldmap-viewer: initial public history

1 files changed, 398 insertions, 0 deletions

diff --git a/export.py b/export.py
new file mode 100644
index 0000000..567048d
--- /dev/null
+++ b/export.py
@@ -0,0 +1,398 @@
+"""Terrain export for games: freeze a square of the map, exactly as the viewer shows it, into engine-neutral files. Generated terrain."""
+from __future__ import annotations
+
+import json
+import math
+import multiprocessing
+import os
+import shutil
+import tempfile
+from datetime import date
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+from scipy.ndimage import map_coordinates
+
+import worldgen_path  # noqa: F401  (mapgen on sys.path)
+import refine
+import rivers as RV
+import serve
+import tiles as T
+from mapgen import render as RN
+from mapgen.sphere import latlon_to_xyz
+
+FORMAT_VERSION = 1
+WATER = ["land", "sea", "lake", "river"]
+
+
+def samples_for(size_km: float, res_m: float) -> int:
+    """The smallest 2^n + 1 covering size_km at res_m (Unity, Godot Terrain3D and most tools take these sizes)."""
+    need = size_km * 1000.0 / res_m + 1
+    n = 1
+    while 2 ** n + 1 < need:
+        n += 1
+    return 2 ** n + 1
+
+
+def grid_latlon(lat0, lon0, n, res_m, radius_km, rows=None):
+    """Latitude/longitude of the n × n samples (or of rows a:b): azimuthal equidistant about (lat0, lon0), x east,
+    y north, row 0 north."""
+    c = (n - 1) / 2
+    x = (np.arange(n) - c) * res_m
+    r = np.arange(n) if rows is None else np.arange(*rows)
+    y = (c - r) * res_m
+    X, Y = np.meshgrid(x, y)
+    d = np.hypot(X, Y) / (radius_km * 1000.0)                  # angular distance from the centre
+    b = np.arctan2(X, Y)                                        # bearing
+    la0, lo0 = math.radians(lat0), math.radians(lon0)
+    la = np.arcsin(np.sin(la0) * np.cos(d) + np.cos(la0) * np.sin(d) * np.cos(b))
+    lo = lo0 + np.arctan2(np.sin(b) * np.sin(d) * np.cos(la0), np.cos(d) - np.sin(la0) * np.sin(la))
+    return np.degrees(la), (np.degrees(lo) + 180.0) % 360.0 - 180.0
+
+
+def to_local(lat0, lon0, lat, lon, radius_km):
+    """Local x east, y north (m) of points (the inverse of grid_latlon)."""
+    la0, lo0 = math.radians(lat0), math.radians(lon0)
+    la, lo = np.radians(lat), np.radians(lon)
+    cd = np.clip(np.sin(la0) * np.sin(la) + np.cos(la0) * np.cos(la) * np.cos(lo - lo0), -1, 1)
+    d = np.arccos(cd)
+    b = np.arctan2(np.sin(lo - lo0) * np.cos(la), np.cos(la0) * np.sin(la) - np.sin(la0) * np.cos(la) * np.cos(lo - lo0))
+    r = d * radius_km * 1000.0
+    return r * np.sin(b), r * np.cos(b)
+
+
+_SRC = None                                                     # the tile source (forked workers inherit it)
+
+
+def _tile(args):
+    return tile_data(_SRC, *args)
+
+
+def tile_data(src, z, x, y) -> dict:
+    """One tile's export layers (256 × 256): heights, water classes, biome, ground, landform, unshaded colours."""
+    f = src.fields(z, x, y, 1, ("z", "water", "cat"))
+    c = lambda a: np.asarray(a)[1:-1, 1:-1]
+    lake, river = c(f["lake"]).astype(bool), c(f["river"]).astype(bool)
+    flat = np.ones(lake.shape)                                  # unshaded colours: the export shades its own grid
+    rgb = RN.relief_rgb(c(f["z"]), flat, c(f["holdridge"]), c(f["ground"]), c(f["ice"]), lake | c(f["river_draw"]),
+                        ~c(f["water"]))
+    water = np.where(c(f["water"]), 1, np.where(lake, 2, np.where(river, 3, 0))).astype(np.uint8)
+    return {"z": c(f["z"]).astype(np.float32), "water": water, "biome": c(f["holdridge"]).astype(np.uint8),
+            "ground": c(f["ground"]).astype(np.uint8), "landform": c(f["landform"]).astype(np.uint8),
+            "rgb": np.asarray(rgb, dtype=np.uint8)}
+
+
+def _rivers(src, rs, lat0, lon0, half_m, res_m):
+    """River centrelines inside the square (local metres), world rivers where the world shows, refined ones inside."""
+    R = src.R
+    centre = latlon_to_xyz(lat0, lon0).reshape(3)
+    radius = half_m * math.sqrt(2) / 1000.0
+    nets = [(src.river_net, False)] + ([(rs.net, True)] if not rs.empty and rs.net is not None else [])
+    out = []
+    for net, refined in nets:
+        if net.tree is None:
+            continue
+        for s in net.candidates(centre, radius):
+            p, _ = net.segment_points(s, centre, radius, res_m / 1000.0)
+            if not len(p):
+                continue
+            lat = np.degrees(np.arcsin(np.clip(p[:, 2], -1, 1)))
+            lon = np.degrees(np.arctan2(p[:, 1], p[:, 0]))
+            x, y = to_local(lat0, lon0, lat, lon, R)
+            keep = (np.abs(x) <= half_m) & (np.abs(y) <= half_m)
+            if not rs.empty:
+                w = rs.weight(p)
+                keep &= (w >= 0.5) if refined else (w < 0.5)
+            runs = np.split(np.arange(len(p)), np.where(~keep)[0])
+            q = (net.half_w[s] / 0.004) ** 2 / 31.7                # back from the hydraulic width (rivers.half_width_km)
+            for run in runs:
+                run = run[keep[run]]
+                if len(run) >= 2:
+                    out.append({"points": [[round(float(a), 1), round(float(b), 1)] for a, b in zip(x[run], y[run])],
+                                "width_m": round(float(2000.0 * net.half_w[s]), 1), "discharge_km3_yr": round(float(q), 3),
+                                "refined": refined})
+    out.sort(key=lambda r: (r["points"][0], r["width_m"]))
+    return out
+
+
+def _hillshade(h, res_m, exag):
+    gy, gx = np.gradient(np.asarray(h, dtype=np.float32), np.float32(res_m))
+    dzdx, dzdn = gx * exag, -gy * exag                          # rows run north → south
+    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))
+    return np.clip((-dzdx * L[0] - dzdn * L[1] + L[2]) / np.sqrt(dzdx ** 2 + dzdn ** 2 + 1.0), 0.0, 1.0)
+
+
+def _hillshade_exaggeration(h, res_m):
+    """Vertical exaggeration for the preview: gentle land steepened until its steeper slopes read (≤ 50×)."""
+    step = max(1, h.shape[0] // 1024)                           # a sample of the slopes is enough
+    hh = np.asarray(h[::step, ::step], dtype=np.float64)
+    gy, gx = np.gradient(hh, res_m * step)
+    p95 = float(np.percentile(np.hypot(gx, gy), 95))
+    return float(np.clip(0.35 / max(p95, 1e-6), 1.0, 50.0))
+
+
+MAX_TILES = 1600        # ≈ an 80 km square at 10 m near the equator (≈ 1 GB of tile data): more is refused
+POLE_LIMIT_DEG = 85.0   # squares reaching closer to a pole are refused (tiles narrow toward the poles)
+
+
+def plan(lat, lon, size_km, res_m, radius_km, pixels=None) -> dict:
+    """Samples, tile zoom and the tile rectangle an export needs (ValueError near a pole)."""
+    n = int(pixels) if pixels else samples_for(size_km, res_m)
+    half = (n - 1) / 2 * res_m
+    reach = math.degrees(half * math.sqrt(2) / (radius_km * 1000.0))
+    if abs(lat) + reach > POLE_LIMIT_DEG:
+        raise ValueError(f"the square reaches within {90 - POLE_LIMIT_DEG:g}° of a pole: move it or make it smaller")
+    # the tile zoom whose pixels are at least as fine as res_m (north-south; east-west pixels are finer still)
+    z = int(min(T.MAX_Z, max(5, math.ceil(math.log2(math.pi * radius_km * 1000.0 / (T.TILE * res_m))))))
+    span = 180.0 / 2 ** z
+    edge = np.r_[0, n - 1]                                      # the square's extremes lie on its edges
+    la1, lo1 = grid_latlon(lat, lon, n, res_m, radius_km, rows=(0, 1))
+    la2, lo2 = grid_latlon(lat, lon, n, res_m, radius_km, rows=(n - 1, n))
+    la3, lo3 = (np.concatenate(v) for v in zip(*(grid_latlon(lat, lon, n, res_m, radius_km, rows=(k, k + 1))
+                                                  for k in range(0, n, max(1, (n - 1) // 64)))))
+    la = np.concatenate([la1.ravel(), la2.ravel(), la3[:, edge].ravel()])
+    lo = np.concatenate([lo1.ravel(), lo2.ravel(), lo3[:, edge].ravel()])
+    lonu = lon + ((lo - lon + 180.0) % 360.0 - 180.0)
+    fx, fy = (lonu + 180.0) / span, (90.0 - la) / span
+    x0, x1 = int(math.floor(fx.min())) - 1, int(math.floor(fx.max())) + 1
+    y0, y1 = max(0, int(math.floor(fy.min())) - 1), min(2 ** z - 1, int(math.floor(fy.max())) + 1)
+    return {"n": n, "z": z, "span": span, "x0": x0, "x1": x1, "y0": y0, "y1": y1,
+            "tiles": (x1 - x0 + 1) * (y1 - y0 + 1)}
+
+
+def _parallel(jobs, fn, threads, put, stop):
+    """fn(job) on daemon threads (never delaying the server's exit), put(job, result) in turn; the first error wins."""
+    import threading
+    it, lock, err = iter(jobs), threading.Lock(), []
+
+    def run():
+        while not err:
+            with lock:
+                job = next(it, None)
+            if job is None:
+                return
+            try:
+                r = fn(job)
+                with lock:
+                    put(job, r)
+                stop()
+            except BaseException as e:                          # noqa: BLE001 — re-raised below
+                err.append(e)
+    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()
+    if err:
+        raise err[0]
+
+
+def export(root: Path, res: int, lat: float, lon: float, size_km: float = 40.0, res_m: float = 10.0,
+           out_dir: Path | None = None, name: str | None = None, pixels: int | None = None,
+           pins_path: Path | None = None, regions_dir: Path | None = None, workers: int = 0, progress=None,
+           src=None, max_tiles: int = MAX_TILES, era: str | None = None) -> Path:
+    """Write the export folder and return it. progress(done, total) per tile. src: a running server's tile source
+    (its world, regions and render workers are used; nothing is loaded again). RuntimeError if the refined regions
+    change meanwhile (the export would mix two maps). era: a built era (default: [eras] default, else the base);
+    a server's src already is its chosen era's."""
+    global _SRC
+    from mapgen import config as C
+    root = Path(root)
+    cfg = C.load(root)[0]
+    out_dir = Path(out_dir or (root / "exports"))
+    name = name or default_name(lat, lon, size_km, res_m)
+    dest = out_dir / name
+    if dest.exists():
+        raise FileExistsError(f"export {dest} exists: pick another name or remove it")
+    own = src is None
+    if own:
+        import refine
+        built = [n for n, _ in refine.world_dirs(root, res, log=lambda m: None)] if (
+            root / "out" / f"r{res}" / "cells.npz").exists() else ["base"]
+        default = (C.load(root)[1].get("eras") or {}).get("default")
+        name_era = era or (default if default in built else built[0])
+        if name_era not in built:
+            raise SystemExit(f"no built era {name_era!r} at res {res} (have: {built})")
+        world = dict(serve.load_worlds(root, res, log=lambda m: None, only=name_era))[name_era]
+        src = T.TileSource(world, int(cfg["build"]["seed"]),
+                           cache_dir=Path(tempfile.mkdtemp(prefix="export-tiles-")),
+                           regions_dir=regions_dir or (root / "out" / f"r{res}" / "regions"))
+    try:
+        return _export(src, own, root, res, lat, lon, size_km, res_m, out_dir, name, dest, pixels, pins_path,
+                       workers, progress, max_tiles, cfg)
+    finally:
+        if own:
+            shutil.rmtree(src.cache_root, ignore_errors=True)
+
+
+def _export(src, own, root, res, lat, lon, size_km, res_m, out_dir, name, dest, pixels, pins_path, workers,
+            progress, max_tiles, cfg):
+    global _SRC
+    world, R = src.w, src.R
+    P = plan(lat, lon, size_km, res_m, R, pixels)
+    if P["tiles"] > max_tiles:
+        raise ValueError(f"{P['tiles']} tiles is too many (≤ {max_tiles}): a coarser resolution or a smaller square")
+    n, z, x0, x1, y0, y1 = P["n"], P["z"], P["x0"], P["x1"], P["y0"], P["y1"]
+    half = (n - 1) / 2 * res_m
+    rs = src.regions                                            # one region set for the whole export
+    changed = lambda: src.regions is not rs
+    W, H = (x1 - x0 + 1), (y1 - y0 + 1)
+    jobs = [(z, xx % 2 ** (z + 1), yy) for yy in range(y0, y1 + 1) for xx in range(x0, x1 + 1)]
+    mosaic = {"z": np.zeros((H * T.TILE, W * T.TILE), np.float32), "rgb": np.zeros((H * T.TILE, W * T.TILE, 3), np.uint8),
+              **{k: np.zeros((H * T.TILE, W * T.TILE), np.uint8) for k in ("water", "biome", "ground", "landform")}}
+    done = [0]
+
+    def put(job, part):
+        zz, xx, yy = job
+        r, c = (yy - y0) * T.TILE, ((xx - x0) % 2 ** (z + 1)) * T.TILE
+        for k, v in part.items():
+            mosaic[k][r:r + T.TILE, c:c + T.TILE] = v
+        done[0] += 1
+        if progress:
+            progress(done[0], len(jobs))
+
+    def check():
+        if changed():
+            raise RuntimeError("the refined regions changed during the export (a region build finished): export again")
+    src.tree()
+    src.river_net
+    src.raster("elevation")
+    src.regions.warm()                                          # region search trees too: forked workers share them
+    if not own:                                                 # a server: several render workers at once
+        pool = src.pool
+        _parallel(jobs, lambda j: src.export_tile(*j), max(1, min(4, pool.alive // 2)) if pool else 1, put, check)
+    elif workers and workers > 0:
+        _SRC = src
+        with multiprocessing.get_context("fork").Pool(workers) as mp:
+            for job, part in zip(jobs, mp.imap(_tile, jobs, chunksize=1)):
+                put(job, part)
+    else:
+        for job in jobs:
+            put(job, tile_data(src, *job))
+    check()
+    # sample the square in row blocks (memory: the mosaic and the outputs, not n² float64 grids)
+    span = P["span"]
+    h = np.empty((n, n), np.float32)
+    cls = {k: np.empty((n, n), np.uint8) for k in ("water", "biome", "ground", "landform")}
+    rgb = np.empty((n, n, 3), np.uint8)
+    for r0 in range(0, n, 256):
+        r1 = min(n, r0 + 256)
+        la, lo = grid_latlon(lat, lon, n, res_m, R, rows=(r0, r1))
+        lonu = lon + ((lo - lon + 180.0) % 360.0 - 180.0)
+        pr = ((90.0 - la) / span - y0) * T.TILE - 0.5            # mosaic pixel positions (pixel centres)
+        pc = ((lonu + 180.0) / span - x0) * T.TILE - 0.5
+        h[r0:r1] = map_coordinates(mosaic["z"], [pr, pc], order=1, mode="nearest")
+        ri = np.clip(np.rint(pr).astype(np.int32), 0, H * T.TILE - 1)
+        ci = np.clip(np.rint(pc).astype(np.int32), 0, W * T.TILE - 1)
+        for k in cls:
+            cls[k][r0:r1] = mosaic[k][ri, ci]
+        rgb[r0:r1] = mosaic["rgb"][ri, ci]
+    del mosaic
+    out_dir.mkdir(parents=True, exist_ok=True)
+    tmp = Path(tempfile.mkdtemp(prefix=f".{name}-", dir=out_dir))   # same file system: the final rename is atomic
+    try:
+        h.astype("<f4").tofile(tmp / "height.f32")
+        lo_, hi_ = float(h.min()), float(h.max())
+        v = np.rint((h - lo_) / max(hi_ - lo_, 1e-9) * 65535).astype(np.uint16)
+        Image.fromarray(v).save(tmp / "height.png")
+        v.astype("<u2").tofile(tmp / "height.r16")
+        del v
+        for k, f in (("water", "water.png"), ("biome", "biome.png"), ("ground", "ground.png"), ("landform", "landform.png")):
+            Image.fromarray(cls[k]).save(tmp / f)
+        exag = _hillshade_exaggeration(h, res_m)
+        hs = _hillshade(h, res_m, exag)
+        land = cls["water"] != 1
+        for r0 in range(0, n, 512):                             # shade the preview in blocks (float32)
+            r1 = min(n, r0 + 512)
+            sh = np.where(land[r0:r1], 0.55 + 0.45 * hs[r0:r1], 0.85 + 0.15 * hs[r0:r1]).astype(np.float32)
+            rgb[r0:r1] = np.clip(rgb[r0:r1] * sh[..., None], 0, 255).astype(np.uint8)
+        del hs
+        Image.fromarray(rgb).save(tmp / "preview.png")
+        (tmp / "rivers.json").write_text(json.dumps(_rivers(src, rs, lat, lon, half, res_m)))
+        pins = []
+        pp = Path(pins_path) if pins_path else root / "places" / "pins.json"
+        for p in (json.loads(pp.read_text()).get("pins", []) if pp.exists() else []):
+            px, py = to_local(lat, lon, float(p["lat"]), float(p["lon"]), R)
+            if abs(px) <= half and abs(py) <= half:
+                pins.append({"name": p.get("name"), "lore": p.get("lore"), "epoch": p.get("epoch"),
+                             "x": round(float(px), 1), "y": round(float(py), 1), "lat": p["lat"], "lon": p["lon"]})
+        (tmp / "pins.json").write_text(json.dumps(pins, ensure_ascii=False, indent=1))
+        leg = world.legends
+        sea = cls["water"] == 1
+        i0 = world.index_of(lat, lon)
+        zones = {k: round(float(world.arrays[k][i0]), 4) for k in
+                 ("gravity_g", "o2_fraction", "po2_bar", "pressure_bar", "fire_reactivity") if k in world.arrays}
+        vents = []
+        for v in serve.vents_of(src)["vents"]:
+            vx, vy = to_local(lat, lon, v["lat"], v["lon"], R)
+            if abs(vx) <= half and abs(vy) <= half:
+                vents.append({"x": round(float(vx), 1), "y": round(float(vy), 1),
+                              **{k: v[k] for k in ("type", "temp_c", "flow", "mineral")}})
+        (tmp / "legend.json").write_text(json.dumps({"water": WATER, "biome": leg["holdridge"], "ground": leg["ground"],
+                                                     "landform": leg["landform"]}, ensure_ascii=False, indent=1))
+        meta = {"format": "worldmap-terrain-export", "version": FORMAT_VERSION, "name": name,
+                "center": {"lat": lat, "lon": lon}, "samples": n, "res_m": res_m, "size_m": (n - 1) * res_m,
+                "projection": {"kind": "azimuthal equidistant", "sphere_radius_km": R, "x": "east", "y": "north",
+                               "origin": "the centre sample", "row_0": "north edge", "grid": "samples on the edges"},
+                "height": {"file": "height.f32 (float32 LE, m)", "min": lo_, "max": hi_,
+                           "u16": "h = min + v / 65535 * (max - min)  (height.png, height.r16 LE)"},
+                "sea_level_m": 0, "center_height_m": round(float(h[n // 2, n // 2]), 3), "zoom": z,
+                "era": world.era,
+                "water": {"sea_fraction": round(float(sea.mean()), 4),
+                          "max_depth_m": round(float(-h[sea].min()), 1) if sea.any() else 0.0,
+                          "center_depth_m": round(float(max(0.0, -h[n // 2, n // 2])), 1)
+                          if sea[n // 2, n // 2] else 0.0},
+                "vents": vents, "zones": zones,
+                "preview": {"file": "preview.png", "exaggeration": round(exag, 2),
+                            "note": "relief colours, hillshade from the north-west with this vertical exaggeration"},
+                "tile_px_m": round(math.pi * R * 1000.0 / (2 ** z * T.TILE), 3),
+                "sources": {"world_build": src.fingerprint, "regions": "" if rs.empty else rs.fingerprint,
+                            "tile_version": T.VERSION, "region_version": T.REGION_VERSION, "refine_model": refine.MODEL,
+                            "seed": int(cfg["build"]["seed"]), "world_res": res},
+                "created": date.today().isoformat(),
+                "notes": ["Generated terrain: build data ≥ ≈ 34 km, refined regions ≥ ≈ 5 km, procedural detail "
+                          "below — plausible, not surveyed.",
+                          "Heights: ground, lake and river surfaces; the sea floor below 0 m."]}
+        (tmp / "meta.json").write_text(json.dumps(meta, ensure_ascii=False, indent=1))
+        (tmp / "README.txt").write_text(README.format(**{**meta, "size_km": meta["size_m"] / 1000, "lo": lo_, "hi": hi_}))
+        check()
+        os.replace(tmp, dest)                                   # complete or not there at all
+    except BaseException:
+        shutil.rmtree(tmp, ignore_errors=True)
+        raise
+    return dest
+
+
+def default_name(lat, lon, size_km, res_m) -> str:
+    """e.g. s13.400-w30.400-40km-10m"""
+    return f"{'n' if lat >= 0 else 's'}{abs(lat):.3f}-{'e' if lon >= 0 else 'w'}{abs(lon):.3f}-{size_km:g}km-{res_m:g}m"
+
+
+README = """Terrain export "{name}" (format {version})
+
+Square of {size_km:g} km around {center[lat]}, {center[lon]}: {samples} x {samples} samples every {res_m:g} m.
+Row 0 is the north edge, samples lie on the edges (vertex grid); x east, y north from the centre sample
+(azimuthal equidistant on a sphere of {projection[sphere_radius_km]:g} km).
+
+height.f32   float32 little-endian metres (ground, lake and river surfaces; sea floor below 0 m; sea level 0 m)
+height.png   16-bit grayscale, height.r16 16-bit little-endian raw: h = {lo:.2f} + v / 65535 * ({hi:.2f} - {lo:.2f}) m
+water.png    0 land, 1 sea, 2 lake, 3 river channel
+biome.png, ground.png, landform.png   class indices, names in legend.json
+preview.png  the map's relief colours with hillshade
+rivers.json  river centrelines [{{points: [[x, y], ...] m, width_m, discharge_km3_yr, refined}}]
+pins.json    map pins in the square (x, y in m)
+meta.json    all of the above, sources and versions; the era, water depth, vents in the square (x, y m; types in the
+             viewer legend) and the centre's gravity / O₂ / air pressure / fire reactivity
+
+Import hints: Unity: Terrain > Import Raw, height.r16, {samples} x {samples}, 16 bit, byte order Windows (little),
+terrain size {size_m:g} x {size_m:g} m, height = max - min ({hi:.2f} - {lo:.2f}), and the terrain's Y position = min
+({lo:.2f}) so heights come out in metres (sea level at Y = 0). Godot (Terrain3D / HTerrain): height.png or height.f32
+(heights in metres). Unreal: height.png (resample to a landscape size such as 4033 or 8129 if asked).
+
+Own engines: height.f32 is a plain row-major grid (x east, row 0 north); subtract center_height_m (meta.json) to put
+the centre at 0 m.
+
+Generated terrain.
+"""