diff options
Diffstat (limited to 'export.py')
| -rw-r--r-- | export.py | 398 |
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. +""" |