diff options
Diffstat (limited to 'tiles.py')
| -rw-r--r-- | tiles.py | 1160 |
1 files changed, 1160 insertions, 0 deletions
diff --git a/tiles.py b/tiles.py new file mode 100644 index 0000000..2b54a26 --- /dev/null +++ b/tiles.py @@ -0,0 +1,1160 @@ +"""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 |