From 3443c1c65e9f1753e1e656b35d08416c1fa298f2 Mon Sep 17 00:00:00 2001 From: godosa Date: Wed, 7 Oct 2026 00:14:38 +0200 Subject: worldmap-viewer: initial public history --- export.py | 398 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 398 insertions(+) create mode 100644 export.py (limited to 'export.py') 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(" 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. +""" -- cgit