"""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(" 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