diff options
Diffstat (limited to 'mapgen')
| -rw-r--r-- | mapgen/__init__.py | 7 | ||||
| -rw-r--r-- | mapgen/check.py | 272 | ||||
| -rw-r--r-- | mapgen/climate.py | 249 | ||||
| -rw-r--r-- | mapgen/config.py | 269 | ||||
| -rw-r--r-- | mapgen/crust.py | 97 | ||||
| -rw-r--r-- | mapgen/elevation.py | 197 | ||||
| -rw-r--r-- | mapgen/environment.py | 158 | ||||
| -rw-r--r-- | mapgen/eras.py | 219 | ||||
| -rw-r--r-- | mapgen/erosion.py | 54 | ||||
| -rw-r--r-- | mapgen/events.py | 87 | ||||
| -rw-r--r-- | mapgen/fields.py | 36 | ||||
| -rw-r--r-- | mapgen/geo.py | 200 | ||||
| -rw-r--r-- | mapgen/graph.py | 477 | ||||
| -rw-r--r-- | mapgen/grid.py | 121 | ||||
| -rw-r--r-- | mapgen/hydrology.py | 231 | ||||
| -rw-r--r-- | mapgen/ice.py | 37 | ||||
| -rw-r--r-- | mapgen/minerals.py | 172 | ||||
| -rw-r--r-- | mapgen/newworld.py | 117 | ||||
| -rw-r--r-- | mapgen/noise.py | 111 | ||||
| -rw-r--r-- | mapgen/ocean.py | 166 | ||||
| -rw-r--r-- | mapgen/pipeline.py | 180 | ||||
| -rw-r--r-- | mapgen/plateaus.py | 179 | ||||
| -rw-r--r-- | mapgen/plates.py | 86 | ||||
| -rw-r--r-- | mapgen/projections.py | 150 | ||||
| -rw-r--r-- | mapgen/render.py | 526 | ||||
| -rw-r--r-- | mapgen/seabed.py | 95 | ||||
| -rw-r--r-- | mapgen/sketch.py | 121 | ||||
| -rw-r--r-- | mapgen/sphere.py | 68 | ||||
| -rw-r--r-- | mapgen/testdata/world.toml | 31 | ||||
| -rw-r--r-- | mapgen/testing.py | 59 | ||||
| -rw-r--r-- | mapgen/viewer_export.py | 39 | ||||
| -rw-r--r-- | mapgen/zones.py | 38 |
32 files changed, 4849 insertions, 0 deletions
diff --git a/mapgen/__init__.py b/mapgen/__init__.py new file mode 100644 index 0000000..63fd969 --- /dev/null +++ b/mapgen/__init__.py @@ -0,0 +1,7 @@ +"""World generator.""" +import os + +# OpenBLAS threads spin ~0.1 s after each call before sleeping; with solves on several threads (and other jobs on the +# machine) that spinning steals the cores doing the work. A short timeout changes nothing in how work is split, so +# results stay the same bits (the thread *count* would not: see pipeline/graph notes). Only before numpy loads. +os.environ.setdefault("OPENBLAS_THREAD_TIMEOUT", "4") diff --git a/mapgen/check.py b/mapgen/check.py new file mode 100644 index 0000000..f3f4639 --- /dev/null +++ b/mapgen/check.py @@ -0,0 +1,272 @@ +"""`mapgen.py check`: physical sanity checks (fail) + environment summary (info).""" +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np + +from . import plateaus as PL, zones as ZN +from .climate import DEFAULTS as CLIMATE_DEFAULTS +from .config import params +from .environment import GROUND_NAMES, LANDFORM_NAMES, REGIONS +from .fields import gravity_mod +from .graph import gradient +from .ice import ICE_NAMES +from .sphere import gc_dist_km, latlon_to_xyz + + +def check_land_fraction(z, area, target, tol=0.02): + f = area[z > 0].sum() / area.sum() + return None if abs(f - target) <= tol else f"land fraction {f:.3f} not within {tol} of target {target}" + + +def check_hypsometry(z, area): + edges = np.arange(-11000, 12500, 500) + h, _ = np.histogram(z, bins=edges, weights=area) + c = (edges[:-1] + edges[1:]) / 2 + ocean = (c >= -7000) & (c <= -1000) + landb = (c >= -500) & (c <= 3000) + io, il = np.flatnonzero(ocean)[np.argmax(h[ocean])], np.flatnonzero(landb)[np.argmax(h[landb])] + valley = h[io:il + 1].min() + if h[io] > 1.5 * valley and h[il] > 1.5 * valley: + return None + return "hypsometry not bimodal (expected separate ocean-floor and continent peaks)" + + +def check_max_elevation(z, gmod, cap=12000.0): + lim = cap * np.clip(1.0 / gmod, 1.0, 3.0) + bad = z > lim + 1.0 + return None if not bad.any() else f"{int(bad.sum())} cells above {cap:.0f} m (outside low-g zones)" + + +def check_drainage(recv, z, endorheic, terminal_ok=None): + roots = recv == np.arange(len(recv)) + bad = roots & (z > 0) & ~endorheic + if bad.any(): + return f"{int(bad.sum())} land cells are sinks outside endorheic basins" + if terminal_ok is not None: + dry = roots & endorheic & ~terminal_ok + if dry.any(): + return f"{int(dry.sum())} endorheic terminals are neither lake nor salt flat" + return None + + +def check_no_uphill(zf, recv, z, endorheic, lake=None): + ar = np.arange(len(recv)) + sel = (z > 0) & ~endorheic & (recv != ar) + if lake is not None: + sel &= ~lake # flat lake surfaces: routing direction is a convention + bad = sel & (zf[recv] >= zf) + return None if not bad.any() else f"{int(bad.sum())} cells route uphill on the filled surface" + + +def _band_means(g, f, edges): + a = np.abs(g.lat) + return [np.sum(f[(a >= lo) & (a < hi)] * g.area_km2[(a >= lo) & (a < hi)]) / + max(g.area_km2[(a >= lo) & (a < hi)].sum(), 1e-9) for lo, hi in zip(edges, edges[1:])] + + +def check_temperature(g, t_mean): + m = _band_means(g, t_mean, list(range(0, 91, 10))) + ok = all(b <= a + 0.5 for a, b in zip(m, m[1:])) and m[0] > m[-1] + 10 + return None if ok else f"zonal-mean temperature does not fall poleward: {np.round(m, 1).tolist()}" + + +def check_rain_bands(g, p_ann, hadley_edge): + """Dry belt ≈ Hadley edge −2..+8°, storm track ≈ edge +12..+25° (scaled from Earth's 30° edge).""" + h = hadley_edge + eq = _band_means(g, p_ann, [0, 10])[0] + sub = _band_means(g, p_ann, [h - 2, h + 8])[0] + mid = _band_means(g, p_ann, [h + 12, h + 25])[0] + if eq > sub and mid > sub: + return None + return (f"no subtropical dry belt: P(0-10)={eq:.0f}, P({h - 2:.0f}-{h + 8:.0f})={sub:.0f}, " + f"P({h + 12:.0f}-{h + 25:.0f})={mid:.0f} mm/yr") + + +def check_rain_shadow(g, data): + z = np.maximum(np.asarray(data["elevation_eroded_m"], dtype=np.float64), 0.0) / 1000.0 + grad = gradient(g, z) + land = z > 0 + ww, lw = [], [] + for s in ("jun", "dec"): + up = np.sum(np.asarray(data[f"wind_{s}"], dtype=np.float64) * grad, axis=1) + p = np.asarray(data[f"P_{s}"]) + ww.append(p[land & (up > 0.02)]) + lw.append(p[land & (up < -0.02)]) + w, l = np.concatenate(ww), np.concatenate(lw) + if len(w) == 0 or len(l) == 0: + return None + return None if w.mean() > l.mean() else f"windward slopes ({w.mean():.0f}) not wetter than lee ({l.mean():.0f})" + + +def _shares(values, area, names, mask): + tot = area[mask].sum() + out = [] + for i, nm in enumerate(names): + s = area[mask & (values == i)].sum() / max(tot, 1e-9) + if s >= 0.005: + out.append(f"{nm} {100 * s:.1f}%") + return ", ".join(out) + + + +def check_plateaus(g, data, plateaus, seed) -> list: + """Spec §8: every plateau's top within top_m (± relief), hidden ones ≤ −500 m, island ones with land.""" + if not plateaus: + return [] + z = np.asarray(data["elevation_eroded_m"], dtype=np.float64) + ocean, ids = np.asarray(data["ocean"]), np.asarray(data["plateau_id"]) + out = [] + for k, p in enumerate(plateaus): + idx = np.flatnonzero(ids == k) + if len(idx) == 0: + out.append(f"plateau {p['name']}: no cells at this resolution") + continue + _, b = PL.semi_axes(p) + core = idx[(1.0 - PL.rho(g.xyz[idx], p, seed, g.radius_km)) * b > PL.MARGIN_KM] + lo, hi = p["top_m"] + if len(core): + med = float(np.median(-z[core])) + if not lo - PL.RELIEF_M <= med <= hi + PL.RELIEF_M: + out.append(f"plateau {p['name']}: median top depth {med:.0f} m outside {lo:.0f}–{hi:.0f} m") + if p.get("islands", False): + isl = [g.cell_index(c["lat"], c["lon"]) for c in PL.features(p, seed, g.radius_km)["cones"] if c["island"]] + if ocean[idx].all() and all(ocean[i] for i in isl): + out.append(f"plateau {p['name']}: no island breaks the surface") + elif z[idx].max() > -500.0: + out.append(f"plateau {p['name']}: a hidden plateau reaches {z[idx].max():.0f} m (limit −500 m)") + return out + + +def check_zones(g, data, max_o2_points=0.23, max_g=0.014, per_km=30.0) -> list: + """Spec §5, the gradual rule: steepest change per day's walk of the O₂ fraction (points) and of gravity.""" + out = [] + for key, scale, lim, unit in (("o2_fraction", 100.0, max_o2_points, "O₂ points"), ("gravity_g", 1.0, max_g, "g")): + if key not in data: + continue + f = np.asarray(data[key], dtype=np.float64) * scale + steep = float(np.max(np.abs(f[g.dst] - f[g.src]) / g.edge_km)) * per_km + if steep > lim + 1e-9: + out.append(f"zones: {key} changes {steep:.3f} {unit} per {per_km:.0f} km (limit {lim})") + return out + + +def edited_areas(xyz, tect, radius_km, plateau_km=300.0, patch_km=1500.0): + """Cells the revision changes by design: plateaus + 300 km, land patches + 1,500 km, low-gravity zones + (their relief).""" + xyz = np.asarray(xyz, dtype=np.float64) + near = lambda c, r: gc_dist_km(xyz, latlon_to_xyz(*c), radius_km) < r + m = np.zeros(len(xyz), bool) + for p in tect.get("plateau", []): + m |= near(p["center"], PL.semi_axes(p)[0] * (1 + PL.EDGE_WARP) / (1 - PL.EDGE_WARP) + plateau_km) + for p in tect.get("land_patch", []): + m |= near(p["center"], p["radius_km"] * 1.3 + patch_km) + for z in tect.get("zone", []): + if z["field"] == "gravity": + m |= near(z["center"], z["radius_km"] / (1 - ZN.WARP)) + return m + + +O2_FIELDS = ("m_o2_zones", "o2_fraction", "po2_bar") # what the base's O₂ zones change (no relief) + + +def o2_areas(xyz, tect, radius_km): + """Cells the base's O₂ zones reach: compare skips only the O₂ fields there.""" + xyz = np.asarray(xyz, dtype=np.float64) + m = np.zeros(len(xyz), bool) + for z in tect.get("zone", []): + if z["field"] == "o2": + m |= gc_dist_km(xyz, latlon_to_xyz(*z["center"]), radius_km) < z["radius_km"] / (1 - ZN.WARP) + return m + + +def compare(old: dict, new: dict, exclude, extra=None) -> list: + """Per field, over the cells outside `exclude` (and outside extra[field] for that field): share changed, largest + change, 99th percentile.""" + keep0 = ~np.asarray(exclude, bool) + lines = [f"cells compared: {int(keep0.sum())} of {len(keep0)} (outside the edited areas)"] + for k in sorted(set(old) & set(new)): + a, b = np.asarray(old[k]), np.asarray(new[k]) + if a.shape != b.shape or a.shape[:1] != keep0.shape or a.dtype.kind not in "biuf" or b.dtype.kind not in "biuf": + continue + keep = keep0 & ~np.asarray(extra[k], bool) if extra and k in extra else keep0 + a, b = a.astype(np.float64)[keep], b.astype(np.float64)[keep] + with np.errstate(invalid="ignore"): + diff = np.where((a == b) | (np.isnan(a) & np.isnan(b)), 0.0, np.abs(b - a)) + if diff.ndim > 1: + diff = diff.reshape(len(diff), -1).max(axis=1) + if not len(diff): + lines.append(f"{k}: no cells to compare") + continue + lines.append(f"{k}: {np.mean(diff > 0):.1%} changed, max {diff.max():.4g}, p99 {np.percentile(diff, 99):.4g}") + return lines + + +def compare_report(root: Path, res: int, old_path: Path) -> int: + from . import config as C + _, tect = C.load(root) + out = root / "out" / f"r{res}" + with np.load(old_path) as z: + old = {k: z[k] for k in z.files} + with np.load(out / "cells.npz") as z: + new = {k: z[k] for k in z.files} + for name, arrays in ((old_path, old), (out / "cells.npz", new)): + if "g_ids" not in arrays: + print(f"error: {name} has no g_ids: not the cells.npz of a world build") + return 2 + if not np.array_equal(old["g_ids"], new["g_ids"]): + print("error: the two worlds have different cells (another resolution?)") + return 2 + radius = float(json.loads((out / "cells_meta.json").read_text())["radius_km"]) + o2 = o2_areas(new["g_xyz"], tect, radius) + for line in compare(old, new, edited_areas(new["g_xyz"], tect, radius), {k: o2 for k in O2_FIELDS}): + print(line) + return 0 + +def info_lines(ctx): + g, d = ctx.grid, ctx.data + z = np.asarray(d["elevation_eroded_m"]) + land = ~np.asarray(d["ocean"]) if "ocean" in d else z > 0 + a = g.area_km2 + sk = np.asarray(d["sk_land"]) > 0.5 + iou = (land & sk).sum() / max((land | sk).sum(), 1) + sites = PL.site_report(g, d, ctx.tect["plateau"]) if ctx.tect.get("plateau") and "plateau_id" in d else [] + return [f"cells {g.n} (H3 res {ctx.res}); land {100 * a[land].sum() / a.sum():.1f}% " + f"({a[land].sum() / 1e6:.0f}M km², Earth: 29%, 149M km²)", + f"sketch land overlap (IoU): {iou:.2f}", + "regions: " + _shares(np.asarray(d["hold_region"]), a, REGIONS, land), + "landforms: " + _shares(np.asarray(d["landform"]), a, LANDFORM_NAMES, land), + "ground: " + _shares(np.asarray(d["ground"]), a, GROUND_NAMES, land), + "ice: " + _shares(np.asarray(d["ice"]), a, ICE_NAMES, np.ones(g.n, bool))] + [f"site: {s}" for s in sites] + + +def run_checks(ctx): + g, d, cfg = ctx.grid, ctx.data, ctx.cfg + z = np.asarray(d["elevation_eroded_m"], dtype=np.float64) + if "ocean" in d: # interior lows below sea level are land: sign-correct z for the land/sea checks + z = np.where(np.asarray(d["ocean"]), np.minimum(z, -0.001), np.maximum(z, 0.001)) + endo = np.asarray(d["endorheic"]) + results = [ + check_land_fraction(z, g.area_km2, cfg["build"]["land_fraction"]), + check_hypsometry(z, g.area_km2), + check_max_elevation(z, gravity_mod(cfg, d["m_gravity_zones"])), + check_drainage(np.asarray(d["recv"]), z, endo, np.asarray(d["lake"]) | np.asarray(d["salt_flat"])), + check_no_uphill(np.asarray(d["z_filled_m"]), np.asarray(d["recv"]), z, endo, np.asarray(d["lake"])), + check_temperature(g, np.asarray(d["T_mean"])), + check_rain_bands(g, np.asarray(d["P_ann"]), params(cfg, "climate", CLIMATE_DEFAULTS)["hadley_edge_deg"]), + check_rain_shadow(g, d), + ] + results += check_plateaus(g, d, ctx.tect.get("plateau", []), ctx.seed) + check_zones(g, d) + return [r for r in results if r], info_lines(ctx) + + +def report(ctx) -> int: + fails, infos = run_checks(ctx) + for line in infos: + print("info:", line) + for f in fails: + print("FAIL:", f) + print("check:", "OK" if not fails else f"{len(fails)} failure(s)") + return 1 if fails else 0 diff --git a/mapgen/climate.py b/mapgen/climate.py new file mode 100644 index 0000000..491675c --- /dev/null +++ b/mapgen/climate.py @@ -0,0 +1,249 @@ +"""Stage `climate`: seasonal insolation → temperature, 3-cell winds + monsoons, moisture transport → rain.""" +from __future__ import annotations + +import numpy as np +from scipy import sparse + +from . import ocean as OC +from .config import params +from .graph import bicgstab_jacobi, distance_to, gradient, nearest_source, pmap, smooth_km +from .grid import rowdot_at +from .sphere import east_north, latlon_to_xyz, rotate_about, tangent_dir + +S0 = 1361.0 +SEA_FREEZE_C = -1.8 # the exported SST never goes below sea water's freezing point (ice-covered sea) +SEASONS = ("jun", "dec", "eq") +DEFAULTS = { + "t_a": -55.2, "t_b": 0.3206, "t_c": -2.974e-4, "land_seasonal": 0.45, "land_summer": 0.9, "continentality_summer_max": 1.5, "ocean_seasonal": 0.15, + "continentality_km": 1500.0, "continentality_max": 1.8, "lapse_c_per_km": 6.5, + "current_c": 4.0, "current_reach_km": 600.0, "current_leak_km": 150.0, "heat_transport_km": 300.0, + "hadley_edge_deg": 20.0, "ferrel_edge_deg": 55.0, "itcz_shift_deg": 8.0, + "trade_u": -6.0, "trade_v": 2.0, "westerly_u": 8.0, "westerly_v": 1.0, "polar_u": -4.0, "polar_v": 1.0, + "monsoon_k": 1.0, "monsoon_length_km": 1500.0, "monsoon_speed_scale": 6000.0, + "coriolis_min_deg": 20.0, "coriolis_span_deg": 50.0, + "base_rate": 0.25, "conv_rate": 2.0, "front_rate": 0.8, "front_lat_deg": 40.0, "front_width_deg": 10.0, "itcz_width_deg": 8.0, "oro_rate": 20.0, "subsidence": 0.2, + "min_rate": 0.05, "cc_per_c": 0.07, "recycle": 0.83, "eddy_k_m2s": 2.2e6, + "eddy_wind_ms": 8.0, + "global_mean_mm": 1000.0, + "lock": False, "lock_at": [0.0, 0.0], "lock_day_c": 120.0, "lock_night_c": -200.0, "lock_wind_ms": 10.0, + "lock_melt_c": 0.0, "lock_melt_width_c": 25.0, +} + + +def t_of_q(q, P): + """Radiative-equilibrium-like surface temperature (°C) from insolation (W/m²); quadratic fit to Earth-like + zonal means for a 20° tilt: equator 27, 45° ≈ 15, 60° ≈ 3, pole ≈ −22 (concave: damps polar-day summers).""" + return P["t_a"] + P["t_b"] * q + P["t_c"] * q * q + + +def declinations(tilt): + return {"jun": tilt, "dec": -tilt, "eq": 0.0} + + +def insolation(lat, decl): + """Daily-mean top-of-atmosphere insolation (W/m²).""" + phi = np.radians(np.clip(lat, -89.9999, 89.9999)) + d = np.radians(decl) + h0 = np.arccos(np.clip(-np.tan(phi) * np.tan(d), -1.0, 1.0)) + return S0 / np.pi * (h0 * np.sin(phi) * np.sin(d) + np.cos(phi) * np.cos(d) * np.sin(h0)) + + +def zonal_mean(g, f, bin_deg=2.0): + b = np.floor((g.lat + 90.0) / bin_deg).astype(np.int64) + s = np.bincount(b, weights=f * g.area_km2) + w = np.bincount(b, weights=g.area_km2) + return (s / np.maximum(w, 1e-12))[b] + + +def current_anomaly(g, land, P): + """Subtropical gyres: cold water off west coasts, warm off east coasts; leaks onto coastal land.""" + if not land.any(): + return np.zeros(g.n) + d, src = nearest_source(g, np.flatnonzero(land)) + e, _ = east_north(g.xyz) + east_comp = np.sum(tangent_dir(g.xyz, g.xyz[np.maximum(src, 0)]) * e, axis=1) + band = np.sin(np.radians((np.clip(np.abs(g.lat), 10.0, 50.0) - 10.0) * 4.5)) + a = -P["current_c"] * np.sign(east_comp) * band * np.exp(-d / P["current_reach_km"]) + a[land] = 0.0 + return smooth_km(g, a, P["current_leak_km"]) + + +def temperatures(g, z, land, P, tilt, cur=None): + decl = declinations(tilt) + Q = {s: insolation(g.lat, d) for s, d in decl.items()} + q_ann = (Q["jun"] + Q["dec"] + 2 * Q["eq"]) / 4 + t_ann = t_of_q(q_ann, P) + dist_ocean = distance_to(g, ~land) if (~land).any() else np.full(g.n, 1e4) + cont = np.clip(1.0 + dist_ocean / P["continentality_km"], 1.0, P["continentality_max"]) + cur = current_anomaly(g, land, P) if cur is None else cur + lapse = -P["lapse_c_per_km"] * np.maximum(z, 0.0) / 1000.0 + def season(s): + raw = t_of_q(Q[s], P) + land_resp = np.where(raw > t_ann, P["land_summer"] * np.minimum(cont, P["continentality_summer_max"]), + P["land_seasonal"] * cont) # land heats faster in summer (dry, low heat capacity) + resp = np.where(land, land_resp, P["ocean_seasonal"]) + return smooth_km(g, t_ann + resp * (raw - t_ann) + cur, P["heat_transport_km"]) + lapse + return dict(zip(SEASONS, pmap(season, SEASONS))), dist_ocean + + +def biotemperature(tmean, trange, n=12): + ph = np.linspace(0.0, 2 * np.pi, n, endpoint=False) + t = np.asarray(tmean)[:, None] + (np.asarray(trange)[:, None] / 2) * np.sin(ph)[None, :] + return np.clip(t, 0.0, 30.0).mean(axis=1) + + +def band_winds(g, itcz_lat, P): + """3-cell surface winds (m/s) relative to the thermal equator.""" + phi = g.lat - itcz_lat + a = np.abs(phi) + sgn = np.where(phi >= 0, 1.0, -1.0) + s1 = 0.5 * (1 + np.tanh((a - P["hadley_edge_deg"]) / 3.0)) + s2 = 0.5 * (1 + np.tanh((a - P["ferrel_edge_deg"]) / 4.0)) + u = (1 - s1) * P["trade_u"] + (s1 - s2) * P["westerly_u"] + s2 * P["polar_u"] + v = sgn * (-(1 - s1) * P["trade_v"] + (s1 - s2) * P["westerly_v"] - s2 * P["polar_v"]) + e, n = east_north(g.xyz) + return u[:, None] * e + v[:, None] * n + + +def monsoon_winds(g, T_s, land, P): + """Thermal lows over hot land / highs over cold land, flow deflected by Coriolis.""" + anom = np.where(land, T_s - zonal_mean(g, T_s), 0.0) + press = -P["monsoon_k"] * smooth_km(g, anom, P["monsoon_length_km"]) + flow = -gradient(g, press) + theta = np.radians(P["coriolis_min_deg"] + P["coriolis_span_deg"] * np.abs(np.sin(np.radians(g.lat)))) + return rotate_about(g.xyz, flow, -np.sign(g.lat) * theta) * P["monsoon_speed_scale"] + + +def _rate(per_1000km, spacing_km): + return 1.0 - np.exp(-np.maximum(per_1000km, 0.0) * spacing_km / 1000.0) + + +def eddy_mixing(g, P): + """kappa · (nbr − I): the per-step eddy mixing with the neighbours — the same for every season (build once).""" + kappa = P["eddy_k_m2s"] / (P["eddy_wind_ms"] * g.spacing_km * 1000.0) + nbr = sparse.csr_matrix((1.0 / g.counts[g.src], (g.src, g.dst)), shape=(g.n, g.n)) + return kappa * (nbr - sparse.identity(g.n, format="csr")) + + +def precipitation(g, wind, T, z, land, itcz_lat, P, mixing=None): + """Steady-state moisture transport along the wind on the cell graph; returns rain (relative units). + mixing: eddy_mixing(g, P), when several seasons share it.""" + t = g.edge_tangents + out = np.maximum(rowdot_at(wind, g.src, t), 0.0) + tot = np.bincount(g.src, weights=out, minlength=g.n) + frac = np.where(tot[g.src] > 0, out / np.maximum(tot[g.src], 1e-12), 0.0) + Tm = sparse.csr_matrix((frac, (g.dst, g.src)), shape=(g.n, g.n)) + stay = (tot <= 0).astype(np.float64) + upslope = np.maximum(np.sum(wind * gradient(g, np.maximum(z, 0.0) / 1000.0), axis=1), 0.0) + conv = P["conv_rate"] * np.exp(-((g.lat - itcz_lat) / P["itcz_width_deg"]) ** 2) * np.clip((T - 10.0) / 20.0, 0, 1) + subs = P["subsidence"] * np.exp(-((np.abs(g.lat - itcz_lat) - P["hadley_edge_deg"]) / 6.0) ** 2) + front = P["front_rate"] * np.exp(-((np.abs(g.lat - itcz_lat) - P["front_lat_deg"]) / P["front_width_deg"]) ** 2) + per = np.maximum(P["base_rate"] + conv + front + P["oro_rate"] * upslope - subs, P["min_rate"]) + r = _rate(per, g.spacing_km) + evap = np.where(land, 0.0, np.exp(P["cc_per_c"] * (np.clip(T, -2.0, 35.0) - 25.0))) + # per advection step (one cell, time h/U): rain out, move downwind, eddy-mix with neighbours + mix = eddy_mixing(g, P) if mixing is None else mixing + eye = sparse.identity(g.n, format="csr") + keep = sparse.diags(1.0 - r) + recyc = sparse.diags(np.where(land, P["recycle"] * r, 0.0)) # land evapotranspiration returns rain + system = (eye - (Tm @ keep + sparse.diags(stay) @ keep) - mix - recyc).tocsr() + W, info = bicgstab_jacobi(system, evap, evap / np.maximum(r, 1e-6), 1.0 / system.diagonal(), 1e-7, 5000) + if info != 0: + raise ValueError(f"precipitation: moisture solve did not converge (info={info})") + return r * np.maximum(W, 0.0) + + +def run_locked(ctx, g, z, land, P) -> dict: + """One face always to the sun: temperature by sun angle, no seasons; surface wind from night to the sun point.""" + sub = latlon_to_xyz(*P["lock_at"]) + mu = np.maximum(g.xyz @ sub, 0.0) + t = P["lock_night_c"] + (P["lock_day_c"] - P["lock_night_c"]) * mu ** 0.25 + lapse = -P["lapse_c_per_km"] * np.maximum(z, 0.0) / 1000.0 + T = smooth_km(g, t, P["heat_transport_km"]) + lapse + dist_ocean = distance_to(g, ~land) if (~land).any() else np.full(g.n, 1e4) + wind = P["lock_wind_ms"] * tangent_dir(g.xyz, np.broadcast_to(sub, g.xyz.shape)) + rain = np.exp(-((T - P["lock_melt_c"]) / P["lock_melt_width_c"]) ** 2) # meltwater and frost in the twilight ring + k = P["global_mean_mm"] / max(np.sum(rain * g.area_km2) / g.area_km2.sum(), 1e-12) + out = {} + for s in SEASONS: + out[f"wind_{s}"] = wind.astype(np.float32) + out[f"P_{s}"] = rain * k + out[f"T_{s}"] = T + out["P_ann"] = rain * k + out["T_mean"] = T + out["T_range"] = np.zeros(g.n) + out["T_min"] = T + out["biotemp"] = biotemperature(T, out["T_range"]) + out["PET"] = 58.93 * out["biotemp"] + out["dist_ocean_km"] = dist_ocean + return out + + +def _still_ocean(g, t_mean): + """No circulation (ocean disabled, locked world): zero currents/upwelling/productivity, SST = T_mean.""" + z = np.zeros(g.n, np.float32) + return {"current": np.zeros((g.n, 3), np.float32), "current_speed": z, "sst": np.asarray(t_mean, np.float32), + "upwelling": z.copy(), "productivity": z.copy()} + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "climate", DEFAULTS) + tilt = float(ctx.cfg["planet"]["tilt_deg"]) + (z,) = ctx.need("elevation_eroded_m") + z = z.astype(np.float64) + water = ctx.data.get("open_water", ctx.data.get("ocean")) # big inland basins are water to the air + land = ~np.asarray(water) if water is not None else z > 0 + O = params(ctx.cfg, "ocean", OC.DEFAULTS) + if P["lock"]: + out = run_locked(ctx, g, z, land, P) + out.update(_still_ocean(g, out["T_mean"])) + return out + sea = ~land + coupled = bool(O["enabled"]) and bool(sea.any()) + T, dist_ocean = temperatures(g, z, land, P, tilt, cur=np.zeros(g.n) if coupled else None) + decl = declinations(tilt) + def season_wind(s): + itcz = P["itcz_shift_deg"] * decl[s] / max(tilt, 1e-9) + wind = band_winds(g, itcz, P) + return wind + monsoon_winds(g, T[s], land, P) if s != "eq" else wind + winds = dict(zip(SEASONS, pmap(season_wind, SEASONS))) + if coupled: # one pass: winds from current-free temperatures, then currents carry heat + day = float(ctx.cfg["planet"]["day_hours"]) + w_ann = (winds["jun"] + winds["dec"] + 2 * winds["eq"]) / 4 + u = OC.currents(g, sea, w_ann, O, day) + T_eq = (T["jun"] + T["dec"] + 2 * T["eq"]) / 4 + T_s = OC.sst(g, sea, u, T_eq, O) + cur = smooth_km(g, np.where(sea, T_s - T_eq, 0.0), P["current_leak_km"]) + T, _ = temperatures(g, z, land, P, tilt, cur=cur) + out = {} + mixing = eddy_mixing(g, P) + rain = dict(zip(SEASONS, pmap(lambda s: precipitation(g, winds[s], T[s], z, land, + P["itcz_shift_deg"] * decl[s] / max(tilt, 1e-9), P, mixing), + SEASONS))) + del mixing + for s in SEASONS: + out[f"wind_{s}"] = winds[s].astype(np.float32) + ann = (rain["jun"] + rain["dec"] + 2 * rain["eq"]) / 4 + k = P["global_mean_mm"] / max(np.sum(ann * g.area_km2) / g.area_km2.sum(), 1e-12) + for s in SEASONS: + out[f"P_{s}"] = rain[s] * k + out[f"T_{s}"] = T[s] + out["P_ann"] = ann * k + out["T_mean"] = (T["jun"] + T["dec"] + 2 * T["eq"]) / 4 + out["T_range"] = np.abs(T["jun"] - T["dec"]) + out["T_min"] = np.minimum(T["jun"], T["dec"]) + out["biotemp"] = biotemperature(out["T_mean"], out["T_range"]) + out["PET"] = 58.93 * out["biotemp"] + out["dist_ocean_km"] = dist_ocean + if coupled: + sst_c = np.where(sea, np.maximum(T_s, SEA_FREEZE_C), out["T_mean"]) + w_up = OC.upwelling(g, sea, w_ann, O, day) + out.update({"current": u.astype(np.float32), + "current_speed": np.linalg.norm(u, axis=1).astype(np.float32), + "sst": sst_c.astype(np.float32), + "upwelling": w_up.astype(np.float32), + "productivity": OC.productivity(g, sea, w_up, z, out["T_range"], sst_c, O).astype(np.float32)}) + else: + out.update(_still_ocean(g, out["T_mean"])) + return out diff --git a/mapgen/config.py b/mapgen/config.py new file mode 100644 index 0000000..a52986a --- /dev/null +++ b/mapgen/config.py @@ -0,0 +1,269 @@ +"""Load and validate map/config/*.toml.""" +from __future__ import annotations + +import tomllib +from pathlib import Path + + +class ConfigError(ValueError): + pass + + +MASK_DEFAULTS = {"o2_range": 0.5, "gravity_range": 0.7} + +REQUIRED = { + "planet": { + "radius_km": (1000.0, 100000.0), + "gravity_g": (0.1, 5.0), + "day_hours": (1.0, 10000.0), + "year_days": (1.0, 100000.0), + "tilt_deg": (0.0, 90.0), + "sea_level_pressure_bar": (0.01, 100.0), + "scale_height_km": (1.0, 100.0), + "o2_fraction": (0.0, 1.0), + }, + "build": { + "seed": (0, 2**31 - 1), + "res_dev": (0, 8), + "res_final": (0, 8), + "raster_width": (64, 32768), + "preview_width": (64, 8192), + "land_fraction": (0.01, 0.99), + }, +} + + +def _check_ranges(cfg: dict, schema: dict, src: str) -> None: + for sec, keys in schema.items(): + if sec not in cfg: + raise ConfigError(f"{src}: missing section [{sec}]") + for k, (lo, hi) in keys.items(): + if k not in cfg[sec]: + raise ConfigError(f"{src}: missing [{sec}].{k}") + v = cfg[sec][k] + if isinstance(v, bool) or not isinstance(v, (int, float)): + raise ConfigError(f"{src}: [{sec}].{k} must be a number, got {v!r}") + if not lo <= v <= hi: + raise ConfigError(f"{src}: [{sec}].{k}={v} outside [{lo}, {hi}]") + + +def _latlon(v, what: str, src: str = "tectonics.toml") -> None: + ok = isinstance(v, list) and len(v) == 2 and all(isinstance(x, (int, float)) for x in v) + if not ok or not (-90 <= v[0] <= 90 and -180 <= v[1] <= 360): + raise ConfigError(f"{src}: {what} must be [lat, lon], got {v!r}") + + +def validate_tectonics(t: dict) -> None: + plates = t.get("plate", []) + if len(plates) < 2: + raise ConfigError("tectonics.toml: need at least 2 [[plate]] tables") + seen = set() + for p in plates: + pid = p.get("id", "?") + for k in ("id", "seed", "kind", "motion"): + if k not in p: + raise ConfigError(f"tectonics.toml: plate {pid} missing {k}") + if pid in seen: + raise ConfigError(f"tectonics.toml: duplicate plate id {pid}") + seen.add(pid) + if p["kind"] not in ("continental", "oceanic"): + raise ConfigError(f"tectonics.toml: plate {pid} kind must be continental|oceanic") + _latlon(p["seed"], f"plate {pid}.seed") + m = p["motion"] + if not (isinstance(m, list) and len(m) == 2 and 0 <= m[1] <= 20): + raise ConfigError(f"tectonics.toml: plate {pid}.motion must be [azimuth_deg, speed 0..20 cm/yr]") + for kind in ("lip", "volcano", "hotspot", "microcontinent"): + for x in t.get(kind, []): + for k in ("name", "center"): + if k not in x: + raise ConfigError(f"tectonics.toml: [[{kind}]] missing {k}") + _latlon(x["center"], f"{kind} {x['name']}.center") + if kind in ("lip", "volcano", "microcontinent") and not x.get("radius_km", 0) > 0: + raise ConfigError(f"tectonics.toml: {kind} {x['name']} needs radius_km > 0") + if kind == "hotspot" and not x.get("length_km", 0) > 0: + raise ConfigError(f"tectonics.toml: hotspot {x['name']} needs length_km > 0") + + +ERAS_FILE = "eras.toml" # events and eras: not part of the base world's inputs key (pipeline.inputs_key) +ZONE_FIELDS = {"gravity_g": (0.05, 5.0), "o2_fraction": (0.0, 1.0), "pressure_bar": (0.01, 10.0), + "fire_reactivity": (0.0, 5.0)} # render.CONTINUOUS holds every value in range +PLATEAU_SPACING_KM = 2500.0 +RESERVED_ERA_KEYS = {"order", "default"} + + +def _num(v, what: str, lo: float, hi: float, src: str = "tectonics.toml") -> None: + if isinstance(v, bool) or not isinstance(v, (int, float)) or not lo <= v <= hi: + raise ConfigError(f"{src}: {what} must be a number in [{lo}, {hi}], got {v!r}") + + +def _pair(v, what: str, lo: float, hi: float, src: str = "tectonics.toml") -> None: + ok = isinstance(v, list) and len(v) == 2 and all(isinstance(x, (int, float)) and not isinstance(x, bool) for x in v) + if not ok or not lo <= v[0] <= v[1] <= hi: + raise ConfigError(f"{src}: {what} must be [low, high] within [{lo}, {hi}], got {v!r}") + + +def _named(items: list, kind: str, src: str = "tectonics.toml") -> set: + seen = set() + for x in items: + if not isinstance(x.get("name"), str): + raise ConfigError(f"{src}: [[{kind}]] needs a name") + if x["name"] in seen: + raise ConfigError(f"{src}: duplicate {kind} name {x['name']!r}") + seen.add(x["name"]) + return seen + + +def _event(e: dict) -> None: + w, src, kind = f"event {e['name']}", ERAS_FILE, e.get("kind") + if kind == "disintegrate": + _latlon(e.get("center"), f"{w}.center", src) + _num(e.get("radius_km"), f"{w}.radius_km", 1.0, 20000.0, src) + _num(e.get("depth_m"), f"{w}.depth_m", 1.0, 20000.0, src) + elif kind == "volcano": + _latlon(e.get("center"), f"{w}.center", src) + _num(e.get("radius_km"), f"{w}.radius_km", 1.0, 5000.0, src) + _num(e.get("peak_m"), f"{w}.peak_m", -11000.0, 12000.0, src) + if e.get("shape", "cone") not in ("cone", "shield", "caldera"): + raise ConfigError(f"{src}: {w}.shape must be cone|shield|caldera, got {e.get('shape')!r}") + if e.get("peak_mode", "above") not in ("above", "absolute"): + raise ConfigError(f"{src}: {w}.peak_mode must be above|absolute, got {e.get('peak_mode')!r}") + elif kind == "zone": + f = e.get("fields") + if not isinstance(f, dict) or not f: + raise ConfigError(f"{src}: {w}.fields must set at least one of {sorted(ZONE_FIELDS)}") + for k, v in f.items(): + if k not in ZONE_FIELDS: + raise ConfigError(f"{src}: {w}.fields: unknown field {k!r} (known: {sorted(ZONE_FIELDS)})") + _num(v, f"{w}.fields.{k}", *ZONE_FIELDS[k], src) + shape = e.get("shape") + if shape == "circle": + _latlon(e.get("center"), f"{w}.center", src) + _num(e.get("radius_km"), f"{w}.radius_km", 1.0, 20000.0, src) + elif shape == "landmass": + _latlon(e.get("seed"), f"{w}.seed", src) + _pair(e.get("reach_km"), f"{w}.reach_km", 0.0, 20000.0, src) + else: + raise ConfigError(f"{src}: {w}.shape must be circle|landmass, got {shape!r}") + if "edge_km" in e: + _num(e["edge_km"], f"{w}.edge_km", 0.001, 5000.0, src) + elif e.get("profile") != "smooth": + raise ConfigError(f'{src}: {w} needs edge_km (a sharp edge) or profile = "smooth"') + else: + raise ConfigError(f"{src}: {w}.kind must be disintegrate|zone|volcano, got {kind!r}") + + +def validate_revision(t: dict, radius_km: float) -> None: + """Plateaus, land patches, zones, events and eras.""" + from .sphere import gc_dist_km, latlon_to_xyz + plateaus = t.get("plateau", []) + _named(plateaus, "plateau") + for p in plateaus: + w = f"plateau {p['name']}" + _latlon(p.get("center"), f"{w}.center") + _num(p.get("area_km2"), f"{w}.area_km2", 1.0e4, 5.0e7) + _num(p.get("elongation", 1.0), f"{w}.elongation", 1.0, 10.0) + _num(p.get("azimuth_deg", 0.0), f"{w}.azimuth_deg", -360.0, 360.0) + _pair(p.get("top_m"), f"{w}.top_m (metres below sea level, shallowest first)", 1.0, 11000.0) + if not isinstance(p.get("islands", False), bool): + raise ConfigError(f"tectonics.toml: {w}.islands must be true or false") + _num(p.get("vent", 1.0), f"{w}.vent", 0.0, 1.0) + for i, a in enumerate(plateaus): + for b in plateaus[i + 1:]: + d = float(gc_dist_km(latlon_to_xyz(*a["center"]), latlon_to_xyz(*b["center"]), radius_km)) + if d < PLATEAU_SPACING_KM: + raise ConfigError(f"tectonics.toml: plateaus {a['name']} and {b['name']} are {d:.0f} km apart " + f"(need ≥ {PLATEAU_SPACING_KM:.0f} km)") + _named(t.get("land_patch", []), "land_patch") + for p in t.get("land_patch", []): + w = f"land_patch {p['name']}" + _latlon(p.get("center"), f"{w}.center") + _num(p.get("radius_km"), f"{w}.radius_km", 1.0, 20000.0) + _num(p.get("strength", 1.0), f"{w}.strength", 0.0, 2.0) + _num(p.get("edge_noise", 0.0), f"{w}.edge_noise", 0.0, 1.0) + _named(t.get("zone", []), "zone") + for z in t.get("zone", []): + w = f"zone {z['name']}" + if z.get("field") not in ("o2", "gravity"): + raise ConfigError(f"tectonics.toml: {w}.field must be o2|gravity, got {z.get('field')!r}") + _latlon(z.get("center"), f"{w}.center") + _num(z.get("radius_km"), f"{w}.radius_km", 1.0, 20000.0) + _num(z.get("v"), f"{w}.v", -1.0, 1.0) + if z.get("profile", "smooth") != "smooth": + raise ConfigError(f"tectonics.toml: {w}.profile must be smooth") + names = _named(t.get("event", []), "event", ERAS_FILE) + for e in t.get("event", []): + _event(e) + eras = t.get("eras") + if eras is None: + return + order = eras.get("order") + if not (isinstance(order, list) and order and all(isinstance(n, str) for n in order) + and len(set(order)) == len(order)) or RESERVED_ERA_KEYS & set(order): + raise ConfigError(f"{ERAS_FILE}: [eras].order must list distinct era names (not order/default), got {order!r}") + if eras.get("default") not in order: + raise ConfigError(f"{ERAS_FILE}: [eras].default {eras.get('default')!r} is not in order {order}") + extra = set(eras) - set(order) - RESERVED_ERA_KEYS + if extra: + raise ConfigError(f"{ERAS_FILE}: [eras] tables {sorted(extra)} are not in order {order}") + for n in order: + e = eras.get(n) + if not isinstance(e, dict) or not isinstance(e.get("label", n), str): + raise ConfigError(f"{ERAS_FILE}: missing [eras.{n}] (label, events)") + evs = e.get("events", []) + if not isinstance(evs, list) or not all(isinstance(x, str) for x in evs): + raise ConfigError(f"{ERAS_FILE}: [eras.{n}].events must be a list of event names") + if "years" in e: + _num(e["years"], f"[eras.{n}].years (since the era before it)", 0.0, 1.0e9, ERAS_FILE) + unknown = [x for x in evs if x not in names] + if unknown: + raise ConfigError(f"{ERAS_FILE}: era {n} names unknown events {unknown}") + + +def era_events(tect: dict, name: str) -> list: + """The events of era `name`: those of every era before it in [eras].order, then its own (cumulative).""" + eras = tect.get("eras") or {} + order = eras.get("order", []) + if name not in order: + raise ConfigError(f"{ERAS_FILE}: unknown era {name!r} (eras: {order})") + by_name = {e["name"]: e for e in tect.get("event", [])} + out = [] + for n in order[: order.index(name) + 1]: + out += [by_name[e] for e in eras.get(n, {}).get("events", [])] + return out + + +def load(root: Path) -> tuple[dict, dict]: + wpath = root / "config" / "world.toml" + tpath = root / "config" / "tectonics.toml" + for p in (wpath, tpath): + if not p.exists(): + raise ConfigError(f"missing {p}") + with open(wpath, "rb") as f: + cfg = tomllib.load(f) + with open(tpath, "rb") as f: + tect = tomllib.load(f) + for k in ("event", "eras"): + if k in tect: + raise ConfigError(f"tectonics.toml: [{k}] belongs in config/{ERAS_FILE} (editing an era must not rebuild " + f"the base world)") + epath = root / "config" / ERAS_FILE + if epath.exists(): + with open(epath, "rb") as f: + e = tomllib.load(f) + unknown = set(e) - {"event", "eras"} + if unknown: + raise ConfigError(f"{ERAS_FILE}: unknown tables {sorted(unknown)} (only [[event]] and [eras])") + tect = {**tect, **e} + _check_ranges(cfg, REQUIRED, "world.toml") + validate_tectonics(tect) + validate_revision(tect, float(cfg["planet"]["radius_km"])) + return cfg, tect + + +def params(cfg: dict, section: str, defaults: dict) -> dict: + """Module defaults overridden by world.toml [section]; unknown keys are errors (typo guard).""" + over = cfg.get(section, {}) + unknown = set(over) - set(defaults) + if unknown: + raise ConfigError(f"world.toml [{section}]: unknown keys {sorted(unknown)}") + return {**defaults, **over} diff --git a/mapgen/crust.py b/mapgen/crust.py new file mode 100644 index 0000000..345a1aa --- /dev/null +++ b/mapgen/crust.py @@ -0,0 +1,97 @@ +"""Stage `crust`: continental vs oceanic crust, age classes, oceanic age.""" +from __future__ import annotations + +import numpy as np + +from .config import params +from . import plateaus as PL +from .graph import distance_to, smooth_km +from .noise import fbm, name_seed, ridged +from .plates import CONV, DIV +from .sphere import azimuth_deg, gc_dist_km, latlon_to_xyz + +OCEANIC, CRATON, PRE_OROGEN, POST_OROGEN, RIFT, LIP, SCAR, VOLCANO = range(8) +AGE_NAMES = ["oceanic", "craton (3g era)", "pre-Lightening orogen", "post-Lightening orogen", + "rift", "Lightening basalt province", "collapse scar", "overshoot volcano"] +DEFAULTS = {"continental_threshold": 0.25, "shelf_km": 450.0, "edge_noise": 0.5, "edge_noise_freq": 4.0, "edge_noise_gain": 0.65, "orogen_width_km": 700.0, "rift_width_km": 200.0, + "pre_orogen_threshold": 0.72, "half_spreading_km_myr": 30.0, "max_ocean_age_myr": 200.0, + "scar_inner": 0.4, "scar_outer": 1.6, "scar_halfwidth_deg": 25.0} + + +def center_dist(g, center): + return gc_dist_km(g.xyz, latlon_to_xyz(*center), g.radius_km) + + +def lip_profile(g, lip, seed): + r = lip["radius_km"] * (1.0 + 0.25 * fbm(g.xyz, name_seed(seed, lip["name"]), 3, 6.0)) + x = center_dist(g, lip["center"]) / r + return np.clip((1.0 - x) / 0.25, 0.0, 1.0) + + +def _scar_azimuths(v, seed): + if "scar_azimuths" in v: + return list(v["scar_azimuths"]) + rng = np.random.default_rng(name_seed(seed, v["name"])) + return list(rng.uniform(0, 360, size=2)) + + +def patch_field(g, patch: dict, seed: int): + """[[land_patch]]: land hint (strength) over a disc whose radius varies ± 25 % · edge_noise (fractal edge).""" + d = center_dist(g, patch["center"]) + r = float(patch["radius_km"]) + out = np.zeros(g.n) + near = d < 1.3 * r + if near.any(): + w = np.clip(2.0 * fbm(g.xyz[near], name_seed(seed, patch["name"]), 6, 12.0), -1.0, 1.0) + out[near] = patch.get("strength", 1.0) * (d[near] < r * (1.0 + 0.25 * patch.get("edge_noise", 0.0) * w)) + return out + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "crust", DEFAULTS) + land, land_hint, btype = ctx.need("sk_land", "m_land_hint", "bnd_type") + rough = (fbm(g.xyz, ctx.seed + 23, 7, P["edge_noise_freq"], gain=P["edge_noise_gain"]) + if P["edge_noise"] > 0 else None) + + def continental_of(raw): + field = smooth_km(g, raw, P["shelf_km"]) + if rough is not None: # fractal margins: bays, peninsulas, offshore continental fragments + f = np.clip(field, 0.0, 1.0) + near = 4.0 * f * (1.0 - f) # strongest at the margin; none deep inland / far offshore + field = field + P["edge_noise"] * near * rough / max(float(rough.std()), 1e-12) + else: + field = np.maximum(field, (raw > 0.5).astype(float)) + cont = field > P["continental_threshold"] + for m in ctx.tect.get("microcontinent", []): + cont |= center_dist(g, m["center"]) < m["radius_km"] + return cont + + raw = np.clip(land + 0.5 * land_hint, 0.0, 1.0) + continental_base = continental_of(raw) # without land patches: its sea level is the world's + patches = sum((patch_field(g, p, ctx.seed) for p in ctx.tect.get("land_patch", [])), np.zeros(g.n)) + continental = continental_of(np.clip(raw + patches, 0.0, 1.0)) if patches.any() else continental_base.copy() + plateau_id = PL.cell_ids(g.xyz, ctx.tect.get("plateau", []), ctx.seed, g.radius_km) + continental |= plateau_id >= 0 # sunken plateaus: continental crust that never rose + d_conv = distance_to(g, btype == CONV) + d_div = distance_to(g, btype == DIV) + age = np.full(g.n, CRATON, np.int8) + age[continental & (ridged(g.xyz, ctx.seed + 21, 4, 4.0) > P["pre_orogen_threshold"])] = PRE_OROGEN + age[continental & (d_div < P["rift_width_km"])] = RIFT + age[continental & (d_conv < P["orogen_width_km"])] = POST_OROGEN + age[~continental] = OCEANIC + for lip in ctx.tect.get("lip", []): + age[lip_profile(g, lip, ctx.seed) > 0] = LIP + for v in ctx.tect.get("volcano", []): + d = center_dist(g, v["center"]) + r = v["radius_km"] + age[d < r] = VOLCANO + az = azimuth_deg(latlon_to_xyz(*v["center"]), g.xyz) + for a in _scar_azimuths(v, ctx.seed): + diff = np.abs((az - a + 180.0) % 360.0 - 180.0) + age[(diff < P["scar_halfwidth_deg"]) & (d > P["scar_inner"] * r) & (d < P["scar_outer"] * r)] = SCAR + oceanic_age = np.minimum(d_div / P["half_spreading_km_myr"], P["max_ocean_age_myr"]) + ocean_age = np.where(continental & (plateau_id < 0), 0.0, oceanic_age) # plateaus: the floor around them + return {"continental": continental, "age_class": age, "ocean_age_myr": ocean_age.astype(np.float32), + "d_conv_km": d_conv, "d_div_km": d_div, "continental_base": continental_base, "plateau_id": plateau_id, + "land_patch": patches.astype(np.float32)} diff --git a/mapgen/elevation.py b/mapgen/elevation.py new file mode 100644 index 0000000..47f0060 --- /dev/null +++ b/mapgen/elevation.py @@ -0,0 +1,197 @@ +"""Stage `elevation`: tectonic relief + Lightening provinces + hotspots + hints; sea level solved.""" +from __future__ import annotations + +import numpy as np + +from .config import params +from .crust import PRE_OROGEN, SCAR, center_dist, lip_profile, name_seed +from .fields import gravity_mod +from .graph import OCEAN_MIN_KM2, distance_to, nearest_source, ocean_mask, smooth_km +from .noise import fbm, ridged +from . import plateaus as PL +from .pipeline import StageError +from .plates import edge_convergence +from .sphere import east_north, gc_dist_km, great_circle_point, latlon_to_xyz + +OVER, SUB, COLLISION = 1, 2, 3 +DEFAULTS = { + "continental_base_m": 400.0, "continental_noise_m": 350.0, + "margin_km": 500.0, "shelf_m": -200.0, "slope_km": 250.0, + "coast_noise_m": 900.0, "coast_band_km": 600.0, "coast_noise_freq": 6.0, + "ridge_depth_m": 2500.0, "age_depth_coeff": 350.0, "abyss_m": 6500.0, + "rate_full_m_yr": 0.05, "min_rate_m_yr": 0.005, + "trench_depth_m": 4000.0, "trench_width_km": 70.0, + "arc_cont_m": 5500.0, "arc_ocean_m": 3500.0, "arc_offset_km": 200.0, "arc_width_km": 110.0, + "collision_peak_m": 9500.0, "collision_width_km": 260.0, "plateau_m": 4500.0, "plateau_km": 800.0, + "rift_depth_m": 1200.0, "rift_shoulder_m": 800.0, + "pre_orogen_m": 1200.0, "lip_m": 1500.0, "lip_step_m": 300.0, + "apron_m": 500.0, "scar_drop_m": 1500.0, + "hotspot_m": 5500.0, "hotspot_spacing_km": 150.0, "hotspot_radius_km": 70.0, + "hint_m": 2500.0, "hint_km": 80.0, "land_hint_m": 800.0, "spire_m": 2500.0, "detail_m": 250.0, + "max_land_m": 12000.0, "min_ocean_m": -11000.0, +} + + +def solve_sea_level(z, area, land_fraction): + order = np.argsort(-z) + cum = np.cumsum(area[order]) / area.sum() + k = min(int(np.searchsorted(cum, land_fraction)), len(z) - 1) + return z - z[order[k]] + + +def solve_sea_level_connected(g, z, land_fraction, min_sea_km2=OCEAN_MIN_KM2, iters=30): + """Shift z so land = everything outside the connected ocean covers land_fraction (interior pits stay land).""" + area, tot = g.area_km2, g.area_km2.sum() + if abs(area[~ocean_mask(g, z, min_sea_km2)].sum() / tot - land_fraction) <= 0.5 * area.min() / tot: + return z # already there: a flat sea floor would pull the bisection onto it + s0 = float(np.asarray(z)[np.argsort(-z)][min(int(np.searchsorted(np.cumsum(area[np.argsort(-z)]) / tot, + land_fraction)), len(z) - 1)]) + lo, hi = s0 - 3000.0, s0 + 3000.0 + for _ in range(iters): + mid = 0.5 * (lo + hi) + if area[~ocean_mask(g, z - mid, min_sea_km2)].sum() / tot > land_fraction: + lo = mid + else: + hi = mid + return z - 0.5 * (lo + hi) + + +def _smoothstep(a, b, x): + t = np.clip((x - a) / (b - a), 0.0, 1.0) + return t * t * (3 - 2 * t) + + +def roles(g, plate, vel, continental, ocean_age, plate_continental, min_rate): + conv, _ = edge_convergence(g, vel) + s, d = g.src, g.dst + e = (plate[s] != plate[d]) & (conv > min_rate) + cs, cd = continental[s[e]], continental[d[e]] + ks, kd = plate_continental[plate[s[e]]], plate_continental[plate[d[e]]] + over_s = np.where(cs & ~cd, True, np.where(~cs & cd, False, + np.where(ks & ~kd, True, np.where(~ks & kd, False, ocean_age[s[e]] < ocean_age[d[e]])))) + role_e = np.where(cs & cd, COLLISION, np.where(over_s, OVER, SUB)).astype(np.int8) + order = np.argsort(conv[e]) + role = np.zeros(g.n, np.int8) + role[s[e][order]] = role_e[order] + return role + + +def _on_plate(g, plate, role_mask, other=None): + d, src = nearest_source(g, np.flatnonzero(role_mask)) + ok = src >= 0 + same = ok & (plate == plate[np.maximum(src, 0)]) + if other is not None: + same |= ok & (plate == other[np.maximum(src, 0)]) + return np.where(same, d, np.inf), np.maximum(src, 0) + + +def hotspot_track(g, h: dict, vel, spacing_km: float): + """(unit vector, km from the active end) along a hotspot chain; the chain follows the plate's motion.""" + c = latlon_to_xyz(*h["center"]) + v = vel[g.cell_index(*h["center"])] + sp = np.linalg.norm(v) + t = v / sp if sp > 0 else east_north(c[None])[0][0] + return [(great_circle_point(c, t, k * spacing_km, g.radius_km), k * spacing_km) + for k in range(int(h["length_km"] // spacing_km) + 1)] + + +def _relief(ctx, g, P, cont): + """Heights before the sea-level solve for the continental mask `cont` → (z, role, d_over, d_sub, d_coll).""" + R, seed = g.radius_km, ctx.seed + (plate, vel, pk, brate, bother, age, ocean_age, d_div, land_hint, sk_mtn, mtn_hint) = ctx.need( + "plate", "vel", "plate_continental", "bnd_rate", "bnd_other", "age_class", "ocean_age_myr", "d_div_km", + "m_land_hint", "sk_mountains", "m_mountain_hint") + f = lambda r: np.clip(r / P["rate_full_m_yr"], 0.3, 1.0) + + # passive-margin profile: interior plateau → coastal ramp → shelf → slope → abyssal floor + d_in = distance_to(g, ~cont) if (~cont).any() else np.full(g.n, np.inf) + d_out = distance_to(g, cont) if cont.any() else np.full(g.n, np.inf) + ramp = _smoothstep(0.0, P["margin_km"], d_in) + z_cont = P["shelf_m"] + (P["continental_base_m"] - P["shelf_m"]) * ramp + z_cont = z_cont + P["continental_noise_m"] * fbm(g.xyz, seed + 31, 5, 3.0) + z_ocean = -np.minimum(P["ridge_depth_m"] + P["age_depth_coeff"] * np.sqrt(ocean_age), P["abyss_m"]) + z_ocean = P["shelf_m"] + (z_ocean - P["shelf_m"]) * _smoothstep(0.0, P["slope_km"], d_out) + z = np.where(cont, z_cont, z_ocean) + # multi-scale coastal noise: sea level cuts it → bays, headlands, drowned valleys, offshore islands + d_edge = np.where(cont, d_in, d_out) + z += P["coast_noise_m"] * np.exp(-d_edge / P["coast_band_km"]) * fbm(g.xyz, seed + 91, 7, P["coast_noise_freq"]) + + role = roles(g, plate, vel, cont, ocean_age, pk, P["min_rate_m_yr"]) + d_over, s_over = _on_plate(g, plate, role == OVER) + d_sub, s_sub = _on_plate(g, plate, role == SUB) + d_coll, s_coll = _on_plate(g, plate, role == COLLISION, other=bother) + + z += -P["trench_depth_m"] * f(brate[s_sub]) * np.exp(-(d_sub / P["trench_width_km"]) ** 2) + arc_h = np.where(cont, P["arc_cont_m"], P["arc_ocean_m"]) + z += arc_h * f(brate[s_over]) * np.exp(-((d_over - P["arc_offset_km"]) / P["arc_width_km"]) ** 2) + w = P["collision_width_km"] + z += P["collision_peak_m"] * f(brate[s_coll]) * np.exp(-(d_coll / w) ** 2) + plateau = _smoothstep(0.5 * w, 1.5 * w, d_coll) * (1 - _smoothstep(0.7 * P["plateau_km"], P["plateau_km"], d_coll)) + z += np.where(cont, P["plateau_m"] * f(brate[s_coll]) * plateau, 0.0) + rift = -P["rift_depth_m"] * np.exp(-(d_div / 60.0) ** 2) + P["rift_shoulder_m"] * np.exp(-((d_div - 120.0) / 60.0) ** 2) + z += np.where(cont, rift, 0.0) + + z += np.where(age == PRE_OROGEN, P["pre_orogen_m"] * ridged(g.xyz, seed + 21, 4, 4.0), 0.0) + for lip in ctx.tect.get("lip", []): + prof = lip_profile(g, lip, seed) + z += np.floor(P["lip_m"] * prof / P["lip_step_m"]) * P["lip_step_m"] + for v in ctx.tect.get("volcano", []): + d = center_dist(g, v["center"]) + r, H = v["radius_km"], v.get("height_m", 7000.0) + z += H * np.clip(1 - d / r, 0, 1) ** 1.5 - 0.35 * H * np.exp(-(d / (0.12 * r)) ** 2) + z += P["apron_m"] * np.exp(-((d - r) / (0.3 * r)) ** 2) * (0.5 + 0.5 * fbm(g.xyz, name_seed(seed, v["name"]), 3, 20.0)) + z -= np.where(age == SCAR, P["scar_drop_m"], 0.0) + + for h in ctx.tect.get("hotspot", []): + peaks = np.zeros(g.n) + for p, s in hotspot_track(g, h, vel, P["hotspot_spacing_km"]): + dk = gc_dist_km(g.xyz, p, R) + hk = P["hotspot_m"] * (1 - s / h["length_km"]) + peaks = np.maximum(peaks, hk * np.clip(1 - dk / P["hotspot_radius_km"], 0, 1) ** 1.2) + z += peaks + + hint = smooth_km(g, np.clip(sk_mtn + mtn_hint, 0, 1), P["hint_km"]) + z += np.where(cont, P["hint_m"] * hint, 0.0) + z += P["land_hint_m"] * land_hint + z += P["detail_m"] * fbm(g.xyz, seed + 41, 6, 12.0) + return z, role, d_over, d_sub, d_coll + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "elevation", DEFAULTS) + continental, sk_land, land_hint, m_grav, m_lock = ctx.need("continental", "sk_land", "m_land_hint", + "m_gravity_zones", "m_lock") + plateau_id = np.asarray(ctx.data.get("plateau_id", np.full(g.n, -1))) + cont = np.asarray(continental) & (plateau_id < 0) # plateaus: ocean until their surfaces are set (§3) + base = np.asarray(ctx.data.get("continental_base", continental)) & (plateau_id < 0) + z, role, d_over, d_sub, d_coll = _relief(ctx, g, P, cont) + + target = ctx.cfg["build"]["land_fraction"] + cont_area = g.area_km2[base].sum() / g.area_km2.sum() + if target > cont_area + 0.005: + raise StageError(f"land_fraction {target} exceeds continental crust area {cont_area:.3f}: sea level would " + f"lift ocean ridges into land; lower [build] land_fraction or raise [crust] shelf_km") + min_sea = ctx.cfg.get("erosion", {}).get("min_sea_km2", OCEAN_MIN_KM2) + sea_ref = np.zeros(g.n, bool) # sea in the world without land patches + if np.array_equal(base, cont): + z = solve_sea_level_connected(g, z, target, min_sea) + else: # land patches (the eastern continent made whole) must not move every other coast: the sea level of the + z_raw = _relief(ctx, g, P, base)[0] # world without them, applied to the world with them + z_ref = solve_sea_level_connected(g, z_raw, target, min_sea) + z = z - float(np.mean(z_raw - z_ref)) + sea_ref = ocean_mask(g, z_ref, min_sea) + patch = smooth_km(g, np.asarray(ctx.data.get("land_patch", np.zeros(g.n)), dtype=np.float64), P["hint_km"]) + z += P["land_hint_m"] * patch # a patch is a land hint in height too (bare crust + # sits below the world's sea level: land) + mult = np.clip(1.0 / gravity_mod(ctx.cfg, m_grav), 1.0, 3.0) + spires = P["spire_m"] * (mult - 1) * ridged(g.xyz, ctx.seed + 71, 4, 40.0) + z = np.where(z > 0, z * mult + spires, z) + target_land = (sk_land + 0.5 * land_hint) > 0.5 + lock = m_lock > 0.5 + z = np.where(lock & target_land, np.maximum(z, 50.0), np.where(lock & ~target_land, np.minimum(z, -50.0), z)) + z = PL.apply(g, z, ctx.tect.get("plateau", []), plateau_id, ctx.seed) # sunken plateaus (after the solve) + z = np.clip(z, P["min_ocean_m"], P["max_land_m"] * mult) + added = ~ocean_mask(g, z, min_sea) & (sea_ref | (plateau_id >= 0)) # land the sea-level solve never saw + return {"elevation_m": z.astype(np.float32), "role": role, "d_over_km": d_over, "d_sub_km": d_sub, + "d_coll_km": d_coll, "land_added": added} diff --git a/mapgen/environment.py b/mapgen/environment.py new file mode 100644 index 0000000..fbd3c8e --- /dev/null +++ b/mapgen/environment.py @@ -0,0 +1,158 @@ +"""Stage `environment`: Holdridge life zone + seasonality tag + lithology + landform + ground.""" +from __future__ import annotations + +import numpy as np + +from .config import params +from .crust import CRATON, LIP, POST_OROGEN, PRE_OROGEN, RIFT, SCAR, VOLCANO +from .graph import nbr_max, nbr_min, steepest_receivers +from . import minerals as MN +from .noise import fbm + +REGIONS = ["polar", "subpolar", "boreal", "cool temperate", "warm temperate", "subtropical", "tropical"] +ZONES = { + "polar": ["desert"], + "subpolar": ["dry tundra", "moist tundra", "wet tundra", "rain tundra"], + "boreal": ["desert", "dry scrub", "moist forest", "wet forest", "rain forest"], + "cool temperate": ["desert", "desert scrub", "steppe", "moist forest", "wet forest", "rain forest"], + "warm temperate": ["desert", "desert scrub", "thorn steppe", "dry forest", "moist forest", "wet forest", + "rain forest"], + "subtropical": ["desert", "desert scrub", "thorn woodland", "dry forest", "moist forest", "wet forest", + "rain forest"], + "tropical": ["desert", "desert scrub", "thorn woodland", "very dry forest", "dry forest", "moist forest", + "wet forest", "rain forest"], +} +THRESH = { # annual precipitation (mm) upper bounds between successive zones + "polar": [], "subpolar": [125, 250, 500], "boreal": [125, 250, 500, 1000], + "cool temperate": [125, 250, 500, 1000, 2000], "warm temperate": [125, 250, 500, 1000, 2000, 4000], + "subtropical": [125, 250, 500, 1000, 2000, 4000], "tropical": [125, 250, 500, 1000, 2000, 4000, 8000], +} +HOLDRIDGE_NAMES = [f"{r} {z}" for r in REGIONS for z in ZONES[r]] + +SEAS_F, SEAS_S, SEAS_W, SEAS_M = 0, 1, 2, 3 +SEASONALITY_NAMES = ["even rainfall", "dry summer", "dry winter / summer monsoon", "tropical monsoon"] + +LI_OCEAN, LI_GRANITE, LI_LIMESTONE, LI_BASALT, LI_ANDESITE, LI_SANDSTONE, LI_METAMORPHIC = range(7) +LITHOLOGY_NAMES = ["oceanic basalt", "granite/gneiss", "limestone", "basalt", "andesite", "sandstone/shale", + "metamorphic"] + +(LF_OCEAN, LF_PLAIN, LF_HILLS, LF_MOUNTAINS, LF_PLATEAU, LF_RIFT, LF_ESCARPMENT, LF_ARC, LF_MASSIF, + LF_BASALT_PLATEAU, LF_DUNES, LF_BADLANDS) = range(12) +LANDFORM_NAMES = ["ocean", "plain", "hills", "mountains", "plateau", "rift valley", "escarpment", + "volcanic arc", "volcanic massif", "basalt plateau", "dunes", "badlands"] + +(GR_NONE, GR_WETLAND, GR_BOG, GR_FLOODPLAIN, GR_DELTA, GR_MANGROVE, GR_SALT_FLAT, GR_PERMAFROST, + GR_KARST) = range(9) +GROUND_NAMES = ["—", "wetland", "bog", "floodplain", "delta", "mangrove", "salt flat", "permafrost", "karst"] + +LEGENDS = {"deposits": MN.LEGEND, "holdridge": HOLDRIDGE_NAMES, "seasonality": SEASONALITY_NAMES, "landform": LANDFORM_NAMES, + "lithology": LITHOLOGY_NAMES, "ground": GROUND_NAMES} + +DEFAULTS = {"relief_km": 100.0, "coastal_km": 60.0, "dry_ratio_s": 3.0, "dry_ratio_w": 4.0, "monsoon_min_mm": 1500.0, "flat_m_per_km": 1.0, + "wet_min_mm": 700.0, "karst_min_mm": 800.0, "small_river_km3_yr": 5.0, "life": True} + + +def holdridge(bio, p, tmin): + region = np.select([bio < 1.5, bio < 3, bio < 6, bio < 12, (bio < 24) & (tmin < 0), bio < 24], + [0, 1, 2, 3, 4, 5], default=6).astype(np.int8) + zone = np.zeros(len(bio), np.int16) + offset = 0 + for ri, r in enumerate(REGIONS): + sel = region == ri + zone[sel] = offset + np.searchsorted(THRESH[r], p[sel], side="right") + offset += len(ZONES[r]) + return zone, region + + +def seasonality(p_jun, p_dec, lat, tmin, p_ann, P): + summer = np.where(lat >= 0, p_jun, p_dec) + winter = np.where(lat >= 0, p_dec, p_jun) + tag = np.full(len(lat), SEAS_F, np.int8) + tag[summer < winter / P["dry_ratio_s"]] = SEAS_S + w = winter < summer / P["dry_ratio_w"] + tag[w] = SEAS_W + tag[w & (tmin >= 18.0) & (p_ann >= P["monsoon_min_mm"])] = SEAS_M + return tag + + +def lithology(g, z, continental, age, d_over, seed): + nz = fbm(g.xyz, seed + 51, 4, 5.0) + lit = np.full(g.n, LI_GRANITE, np.int8) + lit[continental & (z < 600) & (nz > 0.1)] = LI_SANDSTONE + lit[continental & (z < 300) & (nz < -0.1)] = LI_LIMESTONE + lit[np.isin(age, [POST_OROGEN, PRE_OROGEN]) & (z > 2000)] = LI_METAMORPHIC + lit[~continental] = LI_OCEAN + lit[(d_over > 100) & (d_over < 320)] = LI_ANDESITE + lit[np.isin(age, [LIP, VOLCANO, SCAR])] = LI_BASALT + return lit + + +def landform(g, z, age, lit, d_over, p_ann, ocean=None, relief_km=100.0): + hi, lo = np.asarray(z, dtype=np.float64), np.asarray(z, dtype=np.float64) + for _ in range(max(1, int(round(relief_km / g.spacing_km)))): # window ≈ relief_km at any resolution + hi, lo = nbr_max(g, hi), nbr_min(g, lo) + relief = hi - lo + lf = np.full(g.n, LF_PLAIN, np.int8) + lf[relief >= 300] = LF_HILLS + lf[(z > 1000) & (relief < 500)] = LF_PLATEAU + lf[(relief >= 1500) | (z > 2500)] = LF_MOUNTAINS + lf[(p_ann < 250) & (relief < 300) & (lit == LI_SANDSTONE)] = LF_DUNES + lf[(p_ann >= 250) & (p_ann < 500) & (relief >= 200) & (relief < 800) & (lit == LI_SANDSTONE)] = LF_BADLANDS + lf[age == RIFT] = LF_RIFT + lf[age == LIP] = LF_BASALT_PLATEAU + ocean = z <= 0 if ocean is None else ocean + lf[(d_over > 100) & (d_over < 320) & ~ocean] = LF_ARC + lf[age == VOLCANO] = LF_MASSIF + lf[age == SCAR] = LF_ESCARPMENT + lf[ocean] = LF_OCEAN + return lf, relief + + +def ground(g, z, lit, t_mean, t_min, p_ann, dist_ocean, river, strahler, discharge, lake, salt_flat, P, ocean=None): + land = ~(z <= 0 if ocean is None else ocean) + _, slope, _ = steepest_receivers(g, z) + flat = slope < P["flat_m_per_km"] + coastal = land & (dist_ocean < P["coastal_km"]) + wet = land & flat & ~lake & (p_ann > P["wet_min_mm"]) & (discharge > P["small_river_km3_yr"]) + gr = np.zeros(g.n, np.int8) + gr[land & (t_mean < -2.0)] = GR_PERMAFROST + gr[land & (lit == LI_LIMESTONE) & (p_ann > P["karst_min_mm"])] = GR_KARST + gr[wet & (t_mean < 5.0)] = GR_BOG + gr[wet & (t_mean >= 5.0)] = GR_WETLAND + gr[river & (strahler >= 3) & flat] = GR_FLOODPLAIN + gr[coastal & flat & (t_min >= 20.0) & (p_ann > 1000.0)] = GR_MANGROVE + gr[river & (strahler >= 4) & coastal] = GR_DELTA + gr[salt_flat] = GR_SALT_FLAT + if not P["life"]: + gr[np.isin(gr, [GR_WETLAND, GR_BOG, GR_MANGROVE])] = GR_NONE + return gr + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "environment", DEFAULTS) + (z, cont, age, d_over, t_mean, t_min, p_ann, p_jun, p_dec, bio, dist_ocean, river, strahler, q, lake, + salt) = ctx.need("elevation_eroded_m", "continental", "age_class", "d_over_km", "T_mean", "T_min", "P_ann", + "P_jun", "P_dec", "biotemp", "dist_ocean_km", "river", "strahler", "discharge_km3_yr", + "lake", "salt_flat") + z = z.astype(np.float64) + bed = z + if "lake_level_m" in ctx.data: # the ground's shape: lakes at their water, not their beds + lev = np.asarray(ctx.data["lake_level_m"], dtype=np.float64) + z = np.where(np.isfinite(lev), np.maximum(z, lev), z) + ocean = np.asarray(ctx.data["ocean"]) if "ocean" in ctx.data else z <= 0 + zone, region = holdridge(bio, p_ann, t_min) + lit = lithology(g, z, cont, age, d_over, ctx.seed) + lf, relief = landform(g, z, age, lit, d_over, p_ann, ocean, P["relief_km"]) + gr = ground(g, z, lit, t_mean, t_min, p_ann, dist_ocean, river, strahler, q, lake, salt, P, ocean) + coal = (lit == LI_SANDSTONE) & (p_ann > 600) & ~ocean & bool(P["life"]) + get = lambda k, v: np.asarray(ctx.data[k]) if k in ctx.data else np.full(g.n, v) + deposits, main = MN.place(g, {"land": ~ocean & ~lake, "z": bed, "lit": lit, "age": age, "lf": lf, "relief": relief, + "t_mean": t_mean, "p_ann": p_ann, "d_coll": get("d_coll_km", np.inf), + "salt": salt, "endo": get("endorheic", False), "dist_ocean": dist_ocean, + "coal": coal}, ctx.seed) + iron = np.uint32(MN.BIT["bog_iron"] | MN.BIT["ironstone"] | MN.BIT["iron_high"]) + return {"deposits": deposits, "deposit_main": main,"holdridge": zone, "hold_region": region, + "seasonality": seasonality(p_jun, p_dec, g.lat, t_min, p_ann, P), + "landform": lf, "relief_m": relief, "lithology": lit, "ground": gr, + "coal_potential": coal, "iron_potential": (deposits & iron) > 0} diff --git a/mapgen/eras.py b/mapgen/eras.py new file mode 100644 index 0000000..9b0dbc4 --- /dev/null +++ b/mapgen/eras.py @@ -0,0 +1,219 @@ +"""Eras: the base world plus local events. build_era applies an era's events +to the built base, re-runs the stages they affect inside an influence mask (changed cells + mask_km) and keeps every +cell outside it exactly as the base. Results: out/r<res>/eras/<era>/ (as out/r<res>/), previews in +previews/r<res>/eras/<era>/. An era's cache key is the base world's inputs key plus its events (config/eras.toml). Each era's events happen at the end of the era before it; its `years` erode what the events so far +changed; zone events follow the land at the start of their own era.""" +from __future__ import annotations + +import hashlib +import importlib +import json +import re +import shutil +from pathlib import Path + +import numpy as np + +from . import config as C, environment as EN, events as EV, fields as FL, pipeline as P +from .config import params +from .erosion import DEFAULTS as EROSION_DEFAULTS, K_MULT, erode +from .graph import distance_to, ocean_mask + +DEFAULTS = {"mask_km": 1500.0, "climate_blend_km": 300.0, "erosion_steps": 1, "years": 10000.0} +RERUN = ("climate", "hydrology", "seabed", "environment", "ice", "fields") + + +def era_dir(root: Path, res: int, name: str) -> Path: + return Path(root) / "out" / f"r{res}" / "eras" / name + + +def steps(tect: dict, name: str) -> list: + """[(era, its own events, years)] for every era up to and including `name` ([eras] order): an era's events happen + at the end of the era before it; `years` is the time from them to the era's map (None: the default).""" + C.era_events(tect, name) # an unknown era: ConfigError + eras = tect.get("eras") or {} + order = eras["order"] + by_name = {e["name"]: e for e in tect.get("event", [])} + return [(n, [by_name[e] for e in eras.get(n, {}).get("events", [])], eras.get(n, {}).get("years")) + for n in order[: order.index(name) + 1]] + + +def fingerprint(base_key: str, steps_: list) -> str: + content = [[own, yrs] for _, own, yrs in steps_] + return hashlib.sha256((base_key + json.dumps(content, sort_keys=True)).encode()).hexdigest()[:16] + + +def _same(a, b) -> bool: + """Equal arrays; NaN equals NaN in float arrays (a field may hold NaN where it has no value).""" + a, b = np.asarray(a), np.asarray(b) + return np.array_equal(a, b, equal_nan=a.dtype.kind in "fc" and b.dtype.kind in "fc") + + +def merge_lake_ids(base, new, mask): + """Lake ids of an era: the base's outside the mask; inside it a lake the era left as it was keeps its base id and + every other lake gets a fresh id after the base's, so one id never names two lakes.""" + base, new, mask = np.asarray(base), np.asarray(new), np.asarray(mask, bool) + out = base.copy() + nxt = max(int(base.max()) + 1, 0) if len(base) else 0 + counts = np.bincount(base[base >= 0]) if (base >= 0).any() else np.zeros(0, np.int64) + ks = np.unique(new[mask & (new >= 0)]) + idx = np.flatnonzero(np.isin(new, ks)) + idx = idx[np.argsort(new[idx], kind="stable")] + starts = np.searchsorted(new[idx], ks) + ends = np.append(starts[1:], len(idx)) + for s, e in zip(starts.tolist(), ends.tolist()): + cells = idx[s:e] + b = base[cells] + kept = b[0] >= 0 and bool(np.all(b == b[0])) and counts[b[0]] == len(cells) + into = cells[mask[cells]] + if kept: + out[into] = b[0] + else: + out[into] = nxt + nxt += 1 + out[mask & (new < 0)] = -1 + return out.astype(base.dtype) + + +def zone_key(name: str) -> str: + return "zone_" + re.sub(r"[^A-Za-z0-9]", "_", name) + + +def _heights(g, z, events, min_sea, log, name, uncut=lambda z: z): + for ev in events: # heights first, in era order + if ev["kind"] == "disintegrate": + z, z_ref = EV.disintegrate(g.xyz, z, ~ocean_mask(g, uncut(z), min_sea), ev, g.radius_km) + if z_ref is None: + log(f"era {name}: warning: {ev['name']} has no land in its rim ring (0.9–1.0 × radius): " + f"its bowl hangs from 0 m") + else: + log(f"era {name}: {ev['name']} (ground at its rim {z_ref:.0f} m)") + elif ev["kind"] == "volcano": + z = EV.volcano(g.xyz, z, ev, g.radius_km) + return z + + +def build_era(root: Path, res: int, name: str, log=print, base=None, low_memory: bool = False) -> Path: + root = Path(root) + cfg, tect = C.load(root) + events = C.era_events(tect, name) + if not events: + log(f"era {name}: the base world (out/r{res})") + return root / "out" / f"r{res}" + out = era_dir(root, res, name) + base_key = P.inputs_key(root, res) + plan = steps(tect, name) + key = fingerprint(base_key, plan) + label = tect["eras"][name].get("label", name) + meta_path = out / "cells_meta.json" + if meta_path.exists() and (out / "cells.npz").exists(): + meta = json.loads(meta_path.read_text()) + era = meta.get("era", {}) + if era.get("fingerprint") == key: + if (era.get("label"), era.get("name")) != (label, name): # a new label or name needs no rebuild + era["label"], era["name"] = label, name + meta_path.write_text(json.dumps(meta, indent=1)) + log(f"era {name}: cached") + return out + ctx = base if base is not None else P.build(root, res, stop="fields", log=lambda m: None, low_memory=low_memory) + low_memory = low_memory or ctx.low_memory + g, d = ctx.grid, ctx.data + E = params(ctx.cfg, "eras", DEFAULTS) + EP = {**params(ctx.cfg, "erosion", EROSION_DEFAULTS), "steps": E["erosion_steps"]} + kmult = np.vectorize(K_MULT.get)(np.asarray(d["age_class"])).astype(np.float64) + z0 = np.asarray(d["elevation_eroded_m"], dtype=np.float64) + z, ocean = z0.copy(), np.asarray(d["ocean"]) + water = np.asarray(d.get("open_water", ocean)) + cut = np.asarray(d.get("lake_cut_m", np.zeros(g.n)), dtype=np.float64) + uncut = lambda zz: np.where(zz == z0, zz + cut, zz) # sea masks: the ground before lake beds were carved + zones, lands = [], [] + for _, own, yrs in plan: # era by era: its events at the end of the era before + for ev in own: # zones follow the land before their own era's events + if ev["kind"] == "zone": + zones.append(ev) + lands.append(EV.landmass(g, ~water, ev["seed"], ev["name"]) if ev["shape"] == "landmass" else None) + z = _heights(g, z, own, EP["min_sea_km2"], log, name, uncut) + hit = z != z0 + if hit.any(): # the changed ground weathers for the era's years + years = float(E["years"] if yrs is None else yrs) + ze = erode(g, z, kmult, {**EP, "dt_myr": years / 1.0e6 / E["erosion_steps"]}) + z = np.where(hit, np.maximum(np.minimum(z, ze), -11000.0), z) + ocean, water = ocean_mask(g, uncut(z), np.inf), ocean_mask(g, uncut(z), EP["min_sea_km2"]) + hit = z != z0 + weights = [EV.zone_weight(g.xyz, ev, g.radius_km, ctx.seed, None if land is None else g.xyz[land]) + for ev, land in zip(zones, lands)] + changed = hit | (ocean != np.asarray(d["ocean"])) | (water != np.asarray(d.get("open_water", d["ocean"]))) + for w in weights: + changed |= w > 0 + mask = distance_to(g, changed) <= E["mask_km"] if changed.any() else np.zeros(g.n, bool) + + ec = P.Ctx(ctx.root, ctx.cfg, ctx.tect, ctx.res, dict(d), g, low_memory=low_memory) + ec.data.update({"elevation_eroded_m": z.astype(np.asarray(d["elevation_eroded_m"]).dtype), "ocean": ocean, + "open_water": water, "lake_carve": hit}) # lake beds: geology, carved where the ground moved + new_keys = {} + for s in RERUN: + o = importlib.import_module(f"mapgen.{s}").run(ec) + ec.data.update(o) + new_keys[s] = list(o) + blend = np.clip(distance_to(g, ~mask) / E["climate_blend_km"], 0.0, 1.0) if (~mask).any() else np.ones(g.n) + shape = lambda a, v: v.reshape(-1, *([1] * (a.ndim - 1))) + merged = dict(d) + for s, keys in new_keys.items(): + for k in keys: + nv, bv = np.asarray(ec.data[k]), np.asarray(d[k]) + if s == "climate" and nv.dtype.kind == "f": # solved on the whole world, blended in over 300 km + nv = (bv + shape(nv, blend) * (nv - bv)).astype(nv.dtype) + merged[k] = np.where(shape(nv, mask), nv, bv) + for k in ("deposits", "deposit_main", "iron_potential"): # ore is the base's geology: events move ground, the + if k in d: # rank-by-share placement would shift it world-wide + merged[k] = d[k] + if "lake_id" in merged: # one id never names two lakes + merged["lake_id"] = merge_lake_ids(d["lake_id"], ec.data["lake_id"], mask) + merged["elevation_eroded_m"] = np.where(mask, ec.data["elevation_eroded_m"], d["elevation_eroded_m"]) # events and + merged["ocean"] = ocean # lake beds: inside only + if "open_water" in d: + merged["open_water"] = water + H = float(ctx.cfg["planet"]["scale_height_km"]) * 1000.0 + for ev, w, land in zip(zones, weights, lands): # zone events write their fields last + merged.update(EV.apply_zone(merged, w, ev, merged["z_surface_m"], H)) + merged[zone_key(ev["name"])] = w.astype(np.float32) + if land is not None: + merged[zone_key(ev["name"]) + "_land"] = land + if zones: # plants settle into the zone's air and gravity + pl = ctx.cfg["planet"] + merged["plant_height_x"] = FL.plant_height(pl["gravity_g"], merged["gravity_g"]).astype( + np.asarray(d["plant_height_x"]).dtype) + sea_p = np.asarray(merged["pressure_bar"], dtype=np.float64) / FL.pressure(1.0, merged["z_surface_m"], H) + wet = np.sqrt(sea_p / pl["sea_level_pressure_bar"]) # more CO₂ per breath: less water lost per growth + hz, hr = EN.holdridge(merged["biotemp"], np.asarray(merged["P_ann"], dtype=np.float64) * wet, merged["T_min"]) + inside = np.logical_or.reduce([w > 0 for w in weights]) + merged["holdridge"] = np.where(inside, hz, merged["holdridge"]).astype(np.asarray(d["holdridge"]).dtype) + merged["hold_region"] = np.where(inside, hr, merged["hold_region"]).astype(np.asarray(d["hold_region"]).dtype) + merged["era_mask"] = mask + for k, a in d.items(): # guard: nothing outside the mask may differ + a, b = np.asarray(a), np.asarray(merged[k]) + if a.shape[:1] == (g.n,) and not _same(a[~mask], b[~mask]): + raise RuntimeError(f"era {name}: {k} changed outside the influence mask") + + rc = P.Ctx(ctx.root, ctx.cfg, ctx.tect, ctx.res, merged, g, out_dir=out, + preview_dir=root / "previews" / f"r{res}" / "eras" / name, low_memory=low_memory) + importlib.import_module("mapgen.render").run(rc) + meta = json.loads(meta_path.read_text()) + meta["era"] = {"name": name, "label": label, "events": events, "fingerprint": key, "base_key": base_key, + "steps": [[n, [e["name"] for e in own], yrs] for n, own, yrs in plan], + "mask_km": E["mask_km"], "mask_cells": int(mask.sum())} + meta_path.write_text(json.dumps(meta, indent=1)) + log(f"era {name}: {int(mask.sum())} of {g.n} cells inside the influence mask") + return out + + +def build_all(root: Path, res: int, log=print, base=None, low_memory: bool = False) -> list: + """Every configured era with events (the base era is the base build); era folders no longer configured go.""" + _, tect = C.load(Path(root)) + names = [n for n in (tect.get("eras") or {}).get("order", []) if C.era_events(tect, n)] + built = [build_era(root, res, n, log, base, low_memory) for n in names] + top = Path(root) / "out" / f"r{res}" / "eras" + for p in (top.iterdir() if top.is_dir() else []): + if p.is_dir() and p.name not in names: + shutil.rmtree(p, ignore_errors=True) + return built diff --git a/mapgen/erosion.py b/mapgen/erosion.py new file mode 100644 index 0000000..5384988 --- /dev/null +++ b/mapgen/erosion.py @@ -0,0 +1,54 @@ +"""Stage `erosion`: implicit stream-power incision (Braun & Willett 2013) + km-scale hillslope smoothing.""" +from __future__ import annotations + +import numpy as np + +from .config import params +from .crust import CRATON, LIP, OCEANIC, POST_OROGEN, PRE_OROGEN, RIFT, SCAR, VOLCANO +from .elevation import solve_sea_level_connected +from .graph import (accumulate, drop_smooth_cache, ocean_mask, priority_flood, receiver_levels, smooth_km, + steepest_receivers) + +DEFAULTS = {"min_sea_km2": 5.0e6, "steps": 12, "dt_myr": 5.0, "k": 0.02, "m": 0.5, "hillslope_km": 25.0, "sediment_fill": 0.15} +K_MULT = {OCEANIC: 0.0, CRATON: 1.5, PRE_OROGEN: 1.2, POST_OROGEN: 1.0, RIFT: 1.0, LIP: 0.6, SCAR: 0.9, VOLCANO: 0.8} + + +def erode(g, z, kmult, P): + z = np.asarray(z, dtype=np.float64).copy() + ar = np.arange(g.n) + for _ in range(int(P["steps"])): + ocean = ocean_mask(g, z, P["min_sea_km2"]) + zb = np.where(ocean, 0.0, z) # base level = sea level, not the seafloor + zf = priority_flood(g, zb, ocean) + zb = np.where(ocean, zb, zb + P["sediment_fill"] * (zf - zb)) # sediment infills closed basins + recv, _, dist = steepest_receivers(g, zf) + recv = np.where(ocean, ar, recv) + levels = receiver_levels(recv) + area = accumulate(recv, levels, g.area_km2) + F = P["k"] * kmult * P["dt_myr"] * area ** P["m"] / np.where(np.isfinite(dist), dist, 1.0) + F[ocean] = 0.0 + znew = zb.copy() + for lv in levels[1:]: + r = recv[lv] + znew[lv] = np.minimum(zb[lv], (zb[lv] + F[lv] * znew[r]) / (1.0 + F[lv])) + znew = smooth_km(g, znew, P["hillslope_km"], keep=True) # hillslope/sub-grid smoothing, fixed length in km + z = np.where(ocean, z, znew) + drop_smooth_cache(g) + return z + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "erosion", DEFAULTS) + z0, age = ctx.need("elevation_m", "age_class") + kmult = np.vectorize(K_MULT.get)(age).astype(np.float64) + z0 = np.asarray(z0, dtype=np.float64) + added = np.asarray(ctx.data.get("land_added", np.zeros(g.n, bool)), bool) # land the sea-level solve didn't + target = ctx.cfg["build"]["land_fraction"] + g.area_km2[added & (z0 > 0)].sum() / g.area_km2.sum() # see + z = erode(g, z0, kmult, P) + z = solve_sea_level_connected(g, z, target, P["min_sea_km2"]) + z = np.maximum(z, -11000.0) + # the sea is the connected world ocean; a separate basin ≥ min_sea_km2 shaped the coasts above as sea and stays + # open water for the climate (open_water), but its water is a lake (hydrology), not sea + return {"elevation_eroded_m": z.astype(np.float32), "ocean": ocean_mask(g, z, np.inf), + "open_water": ocean_mask(g, z, P["min_sea_km2"])} diff --git a/mapgen/events.py b/mapgen/events.py new file mode 100644 index 0000000..cba874e --- /dev/null +++ b/mapgen/events.py @@ -0,0 +1,87 @@ +"""Era events: pure functions of positions and the event's config, so the era +builder (world cells) and the refinement (fine cells) apply the same event.""" +from __future__ import annotations + +import numpy as np +from scipy.spatial import cKDTree + +from .fields import pressure +from .graph import components +from .noise import fbm, name_seed +from .pipeline import StageError +from .sphere import gc_dist_km, latlon_to_xyz +from .zones import smootherstep + + +def disintegrate(xyz, z, land, ev: dict, radius_km: float): + """(new heights, z_ref): inside radius_km of the centre min(old, z_ref − depth · sqrt(1 − (d/r)²)); z_ref = mean + ground height of the land in the ring 0.9–1.0 r before the event, None when no land lies in that ring (the bowl + then hangs from 0 m). Material vanishes: no ejecta, no rim.""" + z = np.asarray(z, dtype=np.float64) + d = gc_dist_km(np.asarray(xyz, dtype=np.float64), latlon_to_xyz(*ev["center"]), radius_km) + r = float(ev["radius_km"]) + ring = np.asarray(land, bool) & (d >= 0.9 * r) & (d <= r) + z_ref = float(np.mean(z[ring])) if ring.any() else None # None: no land on its rim (the bowl hangs from 0 m) + bowl = (0.0 if z_ref is None else z_ref) - ev["depth_m"] * np.sqrt(np.clip(1.0 - (d / r) ** 2, 0.0, 1.0)) + return np.where(d < r, np.minimum(z, bowl), z), z_ref + + +def volcano(xyz, z, ev: dict, radius_km: float): + """Heights after a volcano (only raises): cone, shield or caldera; peak_m above the old ground (peak_mode "above") + or as an absolute height ("absolute").""" + z = np.asarray(z, dtype=np.float64) + d = gc_dist_km(np.asarray(xyz, dtype=np.float64), latlon_to_xyz(*ev["center"]), radius_km) + x = np.clip(d / ev["radius_km"], 0.0, 1.0) + prof = {"cone": (1.0 - x) ** 1.5, "shield": (1.0 - x * x) ** 2, + "caldera": np.maximum((1.0 - x) ** 1.5 - 0.6 * np.exp(-(x / 0.15) ** 2), 0.0)}[ev.get("shape", "cone")] + peak = float(ev["peak_m"]) + new = z + peak * prof if ev.get("peak_mode", "above") == "above" else z + np.maximum(peak - z, 0.0) * prof + return np.maximum(z, new) + + +def landmass(g, land, seed_latlon, event_name: str): + """The connected land (grid cells) containing the seed point.""" + lab = components(g, land) + i = g.cell_index(*seed_latlon) + if lab[i] < 0: + land = np.flatnonzero(np.asarray(land, bool)) + hint = "" + if len(land): + j = land[np.argmax(g.xyz[land] @ g.xyz[i])] + hint = f"; the nearest land is at [{g.lat[j]:.1f}, {g.lon[j]:.1f}]" + raise StageError(f"event {event_name}: its seed {list(seed_latlon)} lies in the sea in this era{hint} " + f"(move the seed in config/eras.toml)") + return lab == lab[i] + + +def zone_weight(xyz, ev: dict, radius_km: float, seed: int, land_xyz=None): + """0–1 per point: 1 inside, 0 outside; a linear ramp edge_km wide at the edge (sharp), or the smooth profile. + circle: centre + radius_km. landmass: within reach of the landmass cells land_xyz (unit vectors), the reach + varying between reach_km by low-frequency noise.""" + xyz = np.asarray(xyz, dtype=np.float64) + if ev["shape"] == "circle": + d = gc_dist_km(xyz, latlon_to_xyz(*ev["center"]), radius_km) + reach = np.full(len(xyz), float(ev["radius_km"])) + else: + chord, _ = cKDTree(np.asarray(land_xyz, dtype=np.float64)).query(xyz) + d = 2.0 * np.arcsin(np.clip(chord / 2.0, 0.0, 1.0)) * radius_km + lo, hi = ev["reach_km"] + reach = lo + (hi - lo) * np.clip(0.5 + 1.5 * fbm(xyz, name_seed(seed, ev["name"]), 3, 1.5), 0.0, 1.0) + if "edge_km" in ev: + return np.clip((reach - d) / ev["edge_km"], 0.0, 1.0) + return 1.0 - smootherstep(d / reach) + + +def apply_zone(fields: dict, w, ev: dict, z_surface_m, scale_height_m: float) -> dict: + """Fields inside a zone event: absolute values blended by w (gravity_g, o2_fraction, fire_reactivity; pressure_bar + is the sea-level value, falling with height as elsewhere); po2_bar follows. Never height. Cells with w = 0 keep + their values exactly.""" + w = np.asarray(w, dtype=np.float64) + decay = pressure(1.0, z_surface_m, scale_height_m) + new = {} + for k, v in ev["fields"].items(): + base = np.asarray(fields[k], dtype=np.float64) + new[k] = (base / decay * (1.0 - w) + v * w) * decay if k == "pressure_bar" else base * (1.0 - w) + v * w + o2 = new.get("o2_fraction", np.asarray(fields["o2_fraction"], dtype=np.float64)) + new["po2_bar"] = o2 * new.get("pressure_bar", np.asarray(fields["pressure_bar"], dtype=np.float64)) + return {k: np.where(w > 0, v, fields[k]).astype(np.asarray(fields[k]).dtype) for k, v in new.items()} diff --git a/mapgen/fields.py b/mapgen/fields.py new file mode 100644 index 0000000..bf83d3f --- /dev/null +++ b/mapgen/fields.py @@ -0,0 +1,36 @@ +"""Stage `fields` and shared zone modifiers.""" +from __future__ import annotations + +import numpy as np + +from .config import MASK_DEFAULTS, params + + +def gravity_mod(cfg, m_gravity): + M = params(cfg, "masks", MASK_DEFAULTS) + return np.clip(1.0 + M["gravity_range"] * np.asarray(m_gravity, dtype=np.float64), 0.05, None) + + +def o2_mod(cfg, m_o2): + M = params(cfg, "masks", MASK_DEFAULTS) + return np.clip(1.0 + M["o2_range"] * np.asarray(m_o2, dtype=np.float64), 0.0, None) + + +def pressure(p0_bar, z_m, scale_height_m): + """Air pressure at ground height z (bar): sea-level pressure p0 falling with height (below 0 m: p0).""" + return p0_bar * np.exp(-np.maximum(np.asarray(z_m, dtype=np.float64), 0.0) / scale_height_m) + + +def plant_height(g_planet, gravity_g): + """How much taller plants grow than under the planet's normal gravity (height limit ∝ 1 / g).""" + return g_planet / np.asarray(gravity_g, dtype=np.float64) + + +def run(ctx) -> dict: + pl = ctx.cfg["planet"] + zs, m_o2, m_g = ctx.need("z_surface_m", "m_o2_zones", "m_gravity_zones") + p = pressure(pl["sea_level_pressure_bar"], zs, pl["scale_height_km"] * 1000.0) + o2 = pl["o2_fraction"] * o2_mod(ctx.cfg, m_o2) + return {"po2_bar": o2 * p, "gravity_g": pl["gravity_g"] * gravity_mod(ctx.cfg, m_g), "pressure_bar": p, + "o2_fraction": o2, "fire_reactivity": np.ones(len(p)), # zone events (eras) set other values + "plant_height_x": plant_height(pl["gravity_g"], pl["gravity_g"] * gravity_mod(ctx.cfg, m_g))} diff --git a/mapgen/geo.py b/mapgen/geo.py new file mode 100644 index 0000000..152ae66 --- /dev/null +++ b/mapgen/geo.py @@ -0,0 +1,200 @@ +"""GeoJSON exports: rivers, coastlines/lakes (marching squares), plate boundaries.""" +from __future__ import annotations + +import json +from pathlib import Path + +import h3.api.basic_int as h3 +import numpy as np + +BOUNDARY_NAMES = {1: "convergent", 2: "divergent", 3: "transform"} +# marching-squares cases; corner bits TL=8 TR=4 BR=2 BL=1; edges T R B L +_CASES = {1: [("L", "B")], 2: [("B", "R")], 3: [("L", "R")], 4: [("T", "R")], 5: [("T", "R"), ("L", "B")], + 6: [("T", "B")], 7: [("L", "T")], 8: [("L", "T")], 9: [("T", "B")], 10: [("L", "T"), ("B", "R")], + 11: [("T", "R")], 12: [("L", "R")], 13: [("B", "R")], 14: [("L", "B")]} +_EDGE = {"T": (1, 0), "R": (2, 1), "B": (1, 2), "L": (0, 1)} # doubled (dx, dy) inside a 2x2 block + + +def split_antimeridian(coords): + parts, cur = [], [list(coords[0])] + for a, b in zip(coords, coords[1:]): + if abs(b[0] - a[0]) > 180.0: + parts.append(cur) + cur = [] + cur.append(list(b)) + parts.append(cur) + return [p for p in parts if len(p) >= 2] + + +def _feature(geom_type, coords, props): + return {"type": "Feature", "geometry": {"type": geom_type, "coordinates": coords}, "properties": props} + + +def river_lines(g, recv, river, strahler, discharge): + """One LineString per run of equal Strahler order (runs share their junction vertex).""" + n = g.n + ar = np.arange(n) + has_donor = np.zeros(n, bool) + nr = river & (recv != ar) + has_donor[recv[nr]] = True + visited = np.zeros(n, bool) + feats = [] + for s in np.flatnonzero(river & ~has_donor): + cells, c = [int(s)], int(s) + while True: + visited[c] = True + r = int(recv[c]) + if r == c: + break + cells.append(r) + if not river[r] or visited[r]: + break + c = r + body = cells[:-1] if len(cells) > 1 else cells + start = 0 + for i in range(1, len(body) + 1): + if i == len(body) or strahler[body[i]] != strahler[body[i - 1]]: + run = body[start:i] + [cells[i] if i < len(cells) else body[-1]] + coords = [[round(float(g.lon[k]), 4), round(float(g.lat[k]), 4)] for k in run] + props = {"order": int(strahler[body[start]]), "discharge_km3_yr": round(float(discharge[body[i - 1]]), 3)} + feats += [_feature("LineString", p, props) for p in split_antimeridian(coords)] + start = i + return feats + + +def boundary_features(g, plate, bnd_type): + s, d = g.src, g.dst + e = np.flatnonzero((plate[s] != plate[d]) & (s < d)) + lines = {} + for k in e: + a, b = int(g.ids[s[k]]), int(g.ids[d[k]]) + seg = [[round(lng, 4), round(lat, 4)] for lat, lng in h3.directed_edge_to_boundary(h3.cells_to_directed_edge(a, b))] + if abs(seg[0][0] - seg[-1][0]) > 180: + continue + lines.setdefault(BOUNDARY_NAMES.get(int(bnd_type[s[k]]), "transform"), []).append(seg) + return [_feature("MultiLineString", v, {"type": t}) for t, v in sorted(lines.items())] + + +def contours(mask): + m = np.pad(np.asarray(mask, dtype=np.uint8), 1) + code = m[:-1, :-1] * 8 + m[:-1, 1:] * 4 + m[1:, 1:] * 2 + m[1:, :-1] + adj = {} + for case, segs in _CASES.items(): + ys, xs = np.nonzero(code == case) + for (e1, e2) in segs: + for y, x in zip(ys.tolist(), xs.tolist()): + p = (2 * x + _EDGE[e1][0], 2 * y + _EDGE[e1][1]) + q = (2 * x + _EDGE[e2][0], 2 * y + _EDGE[e2][1]) + adj.setdefault(p, []).append(q) + adj.setdefault(q, []).append(p) + seen, rings = set(), [] + for start in adj: + if start in seen: + continue + ring, prev, cur = [start], None, start + seen.add(start) + while True: + nxt = [q for q in adj[cur] if q != prev] + if not nxt: + break + prev, cur = cur, nxt[0] + ring.append(cur) + if cur == start: + break + seen.add(cur) + rings.append([(px / 2.0 - 1.0, py / 2.0 - 1.0) for px, py in ring]) + return rings + + +def _signed_area(ring): + a = np.asarray(ring, dtype=np.float64) + return 0.5 * float(np.sum(a[:-1, 0] * a[1:, 1] - a[1:, 0] * a[:-1, 1])) + + +def _point_in_ring(pt, ring): + a = np.asarray(ring, dtype=np.float64) + x1, y1, x2, y2 = a[:-1, 0], a[:-1, 1], a[1:, 0], a[1:, 1] + cross = (y1 > pt[1]) != (y2 > pt[1]) + xi = x1 + (pt[1] - y1) * (x2 - x1) / np.where(y2 != y1, y2 - y1, 1e-300) + return bool(np.count_nonzero(cross & (pt[0] < xi)) % 2) + + +def _orient(ring, ccw): + return ring if (_signed_area(ring) > 0) == ccw else ring[::-1] + + +def _clip(ring, left, x0=180.0): + """Sutherland–Hodgman clip of a closed ring to x ≤ x0 (left) or x ≥ x0.""" + inside = (lambda p: p[0] <= x0) if left else (lambda p: p[0] >= x0) + cut = lambda p, q: [x0, p[1] + (x0 - p[0]) / (q[0] - p[0]) * (q[1] - p[1])] + pts, out = ring[:-1], [] + for i in range(len(pts)): + cur, prev = pts[i], pts[i - 1] + if inside(cur): + if not inside(prev): + out.append(cut(prev, cur)) + out.append(cur) + elif inside(prev): + out.append(cut(prev, cur)) + return out + [out[0]] if len(out) >= 3 else [] + + +def _split_seam(rings): + """rings[0] exterior + holes in continuous longitude; split at +180 into ≤2 polygons in [−180, 180].""" + if max(p[0] for p in rings[0]) <= 180.0: + return [rings] + polys = [] + for left in (True, False): + ext = _clip(rings[0], left) + if len(ext) < 4: + continue + holes = [h for h in (_clip(r, left) for r in rings[1:]) if len(h) >= 4] + part = [ext] + holes + if not left: + part = [[[x - 360.0, y] for x, y in r] for r in part] + polys.append(part) + return polys + + +def _round(poly): + return [[[round(x, 4), round(y, 4)] for x, y in r] for r in poly] + + +def contour_features(mask, kind): + """Land/lake outlines as RFC 7946 (Multi)Polygons: CCW exteriors, CW holes, split at the antimeridian.""" + H, W = mask.shape + col = int(np.argmin(mask.sum(axis=0))) + rings = [] + for ring in contours(np.roll(mask, -col, axis=1)): + pts = [[(x + col + 0.5) / W * 360.0 - 180.0, max(-90.0, min(90.0, 90.0 - (y + 0.5) / H * 180.0))] + for x, y in ring] + if len(pts) >= 4 and pts[0] == pts[-1]: + rings.append(pts) + depth = [sum(_point_in_ring(r[0], o) for j, o in enumerate(rings) if j != i) for i, r in enumerate(rings)] + feats = [] + for i, r in enumerate(rings): + if depth[i] % 2: + continue + holes = [_orient(rings[j], False) for j in range(len(rings)) + if depth[j] == depth[i] + 1 and _point_in_ring(rings[j][0], r)] + polys = [_round(p) for p in _split_seam([_orient(r, True)] + holes)] + if len(polys) == 1: + feats.append(_feature("Polygon", polys[0], {"kind": kind})) + elif polys: + feats.append(_feature("MultiPolygon", polys, {"kind": kind})) + return feats + + +def _write(path: Path, feats) -> None: + path.write_text(json.dumps({"type": "FeatureCollection", "features": feats}, separators=(",", ":"))) + + +def write_all(ctx, out_dir: Path, land_raster, lake_raster) -> None: + g, d = ctx.grid, ctx.data + gdir = out_dir / "geo" + gdir.mkdir(parents=True, exist_ok=True) + _write(gdir / "rivers.geojson", river_lines(g, d["recv"], d["river"], d["strahler"], d["discharge_km3_yr"])) + _write(gdir / "plate_boundaries.geojson", boundary_features(g, d["plate"], d["bnd_type"])) + step = max(1, land_raster.shape[1] // 2048) + _write(gdir / "coast.geojson", contour_features(land_raster[::step, ::step], "land")) + _write(gdir / "lakes.geojson", contour_features(lake_raster[::step, ::step] & land_raster[::step, ::step], "lake")) diff --git a/mapgen/graph.py b/mapgen/graph.py new file mode 100644 index 0000000..126c71f --- /dev/null +++ b/mapgen/graph.py @@ -0,0 +1,477 @@ +"""Operations on the H3 cell graph (CSR neighbours).""" +from __future__ import annotations + +import heapq + +import numpy as np +from scipy import sparse +from scipy.sparse import csgraph +from scipy.sparse import linalg as splinalg + + +def nbr_mean(g, f): + return np.bincount(g.src, weights=f[g.dst], minlength=g.n) / g.counts + + +def nbr_max(g, f): + out = np.asarray(f, dtype=np.float64).copy() + np.maximum.at(out, g.src, f[g.dst]) + return out + + +def nbr_min(g, f): + out = np.asarray(f, dtype=np.float64).copy() + np.minimum.at(out, g.src, f[g.dst]) + return out + + +def diffuse(g, f, iters: int, alpha: float = 0.5, mask=None): + f = np.asarray(f, dtype=np.float64) + for _ in range(int(iters)): + new = (1.0 - alpha) * f + alpha * nbr_mean(g, f) + f = new if mask is None else np.where(mask, new, f) + return f + + +def laplacian_matrix(g): + """Graph Laplacian (per km²): Δf_i ≈ (4/k_i) Σ_j (f_j − f_i)/d_ij² (exact for a regular hex lattice).""" + w = 4.0 / (g.counts[g.src] * g.edge_km**2) + lap = sparse.csr_matrix((w, (g.src, g.dst)), shape=(g.n, g.n)) + return lap - sparse.diags(np.asarray(lap.sum(axis=1)).ravel()) + + +def workers() -> int: + """Threads for independent jobs (seasons, components): WORLDGEN_THREADS, default 3. Each job computes exactly + what it would alone, so results never depend on it; numpy and scipy's sparse kernels release the GIL.""" + import os + try: + return max(1, int(os.environ.get("WORLDGEN_THREADS", "3"))) + except ValueError: + return 3 + + +def pmap(fn, items) -> list: + """[fn(x) for x in items] on up to workers() threads, in order.""" + items = list(items) + if workers() <= 1 or len(items) <= 1: + return [fn(x) for x in items] + from concurrent.futures import ThreadPoolExecutor + with ThreadPoolExecutor(min(workers(), len(items))) as ex: + return list(ex.map(fn, items)) + + +def smooth_km(g, f, length_km: float, rtol: float = 1e-6, keep: bool = False): + """Resolution-independent smoothing: solve (I − L²Δ) s = f (screened Poisson, decay length ≈ L km). + keep: keep the matrix on g for the next call (loops over one grid; drop_smooth_cache(g) frees it).""" + f = np.asarray(f, dtype=np.float64) + scale = float(np.max(np.abs(f))) if f.size else 0.0 + if length_km <= 0 or scale == 0.0: + return f.copy() + f = f / scale # linear system: solve at unit scale (avoids breakdown) + cache = g.__dict__.get("_screened", {}) + if length_km in cache: + A, inv_diag = cache[length_km] + else: + A = (sparse.identity(g.n, format="csr") - length_km**2 * laplacian_matrix(g)).tocsr() + inv_diag = 1.0 / A.diagonal() + if keep: + g.__dict__.setdefault("_screened", {})[length_km] = (A, inv_diag) + s, info = bicgstab_jacobi(A, f, f.copy(), inv_diag, rtol, 5000) + if info != 0: + raise ValueError(f"smooth_km: solver did not converge (info={info})") + return s * scale + + +def _make_bicg_jit(): + """bicgstab's vector updates fused into single passes (numba optional; WORLDGEN_NO_JIT=1 turns it off). Each + element gets the same operations in the same order as scipy's numpy statements; dot products, norms and sparse + products stay the very calls scipy makes, so the iterates — and the answer — are the same floats. (Threads were + tried and dropped: on a busy machine they wait more than they work.)""" + import os + if os.environ.get("WORLDGEN_NO_JIT"): + return None + try: + import numba + except ImportError: + return None + + @numba.njit(cache=True, nogil=True) + def p_update(p, v, r, omega, beta): # p -= omega*v; p *= beta; p += r + for i in range(len(p)): + p[i] = (p[i] - omega * v[i]) * beta + r[i] + + @numba.njit(cache=True, nogil=True) + def scale(out, d, x): # out = d * x (the Jacobi preconditioner) + for i in range(len(x)): + out[i] = d[i] * x[i] + + @numba.njit(cache=True, nogil=True) + def axpy_neg(r, a, v): # r -= a*v + for i in range(len(r)): + r[i] = r[i] - a * v[i] + + @numba.njit(cache=True, nogil=True) + def x_update(x, alpha, phat, omega, shat): # x += alpha*phat; x += omega*shat + for i in range(len(x)): + x[i] = (x[i] + alpha * phat[i]) + omega * shat[i] + + @numba.njit(cache=True, nogil=True) + def axpy(x, a, v): # x += a*v + for i in range(len(x)): + x[i] = x[i] + a * v[i] + + return p_update, scale, axpy_neg, x_update, axpy + + +_bicg_jit = _make_bicg_jit() + + +def _csr_matvec_into(A): + """mv(x, y): y = A @ x, by scipy's own kernel into a reused buffer (A @ x zero-fills a fresh array and calls the + same csr_matvec; a fresh 16 MB array costs more in page faults than the product).""" + try: + from scipy.sparse import _sparsetools + fn = _sparsetools.csr_matvec + except (ImportError, AttributeError): + fn = None + M, N = A.shape + + def mv(x, y): + if fn is None: + y[:] = A @ x + return + y.fill(0.0) + fn(M, N, A.indptr, A.indices, A.data, x, y) + return mv + + +JIT_MIN_N = 50_000 # smaller systems: scipy (thread start-up outweighs the gain; same answer either way) + + +def bicgstab_jacobi(A, b, x0, inv_diag, rtol, maxiter): + """scipy.sparse.linalg.bicgstab(A, b, x0=x0, rtol=rtol, maxiter=maxiter, M=diag(inv_diag)) for a CSR matrix and + float64 vectors: the same iterates, statement by statement (scipy 1.12+'s pure-Python loop), with the vector + updates fused and everything written into preallocated buffers (memory-bound; fresh 16 MB temporaries cost + more in page faults than in arithmetic). Falls back to scipy without numba.""" + b = np.asarray(b, dtype=np.float64).ravel() + if (_bicg_jit is None or A.shape[0] < JIT_MIN_N or not sparse.isspmatrix_csr(A) + and not isinstance(A, sparse.csr_array) or A.dtype != np.float64): + M = splinalg.LinearOperator(A.shape, matvec=lambda x: inv_diag * x) + return splinalg.bicgstab(A, b, x0=x0, rtol=rtol, maxiter=maxiter, M=M) + p_update, scale, axpy_neg, x_update, axpy = _bicg_jit + mv = _csr_matvec_into(A) + inv_diag = np.ascontiguousarray(inv_diag, dtype=np.float64) + x = np.array(x0, dtype=np.float64) + bnrm2 = np.linalg.norm(b) + atol = max(0.0, float(rtol) * float(bnrm2)) + if bnrm2 == 0: + return b, 0 + rhotol = np.finfo(x.dtype.char).eps ** 2 + omegatol = rhotol + rho_prev, omega, alpha, p, v = None, None, None, None, None + r = b - A @ x if x.any() else b.copy() + rtilde = r.copy() + phat, shat, v, t = np.empty_like(r), np.empty_like(r), np.empty_like(r), np.empty_like(r) + for iteration in range(maxiter): + if np.linalg.norm(r) < atol: + return x, 0 + rho = np.dot(rtilde, r) + if np.abs(rho) < rhotol: + return x, -10 + if iteration > 0: + if np.abs(omega) < omegatol: + return x, -11 + beta = (rho / rho_prev) * (alpha / omega) + p_update(p, v, r, omega, beta) + else: + p = r.copy() + scale(phat, inv_diag, p) + mv(phat, v) + rv = np.dot(rtilde, v) + if rv == 0: + return x, -11 + alpha = rho / rv + axpy_neg(r, alpha, v) + s = r # scipy copies r into s here and reads both unchanged until r -= omega*t + if np.linalg.norm(s) < atol: + axpy(x, alpha, phat) + return x, 0 + scale(shat, inv_diag, s) + mv(shat, t) + omega = np.dot(t, s) / np.dot(t, t) + x_update(x, alpha, phat, omega, shat) + axpy_neg(r, omega, t) + rho_prev = rho + return x, maxiter + + +def drop_smooth_cache(g) -> None: + """Free the matrices smooth_km keeps on g.""" + g.__dict__.pop("_screened", None) + + +def gradient(g, f): + """Tangent-plane gradient (f per km): (2/k) Σ_j (f_j − f_i)/d_ij · t_ij.""" + w = (f[g.dst] - f[g.src]) / g.edge_km + s = np.stack([np.bincount(g.src, weights=w * g.edge_tangents[:, c], minlength=g.n) for c in range(3)], axis=1) + return s * (2.0 / g.counts)[:, None] + + +def _csr(g, weights): + return sparse.csr_matrix((weights, (g.src, g.dst)), shape=(g.n, g.n)) + + +def nearest_source(g, sources, weights=None): + sources = np.asarray(sources, dtype=np.int64) + if len(sources) == 0: + return np.full(g.n, np.inf), np.full(g.n, -9999, dtype=np.int64) + w = g.edge_km if weights is None else weights + dist, _, src = csgraph.dijkstra(_csr(g, w), directed=True, indices=sources, + min_only=True, return_predecessors=True) + return dist, src.astype(np.int64) + + +def distance_to(g, mask): + return nearest_source(g, np.flatnonzero(mask))[0] + + +def priority_flood(g, z, sink_mask, eps: float = 0.01): + """Barnes (2014) priority-flood + ε: every non-sink cell gets a strictly descending path to a sink.""" + sink_mask = np.asarray(sink_mask, dtype=bool) + if not sink_mask.any(): + raise ValueError("priority_flood: no sink cells") + has_open = np.bincount(g.src, weights=(~sink_mask)[g.dst].astype(np.float64), minlength=g.n) > 0 + seeds = np.flatnonzero(sink_mask & has_open) + z = np.asarray(z, dtype=np.float64) + if _flood_jit is not None: + return _flood_jit(z.copy(), sink_mask.copy(), np.asarray(g.nbr_ptr, np.int64), np.asarray(g.nbr_idx, np.int64), + seeds.astype(np.int64), float(eps)) + return _flood_py(z, sink_mask, g.nbr_ptr, g.nbr_idx, seeds, eps) + + +def _flood_py(z, sink_mask, nbr_ptr, nbr_idx, seeds, eps): + zf = np.asarray(z, dtype=np.float64).tolist() + done = np.asarray(sink_mask).tolist() + ptr = np.asarray(nbr_ptr).tolist() + idx = np.asarray(nbr_idx).tolist() + heap = [(zf[i], i) for i in np.asarray(seeds).tolist()] + heapq.heapify(heap) + while heap: + zc, c = heapq.heappop(heap) + for k in range(ptr[c], ptr[c + 1]): + n = idx[k] + if not done[n]: + done[n] = True + if zf[n] < zc + eps: + zf[n] = zc + eps + heapq.heappush(heap, (zf[n], n)) + return np.array(zf) + + +def _make_flood_jit(): + """_flood_py compiled with numba when it is installed (optional: same heap order, same float steps, same result; + WORLDGEN_NO_JIT=1 turns it off).""" + import os + if os.environ.get("WORLDGEN_NO_JIT"): + return None + try: + import numba + except ImportError: + return None + + @numba.njit(cache=True) + def flood(zf, done, ptr, idx, seeds, eps): + heap = [(zf[i], i) for i in seeds] + heapq.heapify(heap) + while len(heap): + zc, c = heapq.heappop(heap) + for k in range(ptr[c], ptr[c + 1]): + n = idx[k] + if not done[n]: + done[n] = True + if zf[n] < zc + eps: + zf[n] = zc + eps + heapq.heappush(heap, (zf[n], n)) + return zf + return flood + + +_flood_jit = _make_flood_jit() + + +def _make_steep_jit(): + """steepest_receivers' per-row maximum, compiled (numba optional, as the flood): the first edge in row order + with the largest slope, slopes computed as numpy does. ok=False (empty row, NaN): use the numpy path.""" + import os + if os.environ.get("WORLDGEN_NO_JIT"): + return None + try: + import numba + except ImportError: + return None + + @numba.njit(cache=True) + def steep(z, ptr, dst, edge): + n = len(ptr) - 1 + first = np.empty(n, np.int64) + best = np.empty(n) + for i in range(n): + a, b = ptr[i], ptr[i + 1] + if a == b: + return first, best, False + bk = a + bs = (z[i] - z[dst[a]]) / edge[a] + if bs != bs: + return first, best, False + for k in range(a + 1, b): + sk = (z[i] - z[dst[k]]) / edge[k] + if sk != sk: + return first, best, False + if sk > bs: + bs, bk = sk, k + first[i], best[i] = bk, bs + return first, best, True + return steep + + +_steep_jit = _make_steep_jit() + + +def steepest_receivers(g, z): + """Per cell: the neighbour of steepest descent (first in neighbour order on ties), the slope and the edge length; + no way down → itself, 0, inf.""" + ptr = np.asarray(g.nbr_ptr, np.int64) + if _steep_jit is not None and g.n: + first, s, ok = _steep_jit(np.asarray(z), ptr, np.asarray(g.nbr_idx, np.int64), g.edge_km) + if ok: + down = s > 0 + recv = np.where(down, g.dst[first], np.arange(g.n)) + return recv, np.where(down, s, 0.0), np.where(down, g.edge_km[first], np.inf) + slope = (z[g.src] - z[g.dst]) / g.edge_km + if g.n == 0 or not (np.diff(ptr) > 0).all() or np.isnan(slope).any(): + order = np.lexsort((-slope, g.src)) # general case (empty rows, NaN) + first = order[ptr[:-1]] + else: # same edge as the stable lexsort, without sorting + smax = np.maximum.reduceat(slope, ptr[:-1]) + cand = np.flatnonzero(slope == smax[g.src]) + rows = g.src[cand] + first = cand[np.concatenate([[True], rows[1:] != rows[:-1]])] + s = slope[first] + down = s > 0 + ar = np.arange(g.n) + recv = np.where(down, g.dst[first], ar) + return recv, np.where(down, s, 0.0), np.where(down, g.edge_km[first], np.inf) + + +def _gather(ptr, arr, sel): + counts = ptr[sel + 1] - ptr[sel] + tot = int(counts.sum()) + if tot == 0: + return arr[:0] + starts = np.repeat(ptr[sel] - np.concatenate([[0], np.cumsum(counts)[:-1]]), counts) + return arr[starts + np.arange(tot)] + + +def receiver_levels(recv): + recv = np.asarray(recv, dtype=np.int64) + n = len(recv) + ar = np.arange(n) + root = recv == ar + donors = ar[~root] + donors = donors[np.argsort(recv[donors], kind="stable")] + dptr = np.concatenate([[0], np.cumsum(np.bincount(recv[donors], minlength=n))]) + levels, frontier, seen = [], ar[root], 0 + while len(frontier): + levels.append(frontier) + seen += len(frontier) + frontier = _gather(dptr, donors, frontier) + if seen != n: + raise ValueError("receiver_levels: cycle in receivers") + return levels + + +def accumulate(recv, levels, w): + """Sum w down the receiver tree. Per level only the receivers are touched (a full bincount per level costs + levels × cells); the sums are added in donor order from 0, as bincount does: the same floats.""" + acc = np.asarray(w, dtype=np.float64).copy() + buf = np.zeros(len(acc)) + if len(levels) > 1: + acc += 0.0 # as the first full-length add did: −0 becomes +0 + for lv in reversed(levels[1:]): + r = recv[lv] + buf[r] = 0.0 + np.add.at(buf, r, acc[lv]) + acc[r] = acc[r] + buf[r] + return acc + + +def _make_components_jit(): + """Connected-component labels as scipy's connected_components numbers them — each component (every node not in + the mask is one by itself) by the order of its lowest node — by union-find over the edges, with no sparse matrix + (numba optional; WORLDGEN_NO_JIT=1 turns it off).""" + import os + if os.environ.get("WORLDGEN_NO_JIT"): + return None + try: + import numba + except ImportError: + return None + + @numba.njit(cache=True, nogil=True) + def labels(n, src, dst, mask): + parent = np.arange(n) + for e in range(src.shape[0]): + a, b = src[e], dst[e] + if mask[a] and mask[b]: + while parent[a] != a: + parent[a] = parent[parent[a]] + a = parent[a] + while parent[b] != b: + parent[b] = parent[parent[b]] + b = parent[b] + if a != b: + if a < b: + parent[b] = a + else: + parent[a] = b + root_lab = np.full(n, -1, np.int32) + out = np.empty(n, np.int32) + count = 0 + for v in range(n): + r = v + while parent[r] != r: + r = parent[r] + if root_lab[r] < 0: + root_lab[r] = count + count += 1 + out[v] = root_lab[r] if mask[v] else -1 + return out + return labels + + +_components_jit = _make_components_jit() + + +def components(g, mask): + mask = np.asarray(mask, dtype=bool) + if _components_jit is not None: + return _components_jit(g.n, g.src, g.dst, mask) + e = mask[g.src] & mask[g.dst] + m = sparse.csr_matrix((np.ones(int(e.sum())), (g.src[e], g.dst[e])), shape=(g.n, g.n)) + _, lab = csgraph.connected_components(m, directed=False) + return np.where(mask, lab, -1) + + +OCEAN_MIN_KM2 = 5.0e6 + + +def ocean_mask(g, z, min_area_km2: float = OCEAN_MIN_KM2): + """The connected world ocean plus any separate basin ≥ min_area_km2; smaller interior lows count as land.""" + wet = np.asarray(z) <= 0 + if not wet.any(): + return wet + lab = components(g, wet) + area = np.bincount(lab[wet], weights=g.area_km2[wet]) + keep = area >= min_area_km2 + keep[np.argmax(area)] = True + return wet & keep[np.maximum(lab, 0)] diff --git a/mapgen/grid.py b/mapgen/grid.py new file mode 100644 index 0000000..8138b5f --- /dev/null +++ b/mapgen/grid.py @@ -0,0 +1,121 @@ +"""H3 hexagonal grid as a CSR cell graph. Stage `grid`.""" +from __future__ import annotations + +from dataclasses import dataclass +from functools import cached_property + +import h3.api.basic_int as h3 +import numpy as np + +from .sphere import latlon_to_xyz, tangent_dir + + +EDGE_BLOCK = 1 << 20 + + +def edge_blocks(n: int, block: int | None = None): + """Slices over n edges, EDGE_BLOCK at a time: per-edge (row-wise) formulas give the same values block by block, + with (block, 3) float64 temporaries instead of (edges, 3) ones (290 MB each at r5, 2 GB at r6).""" + block = block or EDGE_BLOCK + return [slice(a, min(a + block, n)) for a in range(0, n, block)] + + +def rowdot_at(A, idx, t) -> np.ndarray: + """np.sum(A[idx] * t, axis=1) — per edge, A's row at an end dotted with the edge's vector — in edge blocks.""" + out = np.empty(len(idx)) + for s in edge_blocks(len(idx)): + out[s] = np.sum(A[idx[s]] * t[s], axis=1) + return out + + +def cell_parents(cells, res: int) -> np.ndarray: + """h3.cell_to_parent for an array of cells (all of resolution ≥ res), by the index bits: the resolution field + set to res, the digits below it set to 7 (unused). Same ids as h3, without a Python call per cell.""" + cells = np.asarray(cells, dtype=np.uint64) + if len(cells) and int(((cells >> np.uint64(52)) & np.uint64(15)).min()) < res: + raise ValueError(f"cell_parents: a cell is coarser than resolution {res}") + unused = 0 + for r in range(res + 1, 16): + unused |= 7 << ((15 - r) * 3) + out = (cells & np.uint64(~(15 << 52) & (2**64 - 1))) | np.uint64(res << 52) + return out | np.uint64(unused) + + +@dataclass +class Grid: + res: int + radius_km: float + ids: np.ndarray + lat: np.ndarray + lon: np.ndarray + xyz: np.ndarray + area_km2: np.ndarray + nbr_ptr: np.ndarray + nbr_idx: np.ndarray + + @property + def n(self) -> int: + return len(self.ids) + + @cached_property + def counts(self): + return np.diff(self.nbr_ptr) + + @cached_property + def src(self): + return np.repeat(np.arange(self.n), self.counts) + + @property + def dst(self): + return self.nbr_idx + + @cached_property + def edge_km(self): + out = np.empty(len(self.dst)) + for s in edge_blocks(len(self.dst)): + d = np.sum(self.xyz[self.src[s]] * self.xyz[self.dst[s]], axis=1) + out[s] = self.radius_km * np.arccos(np.clip(d, -1.0, 1.0)) + return out + + @cached_property + def edge_tangents(self): + out = np.empty((len(self.dst), 3)) + for s in edge_blocks(len(self.dst)): + out[s] = tangent_dir(self.xyz[self.src[s]], self.xyz[self.dst[s]]) + return out + + @cached_property + def spacing_km(self) -> float: + return float(self.edge_km.mean()) + + def cell_index(self, lat: float, lon: float) -> int: + c = np.uint64(h3.latlng_to_cell(float(lat), float(lon), self.res)) + return int(np.searchsorted(self.ids, c)) + + def to_arrays(self) -> dict: + return {"g_ids": self.ids, "g_lat": self.lat, "g_lon": self.lon, "g_xyz": self.xyz, + "g_area_km2": self.area_km2, "g_nbr_ptr": self.nbr_ptr, "g_nbr_idx": self.nbr_idx} + + @classmethod + def from_arrays(cls, d: dict, res: int, radius_km: float) -> "Grid": + return cls(res, radius_km, d["g_ids"], d["g_lat"], d["g_lon"], d["g_xyz"], d["g_area_km2"], + d["g_nbr_ptr"], d["g_nbr_idx"]) + + +def build_grid(res: int, radius_km: float) -> Grid: + cells: list[int] = [] + for r0 in h3.get_res0_cells(): + cells.extend(h3.cell_to_children(r0, res)) + ids = np.array(sorted(cells), dtype=np.uint64) + ll = np.array([h3.cell_to_latlng(int(c)) for c in ids], dtype=np.float64) + area = np.array([h3.cell_area(int(c), unit="rads^2") for c in ids]) * radius_km**2 + rings = [h3.grid_ring(int(c), 1) for c in ids] + counts = np.array([len(r) for r in rings], dtype=np.int64) + ptr = np.concatenate([[0], np.cumsum(counts)]).astype(np.int64) + flat = np.array([c for r in rings for c in r], dtype=np.uint64) + idx = np.searchsorted(ids, flat).astype(np.int64) + return Grid(res, radius_km, ids, ll[:, 0], ll[:, 1], latlon_to_xyz(ll[:, 0], ll[:, 1]), area, ptr, idx) + + +def run(ctx) -> dict: + return build_grid(ctx.res, float(ctx.cfg["planet"]["radius_km"])).to_arrays() diff --git a/mapgen/hydrology.py b/mapgen/hydrology.py new file mode 100644 index 0000000..d7adde1 --- /dev/null +++ b/mapgen/hydrology.py @@ -0,0 +1,231 @@ +"""Stage `hydrology`: depression filling, lakes vs endorheic basins, discharge, rivers, Strahler order.""" +from __future__ import annotations + +import numpy as np +from scipy import sparse +from scipy.sparse import csgraph + +from .config import params +from .crust import RIFT +from .graph import accumulate, components, distance_to, priority_flood, receiver_levels, steepest_receivers +from .noise import fbm +from .pipeline import StageError + +DEFAULTS = {"fill_eps_m": 0.01, "min_depth_m": 1.0, "river_min_km3_yr": 2.0, + "salt_flat_max_p_mm": 300.0, "salt_flat_fraction": 0.2, "dry_lake_fraction": 0.05, + # lake beds: deepest point = k × area^exp m (Earth-like: ≈95 m at + # 1,000 km², ≈230 m at 20,000 km², ≈610 m at 500,000 km²) × 0.5–2 (noise at the lake), × rift_x in + # rifts, × arid_x for dry terminal lakes; never shallower than lake_min_m (refinement keeps ≥ 15 m) + "lake_depth_k": 12.0, "lake_depth_exp": 0.3, "lake_min_m": 25.0, "lake_max_m": 1800.0, + "lake_rift_x": 2.5, "lake_arid_x": 0.4, "lake_arid_p_mm": 400.0, "lake_shore": 0.25} + + +def strahler(recv, levels, river): + n = len(recv) + order = np.zeros(n, np.int8) + mx = np.zeros(n, np.int8) + cnt = np.zeros(n, np.int16) + for lv in reversed(levels): + cells = lv[river[lv]] + if len(cells) == 0: + continue + order[cells] = np.where(cnt[cells] >= 2, mx[cells] + 1, np.maximum(mx[cells], 1)) + nonroot = cells[recv[cells] != cells] + r, o = recv[nonroot], order[nonroot] + np.maximum.at(mx, r, o) + np.add.at(cnt, r, (o == mx[r]).astype(np.int16)) + return order + + +def _groups(lab): + idx = np.argsort(lab, kind="stable") + ls = lab[idx] + start = np.searchsorted(ls, 0) + idx, ls = idx[start:], ls[start:] + cuts = np.flatnonzero(np.diff(ls)) + 1 + return np.split(idx, cuts) if len(idx) else [] + + +def _leaves(recv, lab, x, limit=100000): + """Does the flow from depression cell x leave its depression for good (not back in over shallow ground)?""" + own, y = lab[x], recv[x] + for _ in range(limit): + if lab[y] == own: + return False + if lab[y] >= 0 or recv[y] == y: + return True + y = recv[y] + return True + + +def _leaves_all(recv, lab, levels, limit=100000): + """_leaves for every cell at once: walking down from recv[x], the first cell that is in a depression or a root + (stop) and how many steps away it is; x leaves unless that cell is in x's own depression (or the walk is longer + than limit steps).""" + n = len(recv) + stop, d = np.arange(n), np.zeros(n, np.int64) + for lv in levels[1:]: # roots first: a cell's receiver is done before it + free = lv[lab[lv] < 0] + stop[free] = stop[recv[free]] + d[free] = d[recv[free]] + 1 + y = recv + return (d[y] >= limit) | (lab[stop[y]] != lab) + + +def _route_inside(g, dep, lab, recv, targets): + """Within each depression, point every cell along the shortest intra-depression path to its target cell.""" + e = dep[g.src] & dep[g.dst] & (lab[g.src] == lab[g.dst]) + m = sparse.csr_matrix((g.edge_km[e], (g.src[e], g.dst[e])), shape=(g.n, g.n)) + _, pred, _ = csgraph.dijkstra(m, directed=False, indices=targets, min_only=True, return_predecessors=True) + out = recv.copy() + inside = dep & (pred >= 0) + out[inside] = pred[inside] + return out + + +def _sweep(recv, levels, water, outlets, cap): + """Accumulate flow upstream→downstream; at each spilling-lake outlet remove up to `cap` (lake evaporation).""" + acc = np.asarray(water, dtype=np.float64).copy() + loss = np.zeros(len(acc)) + is_out = np.zeros(len(acc), bool) + is_out[outlets] = True + cap_cell = np.zeros(len(acc)) + cap_cell[outlets] = cap + buf = np.zeros(len(acc)) # per level only the receivers change (graph.accumulate): same floats + if len(levels) > 1: + acc += 0.0 + for lv in reversed(levels[1:]): + push = acc[lv].copy() + o = is_out[lv] + if o.any(): + cells = lv[o] + lost = np.minimum(acc[cells], cap_cell[cells]) + loss[cells] = lost + push[o] = acc[cells] - lost + r = recv[lv] + buf[r] = 0.0 + np.add.at(buf, r, push) + acc[r] = acc[r] + buf[r] + return acc, loss + + +def lake_levels(z, zf, lab, lake, lake_id, endo): + """Each lake cell's water level (NaN elsewhere): a spilling lake stands at its spill height; a terminal lake at the + lowest ground of its basin it does not cover (the next cell to flood). Carving the beds leaves both unchanged.""" + level = np.full(len(z), np.nan) + m = int(lab.max()) + 1 if len(lab) and lab.max() >= 0 else 0 + dry = np.full(m, np.inf) # per depression: its lowest uncovered ground + sel = (lab >= 0) & ~lake + np.minimum.at(dry, lab[sel], z[sel]) + for c in _groups(np.where(lake, lake_id, -1)): + k = lab[c[0]] + level[c] = dry[k] if endo[c[0]] and k >= 0 and np.isfinite(dry[k]) else zf[c].min() + return level + + +def lake_beds(g, z, level, lake, lake_id, endo, rift, p_ann, seed, P, only=None): + """Lake beds carved below their level (heights unchanged elsewhere): each lake's deepest point from its area + (P lake_*), noise keyed by where the lake lies (not its id), the depth rising from lake_shore × that at the shore + to all of it at the cell farthest from shore; never above the ground (min), never shallower than lake_min_m. + only: carve just the lakes touching these cells (eras: where events changed the ground).""" + z = np.asarray(z, dtype=np.float64).copy() + if not lake.any(): + return z + shore = distance_to(g, ~lake) if (~lake).any() else np.full(g.n, g.spacing_km) + for c in _groups(np.where(lake, lake_id, -1)): + if only is not None and not only[c].any(): + continue + area = g.area_km2[c].sum() + ctr = g.xyz[c].mean(0) + ctr = ctr / max(np.linalg.norm(ctr), 1e-12) + d = P["lake_depth_k"] * area ** P["lake_depth_exp"] * 2.0 ** float(fbm(ctr[None], seed + 4421, 3, 6.0)[0] * 1.4) + if rift[c].mean() > 0.3: + d *= P["lake_rift_x"] + if endo[c[0]] and p_ann[c].mean() < P["lake_arid_p_mm"]: + d *= P["lake_arid_x"] + d = float(np.clip(d, P["lake_min_m"], P["lake_max_m"])) + t = shore[c] / max(shore[c].max(), 1e-9) + prof = np.maximum(d * (P["lake_shore"] + (1 - P["lake_shore"]) * t ** 0.6), P["lake_min_m"]) + z[c] = np.minimum(z[c], level[c] - prof) + return z + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "hydrology", DEFAULTS) + z, p_ann, pet = ctx.need("elevation_eroded_m", "P_ann", "PET") + z = z.astype(np.float64) + ar = np.arange(g.n) + ocean = np.asarray(ctx.data["ocean"]) if "ocean" in ctx.data else z <= 0 + aet = p_ann / np.sqrt(1.0 + (p_ann / np.maximum(pet, 1e-6)) ** 2) # Pike (1964) + runoff = np.where(ocean, 0.0, np.maximum(p_ann - aet, 0.0)) + water = runoff * g.area_km2 * 1e-6 # km³/yr per cell + zf = priority_flood(g, z, ocean, P["fill_eps_m"]) + depth = zf - z + dep = ~ocean & (depth > P["min_depth_m"]) + lab = components(g, dep) + recv1, _, _ = steepest_receivers(g, zf) + recv1 = np.where(ocean, ar, recv1) + levels1 = receiver_levels(recv1) + q1 = accumulate(recv1, levels1, water) + evap_net = np.maximum(pet - p_ann, 0.0) * g.area_km2 * 1e-6 # full-lake evaporation, km³/yr + + groups = _groups(lab) + leaves = _leaves_all(recv1, lab, levels1) if len(groups) else None + outlets = np.zeros(len(groups), np.int64) + terminals = np.zeros(len(groups), np.int64) + cap = np.zeros(len(groups)) + for k, c in enumerate(groups): + ext = c[lab[recv1[c]] != lab[c[0]]] + ext = ext[leaves[ext]] + outlets[k] = ext[np.argmax(q1[ext])] if len(ext) else c[np.argmax(q1[c])] + terminals[k] = c[np.argmin(z[c])] + cap[k] = evap_net[c].sum() + + # 1) every depression drains to its spill outlet; one upstream-first sweep decides spill vs endorheic + recv_a = _route_inside(g, dep, lab, recv1, outlets) if len(groups) else recv1 + acc_a, _ = _sweep(recv_a, receiver_levels(recv_a), water, outlets, cap) + inflow = acc_a[outlets] + spill = inflow >= cap + # 2) endorheic depressions drain to their single lowest cell instead + targets = np.where(spill, outlets, terminals) + recv = _route_inside(g, dep, lab, recv1, targets) if len(groups) else recv1 + recv[terminals[~spill]] = terminals[~spill] + try: + levels = receiver_levels(recv) + except ValueError as e: + raise StageError(f"hydrology: {e}") from e + q, loss = _sweep(recv, levels, water, outlets[spill], cap[spill]) + + lake = np.zeros(g.n, bool) + endo = np.zeros(g.n, bool) + salt = np.zeros(g.n, bool) + lake_id = np.full(g.n, -1, np.int32) + for k, c in enumerate(groups): + if spill[k]: + lake[c] = True + lake_id[c] = k + continue + endo[c] = True + cs = c[np.argsort(z[c])] + nl = int(np.searchsorted(np.cumsum(evap_net[cs]), inflow[k])) + if inflow[k] > 0: + nl = max(nl, 1) # the terminal always holds some water + lake[cs[:nl]] = True + lake_id[cs[:nl]] = k + if nl == 0 or (nl < P["dry_lake_fraction"] * len(cs) and p_ann[c].mean() < P["salt_flat_max_p_mm"]): + a = np.cumsum(g.area_km2[cs]) + ns = max(nl + 1, int(np.searchsorted(a, P["salt_flat_fraction"] * a[-1]))) + salt[cs[nl:ns]] = True + river = ~ocean & ~lake & (q >= P["river_min_km3_yr"]) + level = lake_levels(z, zf, lab, lake, lake_id, endo) + rift = np.asarray(ctx.data["age_class"]) == RIFT if "age_class" in ctx.data else np.zeros(g.n, bool) + zb = lake_beds(g, z, level, lake, lake_id, endo, rift, p_ann, ctx.seed, P, ctx.data.get("lake_carve")) + cut = np.asarray(ctx.data.get("lake_cut_m", 0.0), dtype=np.float64) + (z - zb) # sea masks see the uncut ground + z = zb + depth = zf - z + dtype = np.asarray(ctx.data["elevation_eroded_m"]).dtype + return {"elevation_eroded_m": z.astype(dtype), "lake_level_m": level.astype(np.float32), + "lake_cut_m": np.broadcast_to(cut, z.shape).astype(np.float32), "z_filled_m": zf, "recv": recv, "discharge_km3_yr": q, "runoff_mm": runoff, "aet_mm": aet, + "lake": lake, "lake_id": lake_id, "endorheic": endo, "salt_flat": salt, "river": river, + "strahler": strahler(recv, levels, river), "depression_depth_m": depth, "lake_loss_km3_yr": loss} diff --git a/mapgen/ice.py b/mapgen/ice.py new file mode 100644 index 0000000..2aef33f --- /dev/null +++ b/mapgen/ice.py @@ -0,0 +1,37 @@ +"""Stage `ice`: ice sheets, mountain glaciers, seasonal/perennial sea ice.""" +from __future__ import annotations + +import numpy as np + +from .config import params +from .graph import components, distance_to + +ICE_NONE, ICE_SHEET, ICE_GLACIER, ICE_SEA_SEASONAL, ICE_SEA_PERENNIAL = range(5) +ICE_NAMES = ["none", "ice sheet", "glacier", "seasonal sea ice", "perennial sea ice"] +DEFAULTS = {"melt_summer_c": 0.0, "sheet_min_km2": 250000.0, "sheet_min_p_mm": 100.0, "sea_ice_t_c": -1.8, "sheet_max_m": 3000.0, + "sheet_edge_m": 200.0, "sheet_growth_m_per_km": 3.0} + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "ice", DEFAULTS) + z, t_mean, t_jun, t_dec, p_ann = ctx.need("elevation_eroded_m", "T_mean", "T_jun", "T_dec", "P_ann") + z = z.astype(np.float64) + land = ~np.asarray(ctx.data["ocean"]) if "ocean" in ctx.data else z > 0 + summer = np.where(g.lat >= 0, t_jun, t_dec) + winter = np.where(g.lat >= 0, t_dec, t_jun) + cold = land & (summer < P["melt_summer_c"]) # snow survives the summer → permanent ice + lab = components(g, cold) + area = np.bincount(lab[cold], weights=g.area_km2[cold]) if cold.any() else np.zeros(1) + big = cold & (area[np.maximum(lab, 0)] >= P["sheet_min_km2"]) + sheet = big & (p_ann > P["sheet_min_p_mm"]) + thick = np.zeros(g.n) + if sheet.any(): + d_edge = distance_to(g, ~sheet) + thick = np.where(sheet, np.minimum(P["sheet_max_m"], P["sheet_edge_m"] + P["sheet_growth_m_per_km"] * d_edge), 0.0) + ice = np.zeros(g.n, np.int8) + ice[cold & ~sheet] = ICE_GLACIER + ice[sheet] = ICE_SHEET + ice[~land & (winter < P["sea_ice_t_c"])] = ICE_SEA_SEASONAL + ice[~land & (summer < P["sea_ice_t_c"])] = ICE_SEA_PERENNIAL + return {"ice": ice, "ice_thickness_m": thick, "z_surface_m": z + thick} diff --git a/mapgen/minerals.py b/mapgen/minerals.py new file mode 100644 index 0000000..ab173c9 --- /dev/null +++ b/mapgen/minerals.py @@ -0,0 +1,172 @@ +"""Mineral deposits: where each ore, fuel and industrial mineral can be mined, by +geological rules, placed at an Earth-like share of the land. Everything a society needs to reach steam, rifled guns +and aluminium airships is somewhere; the good stuff (high-grade iron ore and the alloy metals) lies mostly deep in +mountains. Each deposit is a bit of `deposits` (uint32); `deposit_main` names the rarest one in a cell (a map layer). + +A deposit's cells: the land cells (not sea, not under lakes) where its rule allows it, scored by the rule × a district +noise of its own (deposits cluster into mining districts). Half the `share` is taken province by province (≈ 1500 km +cells, each its share of its own land, best-scored first): every region has its local bog iron, coal pits and small +mines where its rocks allow. The rest goes to the best-scored land world-wide (the great districts).""" +from __future__ import annotations + +import numpy as np + +from scipy.spatial import cKDTree + +from .noise import fbm, name_seed + +# name, label, share of land (fraction). Shares are rough Earth-like extents of workable districts, not grades. +DEPOSITS = [ + ("bog_iron", "bog / laterite iron (low grade)", 0.030), + ("ironstone", "ironstone (sedimentary iron)", 0.025), + ("iron_high", "high-grade iron (magnetite, banded)", 0.006), + ("coal", "coal (incl. coking)", 0.050), + ("copper", "copper", 0.012), + ("molybdenum", "molybdenum", 0.002), + ("tin", "tin", 0.003), + ("tungsten", "tungsten", 0.002), + ("lead_zinc", "lead & zinc", 0.010), + ("silver", "silver", 0.005), + ("gold", "gold", 0.008), + ("manganese", "manganese", 0.004), + ("chromium", "chromium", 0.0015), + ("nickel", "nickel", 0.003), + ("cobalt", "cobalt", 0.001), + ("vanadium", "vanadium & titanium", 0.002), + ("platinum", "platinum metals", 0.0005), + ("mercury", "mercury (cinnabar)", 0.0012), + ("sulfur", "sulfur", 0.004), + ("saltpetre", "saltpetre (nitrates)", 0.004), + ("bauxite", "bauxite (aluminium)", 0.010), + ("fluorite", "fluorite / cryolite (flux)", 0.0025), + ("rock_salt", "rock salt & evaporites", 0.018), + ("potash", "potash", 0.004), + ("phosphate", "phosphate", 0.005), + ("oil_gas", "oil & gas", 0.025), + ("helium", "helium (with gas)", 0.0025), + ("diamond", "diamonds (kimberlite)", 0.0004), + ("magnesite", "magnesite / dolomite (magnesium)", 0.003), + ("kaolin", "kaolin / fire clay", 0.010), + ("glass_sand", "glass sand (quartz)", 0.015), + ("graphite", "graphite", 0.002), +] +NAMES = [d[0] for d in DEPOSITS] +BIT = {n: 1 << i for i, n in enumerate(NAMES)} +LEGEND = ["—"] + [d[1] for d in DEPOSITS] +# the good stuff: mostly in deep mountains +DEEP = ("iron_high", "tungsten", "molybdenum", "chromium", "cobalt", "vanadium", "platinum", "manganese", "tin") + + +PROVINCE_KM = 1500.0 # spacing of the provinces that each take their local share +LOCAL = 0.5 # part of a deposit's share taken province by province +LOCAL_MIN = 0.25 # ... only on ground at least this good relative to the world's marginal district + + +def provinces(g, seed): + """Province index of each cell: the nearest of random centres ≈ PROVINCE_KM apart.""" + n = max(1, int(round(4 * np.pi * g.radius_km ** 2 / PROVINCE_KM ** 2))) + c = np.random.default_rng(name_seed(seed, "deposit-provinces")).normal(size=(n, 3)) + c /= np.linalg.norm(c, axis=1)[:, None] + return cKDTree(c).query(g.xyz / np.linalg.norm(g.xyz, axis=1)[:, None])[1] + + +def _rank_select(score, share, land, prov): + """`share` × land cells where score > 0: LOCAL of it the best of each province (its share of its own land, on + ground ≥ LOCAL_MIN × the world's k-th best score), the rest the best world-wide.""" + ok = np.flatnonzero((score > 0) & land) + k = min(len(ok), int(round(share * land.sum()))) + out = np.zeros(len(score), bool) + if k == 0: + return out + cut = LOCAL_MIN * np.partition(score[ok], len(ok) - k)[len(ok) - k] # the k-th best world-wide + ok_l = ok[score[ok] >= cut] + quota = np.floor(LOCAL * share * np.bincount(prov[land], minlength=prov.max() + 1) + 0.5).astype(np.int64) + o = ok_l[np.lexsort((-score[ok_l], prov[ok_l]))] # by province, best first + p = prov[o] + first = np.r_[0, np.flatnonzero(p[1:] != p[:-1]) + 1] + rank = np.arange(len(o)) - np.repeat(first, np.diff(np.r_[first, len(o)])) + local = o[rank < quota[p]] + out[local[np.argsort(-score[local], kind="stable")[:k]]] = True + rest = ok[~out[ok]] + m = k - int(out.sum()) + if m > 0: + out[rest[np.argpartition(-score[rest], m - 1)[:m]]] = True + return out + + +def rules(f): + """Each deposit's suitability (≥ 0, 0 = impossible) from the geology/climate fields in f.""" + from .crust import CRATON, LIP, POST_OROGEN, PRE_OROGEN, RIFT, VOLCANO + from .environment import (LF_ARC, LF_MASSIF, LF_MOUNTAINS, LI_ANDESITE, LI_BASALT, LI_GRANITE, LI_LIMESTONE, + LI_METAMORPHIC, LI_SANDSTONE, LF_DUNES, LF_PLAIN) + z, lit, age, lf, rel = f["z"], f["lit"], f["age"], f["lf"], f["relief"] + t, p = f["t_mean"], f["p_ann"] + b = lambda m: np.asarray(m, dtype=np.float64) + deep = np.clip((z - 800.0) / 2000.0, 0, 1) * np.clip(rel / 1500.0, 0.3, 1.0) # deep-mountain mining ground + deep = np.maximum(deep, 0.6 * b(np.isin(lf, [LF_MOUNTAINS, LF_MASSIF])) * np.clip(z / 2500.0, 0, 1)) + oro = b(np.isin(age, [PRE_OROGEN, POST_OROGEN])) + craton = b(age == CRATON) + rift = b(age == RIFT) + arc = b((lit == LI_ANDESITE) | (lf == LF_ARC)) + volc = b(age == VOLCANO) + arc + gran = b(lit == LI_GRANITE) + meta = b(lit == LI_METAMORPHIC) + lime = b(lit == LI_LIMESTONE) + sand = b(lit == LI_SANDSTONE) + basalt = b((lit == LI_BASALT) | (age == LIP)) + sed = lime + sand + low = b(z < 800) + hot, wet = b(t > 20), b(p > 1200) + arid = b(p < 250) * b(t > 10) + suture = b(f["d_coll"] < 400) * (oro + meta) + return { + "bog_iron": b(p > 700) * low * (b(lf == LF_PLAIN) + 0.5) + hot * wet * low, + "ironstone": sed * b(z < 1200), + "iron_high": deep * (1.0 + craton + meta + oro), + "coal": b(f["coal"]), + "copper": arc * (1 + deep) + 0.4 * sand * rift + 0.3 * basalt, + "molybdenum": arc * deep, + "tin": (gran + meta) * oro * (0.3 + deep), + "tungsten": (gran + meta) * (oro + 0.3) * deep, + "lead_zinc": lime + 0.3 * sand * (1 + rift), + "silver": arc + 0.3 * lime + 0.3 * deep * oro, + "gold": (oro + craton * 0.7 + meta) * (0.3 + deep) + 0.3 * arc, + "manganese": (sed + craton) * (0.2 + deep), + "chromium": suture * deep + 0.3 * craton * basalt * deep, + "nickel": basalt * hot * wet + (craton + suture) * deep, + "cobalt": (basalt * hot * wet + craton * deep + arc * 0.3) * deep, + "vanadium": (basalt + craton) * deep, + "platinum": craton * (basalt + meta + 0.3) * deep, + "mercury": volc * (0.5 + 0.5 * b(f["t_mean"] > -30)), + "sulfur": volc + 0.5 * b(f["salt"]) + 0.3 * sed * arid, + "saltpetre": arid * b(lf != LF_MOUNTAINS) + 0.3 * lime * b(p > 800), + "bauxite": hot * wet * (gran + basalt + 0.3) * b(rel < 800) + 0.4 * lime * b(t > 14) * b(p > 600), + "fluorite": (gran + lime * 0.4) * (rift + oro + 0.2), + "rock_salt": b(f["salt"]) + 0.6 * b(f["endo"]) + 0.3 * sed * b(p < 500), + "potash": b(f["salt"]) + 0.3 * b(f["endo"]) + 0.2 * sed * b(p < 400), + "phosphate": sed * low * b(f["dist_ocean"] < 400) + 0.3 * lime, + "oil_gas": sed * low * (1 + b(f["d_coll"] < 1200) + rift), + "helium": sed * low * (craton + gran), + "diamond": craton * (1 + deep), + "magnesite": (lime + meta * 0.5) * (0.3 + deep) + 0.3 * b(f["salt"]), + "kaolin": (gran + sand * 0.4) * b(p > 900) * b(t > 8), + "glass_sand": sand * (1 + b(lf == LF_DUNES) + b(f["dist_ocean"] < 150)), + "graphite": meta * (0.3 + deep), + } + + +def place(g, f, seed): + """(deposits uint32 bitmask, deposit_main uint8 legend index) over the grid.""" + land = np.asarray(f["land"], bool) + R = rules(f) + prov = provinces(g, seed) + bits = np.zeros(g.n, np.uint32) + main = np.zeros(g.n, np.uint8) + best = np.full(g.n, np.inf) + for i, (name, _, share) in enumerate(DEPOSITS): + nz = 0.5 + 0.5 * fbm(g.xyz, name_seed(seed, "deposit:" + name), 3, 6.0) # mining districts + sel = _rank_select(R[name] * (0.25 + nz) ** 2, share, land, prov) + bits[sel] |= np.uint32(1 << i) + rarer = sel & (share < best) + main[rarer], best[rarer] = i + 1, share + return bits, main diff --git a/mapgen/newworld.py b/mapgen/newworld.py new file mode 100644 index 0000000..80109b0 --- /dev/null +++ b/mapgen/newworld.py @@ -0,0 +1,117 @@ +"""`mapgen.py --world DIR new-world`: a random starting world (sketch/*.png continents, config/tectonics.toml plates), +so a build runs without drawing anything. Used for the moons; edit or redraw afterwards.""" +from __future__ import annotations + +from pathlib import Path + +import numpy as np +from PIL import Image +from scipy import ndimage + +from .noise import fbm, ridged +from .sphere import latlon_to_xyz + +W, H = 2000, 1000 + + +def _spread(rng, n, min_deg, avoid=(), max_lat=70.0, tries=4000): + pts = list(avoid) + out = [] + for _ in range(tries): + if len(out) == n: + break + lat = np.degrees(np.arcsin(rng.uniform(np.sin(np.radians(-max_lat)), np.sin(np.radians(max_lat))))) + lon = rng.uniform(-180, 180) + p = latlon_to_xyz(lat, lon) + if all(np.degrees(np.arccos(np.clip(p @ latlon_to_xyz(*q), -1, 1))) >= min_deg for q in pts): + pts.append((lat, lon)) + out.append((float(lat), float(lon))) + return out + + +def _centroid(mask, LAT, LON): + w = np.cos(np.radians(LAT)) * mask + p = (latlon_to_xyz(LAT, LON) * w[..., None]).sum(axis=(0, 1)) + p /= np.linalg.norm(p) + return float(np.degrees(np.arcsin(p[2]))), float(np.degrees(np.arctan2(p[1], p[0]))) + + +def make(root: Path, seed: int, continents: int = 6, land_share: float = 0.27, toward=None, ridge: float = 0.86, + force: bool = False) -> None: + root = Path(root) + sk, cfg = root / "sketch", root / "config" + targets = [sk / "land.png", cfg / "tectonics.toml"] + if not force and any(p.exists() for p in targets): + raise SystemExit("new-world: sketch/ or config/tectonics.toml already exist (use --force to overwrite)") + rng = np.random.default_rng(seed) + w, h = W // 2, H // 2 + LAT, LON = np.meshgrid(90.0 - (np.arange(h) + 0.5) * 180.0 / h, -180.0 + (np.arange(w) + 0.5) * 360.0 / w, + indexing="ij") + xyz = latlon_to_xyz(LAT, LON).reshape(-1, 3) + + centres = _spread(rng, continents, 35.0, max_lat=55.0) + field = np.zeros(len(xyz)) + for lat, lon in centres: + r = np.radians(rng.uniform(18.0, 34.0)) + d = np.arccos(np.clip(xyz @ latlon_to_xyz(lat, lon), -1, 1)) + field = np.maximum(field, np.clip(1.0 - d / r, 0.0, None) * rng.uniform(0.8, 1.2)) + field += 0.55 * fbm(xyz, seed + 7, 6, 2.2) + if toward is not None: + field += 1.2 * (xyz @ latlon_to_xyz(*toward)) + area = np.cos(np.radians(LAT)).ravel() + order = np.argsort(-field) + cut = field[order][np.searchsorted(np.cumsum(area[order]) / area.sum(), land_share)] + land = (field > cut).reshape(h, w) + + mount = (ridged(xyz, seed + 11, 5, 3.0) > ridge).reshape(h, w) & ndimage.binary_erosion(land, iterations=6) + + def save(name, m): + Image.fromarray((m * 255).astype(np.uint8), "L").resize((W, H), Image.BILINEAR).point( + lambda v: 255 if v > 127 else 0).save(sk / f"{name}.png") + + sk.mkdir(parents=True, exist_ok=True) + cfg.mkdir(parents=True, exist_ok=True) + save("land", land) + save("mountains", mount) + for name in ("desert", "rainforest", "trench"): + save(name, np.zeros_like(land)) + + lab, n = ndimage.label(land) + for a, b in zip(lab[:, 0], lab[:, -1]): + if a and b and a != b: + lab[lab == b] = a + sizes = sorted(((int((lab == i).sum()), i) for i in np.unique(lab) if i), reverse=True) + plates, seeds = [], [] + for k, (px, i) in enumerate(sizes[:continents]): + m = lab == i + if px > 0.035 * m.size: # a big landmass: two plates, a collision belt between + lat, lon = _centroid(m, LAT, LON) + ang = rng.uniform(0, np.pi) + side = (np.sin(ang) * (LAT - lat) + np.cos(ang) * ((LON - lon + 180) % 360 - 180) * np.cos(np.radians(lat))) > 0 + parts = [m & side, m & ~side] + else: + parts = [m] + for j, part in enumerate(parts): + if part.sum() < 50: + continue + lat, lon = _centroid(part, LAT, LON) + seeds.append((lat, lon)) + plates.append((f"continent-{k + 1}{'ab'[j] if len(parts) > 1 else ''}", lat, lon, "continental", + rng.uniform(0, 360), rng.uniform(2.0, 5.0))) + for k, (lat, lon) in enumerate(_spread(rng, 14 - len(plates) // 2, 28.0, seeds, max_lat=85.0)): + plates.append((f"ocean-{k + 1}", lat, lon, "oceanic", rng.uniform(0, 360), rng.uniform(3.0, 8.0))) + + ocean = np.argwhere(~ndimage.binary_dilation(land, iterations=10)) + def sea_point(): + y, x = ocean[rng.integers(len(ocean))] + return float(LAT[y, x]), float(LON[y, x]) + f = lambda v: f"[{v[0]:.1f}, {v[1]:.1f}]" + t = [f"# Generated by `mapgen.py new-world --seed {seed}`. Edit freely (see README.md → Tectonics).", + "# seed = [lat, lon] where the plate grows from; motion = [azimuth° clockwise from north, speed cm/yr].", ""] + for pid, lat, lon, kind, az, sp in plates: + t += ["[[plate]]", f'id = "{pid}"', f"seed = {f((lat, lon))}", f'kind = "{kind}"', f"motion = [{az:.0f}.0, {sp:.1f}]", ""] + if len(ocean): + t += ["[[hotspot]]", 'name = "island-chain"', f"center = {f(sea_point())}", "length_km = 600.0", ""] + (cfg / "tectonics.toml").write_text("\n".join(t)) + + print(f"new-world: {len(plates)} plates, sketch/ and config/tectonics.toml written (seed {seed})") diff --git a/mapgen/noise.py b/mapgen/noise.py new file mode 100644 index 0000000..e52b7f8 --- /dev/null +++ b/mapgen/noise.py @@ -0,0 +1,111 @@ +"""Seeded, vectorized 3D value noise and fractal sums on points (N, 3).""" +from __future__ import annotations + +import zlib + +import numpy as np + +_MASK = (1 << 64) - 1 +_K1 = np.uint64(0x9E3779B185EBCA87) +_K2 = np.uint64(0xC2B2AE3D27D4EB4F) +_K3 = np.uint64(0x165667B19E3779F9) +_M1 = np.uint64(0x94D049BB133111EB) + + +def _hash(ix, iy, iz, seed: int): + s = np.uint64((seed * 0x27D4EB2F165667C5 + 0x632BE59BD9B4E019) & _MASK) + h = ix.astype(np.uint64) * _K1 ^ iy.astype(np.uint64) * _K2 ^ iz.astype(np.uint64) * _K3 ^ s + h ^= h >> np.uint64(31) + h *= _M1 + h ^= h >> np.uint64(29) + return (h >> np.uint64(11)).astype(np.float64) / float(1 << 53) * 2.0 - 1.0 + + +def value_noise(p, seed: int): + if _noise_jit is not None: + p = np.asarray(p) + if p.dtype == np.float64 and p.ndim == 2 and p.shape[1] == 3: + s = np.uint64((seed * 0x27D4EB2F165667C5 + 0x632BE59BD9B4E019) & _MASK) + return _noise_jit(np.ascontiguousarray(p), s) + return _value_noise(p, seed) + + +def _value_noise(p, seed: int): + f = np.floor(p) + i = f.astype(np.int64) + t = p - f + t = t * t * (3.0 - 2.0 * t) + out = np.zeros(len(p)) + for dx in (0, 1): + wx = t[:, 0] if dx else 1.0 - t[:, 0] + for dy in (0, 1): + wy = t[:, 1] if dy else 1.0 - t[:, 1] + for dz in (0, 1): + wz = t[:, 2] if dz else 1.0 - t[:, 2] + out += wx * wy * wz * _hash(i[:, 0] + dx, i[:, 1] + dy, i[:, 2] + dz, seed) + return out + + +def _make_noise_jit(): + """value_noise compiled (numba optional): per point the same integer hash and the same float steps in the same + order, so the same values; one call costs microseconds instead of a dozen numpy passes (river meanders call it + on a few points at a time). WORLDGEN_NO_JIT=1 turns it off.""" + import os + if os.environ.get("WORLDGEN_NO_JIT"): + return None + try: + import numba + except ImportError: + return None + K1, K2, K3, M1 = _K1, _K2, _K3, _M1 + + @numba.njit(cache=True) + def noise(p, s): + n = p.shape[0] + out = np.empty(n) + for k in range(n): + fx, fy, fz = np.floor(p[k, 0]), np.floor(p[k, 1]), np.floor(p[k, 2]) + ix, iy, iz = np.int64(fx), np.int64(fy), np.int64(fz) + tx, ty, tz = p[k, 0] - fx, p[k, 1] - fy, p[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) + acc = 0.0 + for dx in range(2): + wx = tx if dx else 1.0 - tx + for dy in range(2): + wy = ty if dy else 1.0 - ty + for dz in range(2): + wz = tz if dz else 1.0 - tz + h = (np.uint64(ix + dx) * K1) ^ (np.uint64(iy + dy) * K2) ^ (np.uint64(iz + dz) * K3) ^ s + h ^= h >> np.uint64(31) + h *= M1 + h ^= h >> np.uint64(29) + v = np.float64(h >> np.uint64(11)) / 9007199254740992.0 * 2.0 - 1.0 + acc += wx * wy * wz * v + out[k] = acc + return out + return noise + + +_noise_jit = _make_noise_jit() + + +def fbm(xyz, seed: int, octaves: int = 5, freq: float = 2.0, lacunarity: float = 2.0, gain: float = 0.5): + total = np.zeros(len(xyz)) + amp, norm, f = 1.0, 0.0, freq + for o in range(octaves): + total += amp * value_noise(xyz * f, seed + 1013 * o) + norm += amp + amp *= gain + f *= lacunarity + return total / norm + + +def ridged(xyz, seed: int, octaves: int = 5, freq: float = 2.0): + return 1.0 - np.abs(fbm(xyz, seed, octaves, freq)) + + +def name_seed(seed: int, name: str) -> int: + """A per-feature noise seed: the build seed plus a hash of the feature's name (stable across runs).""" + return seed + zlib.crc32(name.encode()) % 100000 diff --git a/mapgen/ocean.py b/mapgen/ocean.py new file mode 100644 index 0000000..67dfaf7 --- /dev/null +++ b/mapgen/ocean.py @@ -0,0 +1,166 @@ +"""Ocean circulation for stage `climate`: wind-driven surface currents, sea-surface temperature, Ekman upwelling +and a productivity index. + +Currents: the Stommel model on the sphere, r∇²ψ + βψ_x = k·curl(τ/ρH), solved for the stream function ψ over +every cell; land is the same fluid with `land_friction`× the friction (Brinkman penalisation), so coasts block the +flow and islands need no special treatment. u = k × ∇ψ. Western boundary currents come out ≈ r/β wide. +""" +from __future__ import annotations + +import numpy as np +from scipy import sparse +from scipy.sparse import linalg as splinalg + +from .graph import bicgstab_jacobi, distance_to, gradient, pmap, smooth_km +from .grid import rowdot_at +from .sphere import east_north + +DEFAULTS = { + "enabled": True, "friction_days": 5.0, "layer_m": 150.0, "stress_k": 1.0, "land_friction": 1000.0, + "rho_air": 1.2, "drag": 1.3e-3, "direct_max_cells": 500000, + "relax_days": 300.0, "kappa_m2s": 1000.0, + "ekman_min_lat": 3.0, "upwell_coast_km": 100.0, + "prod_upwell": 0.6, "prod_shelf": 0.4, "prod_mix": 0.3, "upwell_ref_m_yr": 100.0, + "prod_front": 0.4, "front_min": 0.5, "front_ref": 1.0, +} +RHO_W = 1025.0 +YEAR_S = 3.15576e7 + + +def omega(day_hours): + return 2.0 * np.pi / (float(day_hours) * 3600.0) + + +def divergence(g, V): + """Graph divergence of a tangent field (V per km): (2/k_i) Σ_j V_j·t_ij / d_ij (exact for linear V on a hex + lattice; the V_i terms cancel).""" + w = rowdot_at(V, g.dst, g.edge_tangents) / g.edge_km + return np.bincount(g.src, weights=w, minlength=g.n) * (2.0 / g.counts) + + +def _edge_operator(g, w): + """Sparse operator f ↦ Σ_j w_ij (f_j − f_i).""" + A = sparse.csr_matrix((w, (g.src, g.dst)), shape=(g.n, g.n)) + return A - sparse.diags(np.asarray(A.sum(axis=1)).ravel()) + + +def wind_stress(wind, P): + """Bulk formula τ = k·ρ_air·C_d·|w|·w (N/m²).""" + wind = np.asarray(wind, dtype=np.float64) + return P["stress_k"] * P["rho_air"] * P["drag"] * np.linalg.norm(wind, axis=1, keepdims=True) * wind + + +def streamfunction(g, ocean, wind, P, day_hours): + """ψ (m²/s) of the wind-driven surface flow; one direct sparse solve over every cell.""" + Om = omega(day_hours) + R = g.radius_km * 1e3 + beta = 2.0 * Om * np.cos(np.radians(g.lat)) / R # 1/(m·s) + beta30 = 2.0 * Om * np.cos(np.radians(30.0)) / R + r = max(1.0 / (P["friction_days"] * 86400.0), beta30 * g.spacing_km * 1e3) # boundary layer ≥ one cell + F = np.where(ocean[:, None], wind_stress(wind, P), 0.0) / (RHO_W * P["layer_m"]) + curl = divergence(g, np.cross(F, g.xyz)) * 1e-3 # k·curl F = ∇·(F×k), 1/s² + fr = np.where(ocean, 1.0, P["land_friction"]) + fe = 2.0 / (1.0 / fr[g.src] + 1.0 / fr[g.dst]) # harmonic mean: coasts count as land + L = _edge_operator(g, 4.0 / (g.counts[g.src] * g.edge_km ** 2) * fe) # ∇·(fr∇), per km² + e, _ = east_north(g.xyz) + Dx = _edge_operator(g, (2.0 / g.counts[g.src]) * rowdot_at(e, g.src, g.edge_tangents) / g.edge_km) + A = (L + sparse.diags(beta / r * 1e3) @ Dx).tolil() # β/r·1e3: per km + b = 1e6 * curl / r + # gauge: ψ is defined up to a constant; pin it deep inland, where land friction soaks up the solve's + # compatibility residual (pinned at sea it would leave a point vortex there) + k = int(np.argmax(distance_to(g, ocean))) if (~ocean).any() else 0 + A[k, :] = 0.0 + A[k, k] = 1.0 + b[k] = 0.0 + return splinalg.spsolve(A.tocsc(), b) + + +def velocity(g, psi): + """u = k × ∇ψ (m/s), tangent 3-vectors.""" + return np.cross(g.xyz, gradient(g, psi) * 1e-3) + + +def _coarse_currents(g, ocean, wind, P, day_hours): + """Grids too big for a direct solve: solve on the parent H3 resolution, carry u down, smooth.""" + from .grid import build_grid, cell_parents + cg = build_grid(g.res - 1, g.radius_km) + parent = np.searchsorted(cg.ids, cell_parents(g.ids, g.res - 1)) + cnt = np.maximum(np.bincount(parent, minlength=cg.n), 1) + c_ocean = np.bincount(parent, weights=ocean.astype(float), minlength=cg.n) / cnt > 0.5 + c_wind = np.stack([np.bincount(parent, weights=wind[:, k], minlength=cg.n) / cnt for k in range(3)], axis=1) + cu = currents(cg, c_ocean, c_wind, P, day_hours) + u = np.stack(pmap(lambda k: smooth_km(g, cu[parent, k], cg.spacing_km), range(3)), axis=1) + return u - np.sum(u * g.xyz, axis=1, keepdims=True) * g.xyz # back into the tangent plane + + +def currents(g, ocean, wind, P, day_hours): + """Surface current (n,3) m/s; 0 on land.""" + ocean = np.asarray(ocean, bool) + wind = np.asarray(wind, dtype=np.float64) + if not ocean.any(): + return np.zeros((g.n, 3)) + if g.n > P["direct_max_cells"]: + u = _coarse_currents(g, ocean, wind, P, day_hours) + else: + u = velocity(g, streamfunction(g, ocean, wind, P, day_hours)) + return np.where(ocean[:, None], u, 0.0) + + +def _masked_edges(g, mask, w): + return _edge_operator(g, np.where(mask[g.src] & mask[g.dst], w, 0.0)) + + +def sst(g, ocean, u, T_eq, P): + """Annual-mean sea-surface temperature (°C): steady u·∇T − κ∇²T = λ(T_eq − T) over the ocean (upwind + advection, no flux into land); land keeps T_eq.""" + ocean = np.asarray(ocean, bool) + T_eq = np.asarray(T_eq, dtype=np.float64) + lam = 1.0 / (P["relax_days"] * 86400.0) + up = (4.0 / g.counts[g.src]) * np.maximum(-rowdot_at(u, g.src, g.edge_tangents), 0.0) / (g.edge_km * 1e3) + adv = -_masked_edges(g, ocean, up) # Σ_upstream c_ij (T_i − T_j), 1/s + dif = _masked_edges(g, ocean, 4.0 / (g.counts[g.src] * g.edge_km ** 2)) * (P["kappa_m2s"] * 1e-6) + A = (adv - dif) / lam + sparse.identity(g.n) + A = (sparse.diags(ocean.astype(float)) @ A + sparse.diags((~ocean).astype(float))).tocsr() + T, info = bicgstab_jacobi(A, T_eq, T_eq.copy(), 1.0 / A.diagonal(), 1e-9, 5000) + if info != 0: + raise ValueError(f"sst: solver did not converge (info={info})") + return T + + +def upwelling(g, ocean, wind, P, day_hours): + """Ekman pumping (m/yr, + up): w = ∇·M, M = τ×k/(ρf), |f| floored at `ekman_min_lat`. Transport pointing off a + coast leaves the coast cell (land carries none), so coastal upwelling needs no separate rule. Smoothed over + `upwell_coast_km`; 0 on land.""" + ocean = np.asarray(ocean, bool) + Om = omega(day_hours) + fmin = 2.0 * Om * np.sin(np.radians(P["ekman_min_lat"])) + f = np.where(g.lat >= 0, 1.0, -1.0) * np.maximum(np.abs(2.0 * Om * np.sin(np.radians(g.lat))), fmin) + tau = np.where(ocean[:, None], wind_stress(wind, P), 0.0) + M = np.cross(tau, g.xyz) / (RHO_W * f[:, None]) # m²/s + w = np.where(ocean, divergence(g, M) * 1e-3 * YEAR_S, 0.0) + return np.where(ocean, smooth_km(g, w, P["upwell_coast_km"]), 0.0) + + +def sst_front(g, ocean, sst_c): + """|∇SST| over the sea (°C per 100 km); edges to land count as flat, so coasts are no front.""" + sst_c = np.asarray(sst_c, dtype=np.float64) + both = ocean[g.src] & ocean[g.dst] + w = np.where(both, sst_c[g.dst] - sst_c[g.src], 0.0) / g.edge_km + grad = np.stack([np.bincount(g.src, weights=w * g.edge_tangents[:, c], minlength=g.n) for c in range(3)], axis=1) + return np.where(ocean, np.linalg.norm(grad, axis=1) * (2.0 / g.counts) * 100.0, 0.0) + + +def productivity(g, ocean, w, z, t_range, sst_c, P): + """0–1 sea productivity: upwelling (saturating), shallow shelf, winter mixing and SST fronts (where warm and cold + currents meet, e.g. a Brazil–Malvinas confluence; gradients under `front_min` °C/100 km add nothing), dimmed + toward the poles.""" + ocean = np.asarray(ocean, bool) + depth = np.maximum(-np.asarray(z, dtype=np.float64), 0.0) + x = np.maximum(w, 0.0) / P["upwell_ref_m_yr"] + shelf = np.clip((1000.0 - depth) / 800.0, 0.0, 1.0) + mix = np.clip(np.asarray(t_range) / 20.0, 0.0, 1.0) * np.clip((20.0 - np.asarray(sst_c)) / 20.0, 0.0, 1.0) + light = 0.3 + 0.7 * np.cos(np.radians(g.lat)) + xf = np.maximum(sst_front(g, ocean, sst_c) - P["front_min"], 0.0) / P["front_ref"] + p = (P["prod_upwell"] * x / (1.0 + x) + P["prod_shelf"] * shelf + P["prod_mix"] * mix + + P["prod_front"] * xf / (1.0 + xf)) * light + return np.where(ocean, np.clip(p, 0.0, 1.0), 0.0) diff --git a/mapgen/pipeline.py b/mapgen/pipeline.py new file mode 100644 index 0000000..fca2fa1 --- /dev/null +++ b/mapgen/pipeline.py @@ -0,0 +1,180 @@ +"""Stage runner with an input-hash cache.""" +from __future__ import annotations + +import hashlib +import importlib +import os +import time +from dataclasses import dataclass, field +from pathlib import Path + +import numpy as np + +from . import config as C +from .grid import Grid + +STAGES = ["grid", "sketch", "plates", "crust", "elevation", "erosion", "climate", + "hydrology", "seabed", "environment", "ice", "fields", "render"] + + +class StageError(RuntimeError): + pass + + +@dataclass +class Ctx: + root: Path + cfg: dict + tect: dict + res: int + data: dict = field(default_factory=dict) + grid: Grid | None = None + out_dir: Path | None = None # render: out/r<res> unless set (an era writes out/r<res>/eras/<era>) + preview_dir: Path | None = None + low_memory: bool = False # trade speed for a lower memory peak; never changes the results + + @property + def seed(self) -> int: + return int(self.cfg["build"]["seed"]) + + def need(self, *keys): + missing = [k for k in keys if k not in self.data] + if missing: + raise StageError(f"missing inputs {missing}: run the stage that produces them first") + return [self.data[k] for k in keys] + + +def inputs_key(root: Path, res: int) -> str: + h = hashlib.sha256(str(res).encode()) + files = [] + for sub, pat in (("config", "*.toml"), ("masks", "*.png"), ("sketch", "*.png")): + files += sorted(p for p in (root / sub).glob(pat) if p.name != C.ERAS_FILE) # eras: their own key + files += sorted(Path(__file__).resolve().parent.glob("*.py")) + for p in files: + h.update(p.name.encode()) + h.update(p.read_bytes()) + return h.hexdigest()[:16] + + +ALIGN = 64 + + +def save_npz_aligned(path, /, **arrays) -> None: + """np.savez(path, **arrays), but each array's data starts at a multiple of ALIGN bytes in the file (the local zip + header gets a padding extra field, as zipalign does): an ordinary .npz for np.load, whose arrays npz_maps can + map in place. Alignment matters for more than speed: numpy sums 8-byte-misaligned data in buffered chunks, + i.e. in another order — mapped misaligned inputs would change the last bits of results.""" + import io + import struct + import zipfile + from numpy.lib import format as F + with zipfile.ZipFile(path, "w", compression=zipfile.ZIP_STORED, allowZip64=True) as zf: + for name, v in arrays.items(): + v = np.asarray(v) + if v.dtype.hasobject: + raise ValueError(f"save_npz_aligned: {name} has object dtype") + head = io.BytesIO() + d = F.header_data_from_array_1_0(v) + try: + F.write_array_header_1_0(head, d) + except ValueError: + head = io.BytesIO() + F.write_array_header_2_0(head, d) + info = zipfile.ZipInfo(f"{name}.npy", date_time=(1980, 1, 1, 0, 0, 0)) + info.compress_type = zipfile.ZIP_STORED + start = zf.fp.tell() + 30 + len(info.filename.encode()) + 20 + len(head.getvalue()) # 20: zip64 field + pad = -start % ALIGN + if 0 < pad < 4: # an extra field is at least its 4-byte header + pad += ALIGN + if pad: + info.extra = struct.pack("<HH", 0xA1A1, pad - 4) + bytes(pad - 4) + with zf.open(info, "w", force_zip64=True) as m: + m.write(head.getvalue()) + _write_data(m, v) + + +def _write_data(m, v) -> None: + """The array's bytes (C order, or Fortran order for an F-contiguous array, as np.save) in 16 MB pieces.""" + flat = v.T.reshape(-1) if (v.flags.f_contiguous and not v.flags.c_contiguous) else np.ascontiguousarray(v).reshape(-1) + step = max(1, (16 << 20) // max(v.itemsize, 1)) + for i in range(0, flat.size, step): + m.write(flat[i:i + step].tobytes()) + + +def npz_maps(f: Path) -> dict: + """The arrays of an uncompressed .npz (np.savez) mapped from the file, copy-on-write: plain writable arrays with + the same values, whose pages the OS reads on use and can drop again (writes stay private, the file never + changes). Members that can't be mapped (compressed, object dtype) are read into memory as np.load would.""" + import mmap + import struct + import zipfile + out = {} + with open(f, "rb") as fh, zipfile.ZipFile(fh) as zf: + mm = mmap.mmap(fh.fileno(), 0, access=mmap.ACCESS_COPY) + for info in zf.infolist(): + name = info.filename[:-4] if info.filename.endswith(".npy") else info.filename + ok = info.compress_type == zipfile.ZIP_STORED + if ok: + fh.seek(info.header_offset) + local = fh.read(30) + n_name, n_extra = struct.unpack("<HH", local[26:30]) + fh.seek(info.header_offset + 30 + n_name + n_extra) + version = np.lib.format.read_magic(fh) + read_header = {(1, 0): np.lib.format.read_array_header_1_0, + (2, 0): np.lib.format.read_array_header_2_0}.get(version) + ok = read_header is not None + if ok: + shape, fortran, dtype = read_header(fh, max_header_size=1 << 20) + ok = not dtype.hasobject + ok = fh.tell() % ALIGN == 0 # misaligned: read (see save_npz_aligned) + if ok: + count = int(np.prod(shape, dtype=np.int64)) + a = np.frombuffer(mm, dtype=dtype, count=count, offset=fh.tell()) if count else np.empty(0, dtype) + out[name] = a.reshape(shape, order="F" if fortran else "C") + else: + with zf.open(info) as m: + out[name] = np.lib.format.read_array(m, allow_pickle=False) + return out + + +def _after(ctx: Ctx, name: str) -> None: + if name == "grid": + ctx.grid = Grid.from_arrays(ctx.data, ctx.res, float(ctx.cfg["planet"]["radius_km"])) + + +def build(root: Path, res: int, start: str | None = None, stop: str | None = None, log=print, + low_memory: bool = False) -> Ctx: + cfg, tect = C.load(root) + ctx = Ctx(root, cfg, tect, res, low_memory=low_memory) + cache = root / "out" / "cache" / f"r{res}" + cache.mkdir(parents=True, exist_ok=True) + key = inputs_key(root, res) + stages = STAGES[: STAGES.index(stop) + 1] if stop else STAGES + forced = False + for name in stages: + forced = forced or name == start + f = cache / f"{name}.npz" + t0 = time.time() + if not forced and f.exists(): + with np.load(f, allow_pickle=False) as z: + hit = str(z["_key"]) == key + if hit and not ctx.low_memory: + ctx.data.update({k: z[k] for k in z.files if k != "_key"}) + if hit: + if ctx.low_memory: + ctx.data.update({k: v for k, v in npz_maps(f).items() if k != "_key"}) + _after(ctx, name) + log(f"{name}: cached") + continue + forced = True + out = importlib.import_module(f"mapgen.{name}").run(ctx) + tmp = f.with_name(f"{f.stem}.tmp-{os.getpid()}.npz") + save_npz_aligned(tmp, _key=np.array(key), **out) + os.replace(tmp, f) # never truncated in place: maps of the old file stay valid + if ctx.low_memory: # fields kept on disk from here on (the cache file just written) + maps = npz_maps(f) + out = {k: maps.get(k, v) if isinstance(v, np.ndarray) else v for k, v in out.items()} + ctx.data.update(out) + _after(ctx, name) + log(f"{name}: {time.time() - t0:.1f}s") + return ctx diff --git a/mapgen/plateaus.py b/mapgen/plateaus.py new file mode 100644 index 0000000..9647d04 --- /dev/null +++ b/mapgen/plateaus.py @@ -0,0 +1,179 @@ +"""Sunken plateaus: Kerguelen-type continental crust that never rose. Outline, +surface and volcanic field are functions of position and the plateau's config (deterministic per name), so the world +build, the seabed pass and the viewer agree at any resolution.""" +from __future__ import annotations + +import numpy as np + +from .graph import distance_to +from .noise import fbm, name_seed +from .sphere import east_north, gc_dist_km, great_circle_point, latlon_to_xyz, tangent_dir, xyz_to_latlon +from .zones import smootherstep + +EDGE_WARP = 0.15 # outline radius ± 15 % (fractal noise): bays and lobes, like a small continent +MARGIN_KM = 150.0 # the surface blends down to the surrounding sea floor over this distance inside the outline +RELIEF_M = 400.0 # internal relief (ridges, basins) +HIDDEN_MAX_M = -550.0 # hidden plateaus: every cell at least this deep (checked at −500 m after erosion's re-solve) +SITE_RULES = (("land", 400.0), ("a plate boundary", 150.0), ("the trench", 200.0)) # km from the outline + + +def semi_axes(p: dict): + """(long, short) semi-axes (km) of the plateau's ellipse: π·a·b = area_km2, a / b = elongation.""" + e = float(p.get("elongation", 1.0)) + b = float(np.sqrt(p["area_km2"] / (np.pi * e))) + return e * b, b + + +def _near(xyz, p, radius_km, pad_km=0.0): + a, _ = semi_axes(p) + lim = min(np.pi, (a * (1.0 + EDGE_WARP) / (1.0 - EDGE_WARP) + pad_km) / radius_km) + return xyz @ latlon_to_xyz(*p["center"]) > np.cos(lim) + + +def rho(xyz, p: dict, seed: int, radius_km: float): + """Normalised radius of points: 0 at the centre, < 1 inside the noise-warped elliptical outline.""" + xyz = np.asarray(xyz, dtype=np.float64) + c = latlon_to_xyz(*p["center"]) + e, n = east_north(c[None]) + ang = np.arccos(np.clip(xyz @ c, -1.0, 1.0)) * radius_km # km from the centre along the surface + t = tangent_dir(np.broadcast_to(c, xyz.shape), xyz) + x, y = ang * (t @ e[0]), ang * (t @ n[0]) # east, north (azimuthal equidistant) + az = np.radians(p.get("azimuth_deg", 0.0)) + along, across = x * np.sin(az) + y * np.cos(az), x * np.cos(az) - y * np.sin(az) + a, b = semi_axes(p) + r = np.sqrt((along / a) ** 2 + (across / b) ** 2) + w = np.clip(2.0 * fbm(xyz, name_seed(seed, p["name"]), 4, 25.0), -1.0, 1.0) + return r / (1.0 + EDGE_WARP * w) + + +def cell_ids(xyz, plateaus: list, seed: int, radius_km: float): + """Plateau index per point (−1 outside every plateau).""" + xyz = np.asarray(xyz, dtype=np.float64) + ids = np.full(len(xyz), -1, np.int16) + for k, p in enumerate(plateaus): + idx = np.flatnonzero(_near(xyz, p, radius_km)) + if len(idx): + ids[idx[rho(xyz[idx], p, seed, radius_km) < 1.0]] = k + return ids + + +def _place(rng, p, seed, radius_km, origin, max_km, inside=0.85, tries=30): + """A random point within max_km of origin inside the outline (rho < inside); None if none is found.""" + e, n = east_north(origin[None]) + for _ in range(tries): + az, dist = rng.uniform(0.0, 2.0 * np.pi), rng.uniform(0.0, max_km) + q = great_circle_point(origin, np.sin(az) * e[0] + np.cos(az) * n[0], dist, radius_km) + if rho(q[None], p, seed, radius_km)[0] < inside: + return q + return None + + +def _ll(q): + lat, lon = xyz_to_latlon(np.asarray(q, dtype=np.float64)) + return round(float(lat), 4), round(float(lon), 4) + + +def features(p: dict, seed: int, radius_km: float) -> dict: + """The plateau's volcanic field at world scale (deterministic per name): its own hotspot point, 6–20 cones and + 1–3 calderas (none when vent = 0); island plateaus lift 3–6 of the cones to +0.6…+1.5 km.""" + rng = np.random.default_rng(name_seed(seed, p["name"]) + 7) + c = latlon_to_xyz(*p["center"]) + _, b = semi_axes(p) + vent = float(p.get("vent", 1.0)) + hot = _place(rng, p, seed, radius_km, c, 0.3 * b) + hot = c if hot is None else hot + cones, calderas = [], [] + if vent > 0: + for _ in range(int(rng.integers(6, 21))): + q = _place(rng, p, seed, radius_km, hot, 0.6 * b) + r_km, h = float(rng.uniform(15.0, 40.0)), float(rng.uniform(500.0, 2000.0)) * vent + if q is not None: + lat, lon = _ll(q) + cones.append({"lat": lat, "lon": lon, "radius_km": r_km, "height_m": h, "island": False}) + for _ in range(int(rng.integers(1, 4))): + q = _place(rng, p, seed, radius_km, hot, 0.5 * b) + r_km = float(rng.uniform(20.0, 50.0)) + if q is not None: + lat, lon = _ll(q) + calderas.append({"lat": lat, "lon": lon, "radius_km": r_km, "rim_m": 300.0 * vent, + "floor_m": -400.0 * vent}) + if p.get("islands", False) and cones: + k = min(len(cones), int(rng.integers(3, 7))) + for i in rng.choice(len(cones), size=k, replace=False): + cones[int(i)].update(island=True, peak_m=float(rng.uniform(600.0, 1500.0)), + radius_km=float(rng.uniform(40.0, 70.0))) + return {"name": p["name"], "center": [float(v) for v in p["center"]], "hotspot": list(_ll(hot)), + "cones": cones, "calderas": calderas} + + +def surface(xyz, p: dict, feat: dict, seed: int, radius_km: float, z_floor): + """Heights of points inside the outline: the plateau top (depth across top_m by low-frequency noise, ± RELIEF_M) + plus its cones and calderas, blended down to z_floor over MARGIN_KM inside the edge.""" + xyz = np.asarray(xyz, dtype=np.float64) + s = name_seed(seed, p["name"]) + t = np.clip(0.5 + 1.5 * fbm(xyz, s + 1, 3, 8.0), 0.0, 1.0) + top0, top1 = p["top_m"] + z = -(top0 + (top1 - top0) * t) + RELIEF_M * np.clip(2.5 * fbm(xyz, s + 2, 5, 40.0), -1.0, 1.0) + for c in feat["cones"]: + if not c["island"]: + d = gc_dist_km(xyz, latlon_to_xyz(c["lat"], c["lon"]), radius_km) + z = z + c["height_m"] * np.clip(1.0 - d / c["radius_km"], 0.0, 1.0) ** 1.5 + for c in feat["calderas"]: + d = gc_dist_km(xyz, latlon_to_xyz(c["lat"], c["lon"]), radius_km) + r = c["radius_km"] + z = z + c["rim_m"] * np.exp(-((d - r) / (0.25 * r)) ** 2) + c["floor_m"] * (1.0 - smootherstep(d / r)) + _, b = semi_axes(p) + w = smootherstep((1.0 - rho(xyz, p, seed, radius_km)) * b / MARGIN_KM) + return np.asarray(z_floor, dtype=np.float64) + (z - z_floor) * w + + +def islands(xyz, feat: dict, radius_km: float, z): + """Island cones lift the ground to their peak_m (small volcanic islands); only raises.""" + z = np.asarray(z, dtype=np.float64) + for c in feat["cones"]: + if c["island"]: + d = gc_dist_km(np.asarray(xyz, dtype=np.float64), latlon_to_xyz(c["lat"], c["lon"]), radius_km) + f = np.clip(d / c["radius_km"], 0.0, 1.0) + z = np.where(d < c["radius_km"], np.maximum(z, c["peak_m"] - (c["peak_m"] - z) * f ** 1.2), z) + return z + + +def apply(g, z, plateaus: list, plateau_id, seed: int): + """Grid heights with the plateaus set (after the world's sea-level solve): surfaces and volcanic fields; hidden + plateaus clamped to HIDDEN_MAX_M; islands, each island's nearest cell raised to its peak (so every resolution + keeps it).""" + z = np.asarray(z, dtype=np.float64).copy() + R = g.radius_km + for k, p in enumerate(plateaus): + idx = np.flatnonzero(np.asarray(plateau_id) == k) + if len(idx) == 0: + continue + feat = features(p, seed, R) + zk = surface(g.xyz[idx], p, feat, seed, R, z[idx]) + if p.get("islands", False): + z[idx] = islands(g.xyz[idx], feat, R, zk) + for c in feat["cones"]: + if c["island"]: + i = g.cell_index(c["lat"], c["lon"]) + z[i] = max(z[i], c["peak_m"]) + else: + z[idx] = np.minimum(zk, HIDDEN_MAX_M) + return z + + +def site_report(g, data: dict, plateaus: list) -> list: + """Plateaus that no longer fit their site on this world, as warnings: outline ≥ 400 km from land, + ≥ 150 km from plate boundaries, ≥ 200 km from the sketch's trench. Reported, never moved.""" + ids = np.asarray(data["plateau_id"]) + masks = {"land": ~np.asarray(data["ocean"]) & (ids < 0), "a plate boundary": np.asarray(data["bnd_type"]) > 0, + "the trench": np.asarray(data.get("sk_trench", np.zeros(g.n))) > 0.5} + out = [] + for what, lim in SITE_RULES: + if not masks[what].any(): + continue + d = distance_to(g, masks[what]) + for k, p in enumerate(plateaus): + m = ids == k + if m.any() and float(d[m].min()) < lim: + out.append(f"plateau {p['name']}: {float(d[m].min()):.0f} km from {what} (rule ≥ {lim:.0f} km)") + return out diff --git a/mapgen/plates.py b/mapgen/plates.py new file mode 100644 index 0000000..af798bf --- /dev/null +++ b/mapgen/plates.py @@ -0,0 +1,86 @@ +"""Stage `plates`: grow plates from seeds, assign Euler motions, classify boundaries.""" +from __future__ import annotations + +import numpy as np +from scipy.spatial import cKDTree + +from .config import params +from .graph import nearest_source +from .grid import edge_blocks +from .noise import fbm +from .pipeline import StageError +from .sphere import motion_to_omega, velocity + +NONE, CONV, DIV, TRANS = 0, 1, 2, 3 +DEFAULTS = {"land_cost": 0.25, "noise": 0.35, "noise_freq": 3.0, "warp_km": 1000.0, "warp_freq": 1.5, + "min_rate_m_yr": 0.005} + + +def edge_convergence(g, vel): + """Per directed edge s→d: closing speed (m/yr, >0 converging) and tangential slip speed.""" + t, src, dst = g.edge_tangents, g.src, g.dst + along, tang = np.empty(len(dst)), np.empty(len(dst)) + for s in edge_blocks(len(dst)): # per edge: the same values, small temporaries + rel = vel[dst[s]] - vel[src[s]] + along[s] = np.sum(rel * t[s], axis=1) + tang[s] = np.linalg.norm(rel - along[s][:, None] * t[s], axis=1) + return -along, tang + + +def classify_boundaries(g, plate, vel, min_rate): + conv, tang = edge_convergence(g, vel) + s, d = g.src, g.dst + b = plate[s] != plate[d] + cnt = np.bincount(s[b], minlength=g.n) + c_mean = np.bincount(s[b], weights=conv[b], minlength=g.n) / np.maximum(cnt, 1) + t_mean = np.bincount(s[b], weights=tang[b], minlength=g.n) / np.maximum(cnt, 1) + other = np.full(g.n, -1, np.int16) + other[s[b]] = plate[d[b]] + typ = np.zeros(g.n, np.int8) + isb = cnt > 0 + typ[isb] = TRANS + typ[isb & (c_mean > min_rate)] = CONV + typ[isb & (c_mean < -min_rate)] = DIV + rate = np.where(typ == CONV, c_mean, np.where(typ == DIV, -c_mean, np.where(isb, t_mean, 0.0))) + return typ, rate.astype(np.float32), other + + +def _warp_labels(g, plate, seeds, P, seed, continental_kind): + """Domain warp: each cell takes the label found at a noise-displaced point → meandering boundaries. + Only same-kind swaps (cont↔cont, ocean↔ocean): ocean–continent margins stay where growth put them.""" + if P["warp_km"] <= 0: + return plate + amp = P["warp_km"] / g.radius_km + disp = np.stack([fbm(g.xyz, seed + 101 + k, 3, P["warp_freq"]) for k in range(3)], axis=1) + disp /= max(float(disp.std()), 1e-12) # warp_km = RMS displacement per component + q = g.xyz + amp * disp + q /= np.linalg.norm(q, axis=1, keepdims=True) + _, j = cKDTree(g.xyz).query(q) + out = plate[j].astype(np.int16) + keep = continental_kind[out] != continental_kind[plate] + out[keep] = plate[keep] + out[seeds] = np.arange(len(seeds), dtype=np.int16) + return out + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "plates", DEFAULTS) + land, land_hint = ctx.need("sk_land", "m_land_hint") + plates = ctx.tect["plate"] + seeds = np.array([g.cell_index(*p["seed"]) for p in plates], dtype=np.int64) + if len(np.unique(seeds)) < len(seeds): + raise StageError("two plate seeds fall in the same cell; move one in tectonics.toml") + hint = np.clip(land + 0.5 * land_hint, 0.0, 1.0) + nz = fbm(g.xyz, ctx.seed + 11, octaves=4, freq=P["noise_freq"]) + mult = np.where(hint[g.dst] > 0.5, P["land_cost"], 1.0) * (1.0 + P["noise"] * nz[g.dst]) + _, src = nearest_source(g, seeds, g.edge_km * np.maximum(mult, 0.05)) + order = np.argsort(seeds) + plate = order[np.searchsorted(seeds[order], src)].astype(np.int16) + kinds = np.array([p["kind"] == "continental" for p in plates]) + plate = _warp_labels(g, plate, seeds, P, ctx.seed, kinds) + omegas = np.array([motion_to_omega(*p["seed"], *p["motion"], g.radius_km) for p in plates]) + vel = velocity(g.xyz, omegas[plate], g.radius_km) + btype, brate, bother = classify_boundaries(g, plate, vel, P["min_rate_m_yr"]) + return {"plate": plate, "vel": vel, "plate_continental": kinds, + "bnd_type": btype, "bnd_rate": brate, "bnd_other": bother} diff --git a/mapgen/projections.py b/mapgen/projections.py new file mode 100644 index 0000000..67efb19 --- /dev/null +++ b/mapgen/projections.py @@ -0,0 +1,150 @@ +"""Map projections of the equirectangular relief: Equal Earth, Mollweide (equal-area) and orthographic globes.""" +from __future__ import annotations + +import numpy as np +from PIL import Image, ImageDraw + +from .sketch import sample_equirect +from .sphere import east_north, latlon_to_xyz, xyz_to_latlon + +BACKGROUND = (16, 18, 24) +_A1, _A2, _A3, _A4 = 1.340264, -0.081106, 0.000893, 0.003796 +_M = np.sqrt(3.0) / 2.0 +EE_XMAX = 2.0 * np.sqrt(3.0) * np.pi / (3.0 * _A1) +EE_YMAX = _A1 * np.pi / 3 + _A2 * (np.pi / 3) ** 3 + _A3 * (np.pi / 3) ** 7 + _A4 * (np.pi / 3) ** 9 +MW_XMAX, MW_YMAX = 2.0 * np.sqrt(2.0), np.sqrt(2.0) + + +def equal_earth_forward(lat, lon): + t = np.arcsin(_M * np.sin(np.radians(lat))) + d = _A1 + 3 * _A2 * t**2 + 7 * _A3 * t**6 + 9 * _A4 * t**8 + x = 2 * np.sqrt(3.0) * np.radians(lon) * np.cos(t) / (3 * d) + y = _A1 * t + _A2 * t**3 + _A3 * t**7 + _A4 * t**9 + return x, y + + +def equal_earth_inverse(x, y): + t = np.asarray(y, dtype=np.float64) / _A1 + for _ in range(12): + f = _A1 * t + _A2 * t**3 + _A3 * t**7 + _A4 * t**9 - y + t = t - f / (_A1 + 3 * _A2 * t**2 + 7 * _A3 * t**6 + 9 * _A4 * t**8) + d = _A1 + 3 * _A2 * t**2 + 7 * _A3 * t**6 + 9 * _A4 * t**8 + lat = np.degrees(np.arcsin(np.clip(np.sin(t) / _M, -1, 1))) + lon = np.degrees(3 * np.asarray(x) * d / (2 * np.sqrt(3.0) * np.cos(t))) + return lat, lon + + +def mollweide_forward(lat, lon): + phi = np.radians(lat) + t = phi.copy() if isinstance(phi, np.ndarray) else np.array(phi, dtype=np.float64) + for _ in range(30): + f = 2 * t + np.sin(2 * t) - np.pi * np.sin(phi) + t = t - f / np.maximum(2 + 2 * np.cos(2 * t), 1e-12) + return MW_XMAX / np.pi * np.radians(lon) * np.cos(t), MW_YMAX * np.sin(t) + + +def mollweide_inverse(x, y): + t = np.arcsin(np.clip(np.asarray(y) / MW_YMAX, -1, 1)) + lat = np.degrees(np.arcsin(np.clip((2 * t + np.sin(2 * t)) / np.pi, -1, 1))) + lon = np.degrees(np.pi * np.asarray(x) / (MW_XMAX * np.maximum(np.cos(t), 1e-12))) + return lat, lon + + +def orthographic_inverse(x, y, lat0, lon0): + """Unit-disc view coords → lat/lon on the visible hemisphere centred on (lat0, lon0); ok = inside the disc.""" + x, y = np.asarray(x, dtype=np.float64), np.asarray(y, dtype=np.float64) + rho2 = x**2 + y**2 + ok = rho2 <= 1.0 + z = np.sqrt(np.clip(1.0 - rho2, 0.0, 1.0)) + c = latlon_to_xyz(np.array([lat0]), np.array([lon0])) + e, n = east_north(c) + p = x[..., None] * e[0] + y[..., None] * n[0] + z[..., None] * c[0] + lat, lon = xyz_to_latlon(p) + return lat, lon, ok + + +def _sample_rgb(img, lat, lon): + """Bilinear RGB samples, read straight from the uint8 channels (each sampled value promotes to float64 exactly, + so no full-size float copy of a big image is needed).""" + return np.stack([sample_equirect(img[..., k], lat, lon) for k in range(3)], axis=-1) + + +def reproject(img, proj, width): + """Equirectangular RGB → equal-area world map ('equal_earth' | 'mollweide').""" + xmax, ymax, inv = {"equal_earth": (EE_XMAX, EE_YMAX, equal_earth_inverse), + "mollweide": (MW_XMAX, MW_YMAX, mollweide_inverse)}[proj] + height = int(round(width * ymax / xmax)) + xs = ((np.arange(width) + 0.5) / width * 2 - 1) * xmax + ys = (1 - (np.arange(height) + 0.5) / height * 2) * ymax + X, Y = np.meshgrid(xs, ys) + lat, lon = inv(X, Y) + ok = np.isfinite(lat) & np.isfinite(lon) & (np.abs(lon) <= 180.0) + rgb = _sample_rgb(img, np.where(ok, lat, 0.0), np.where(ok, lon, 0.0)) + out = np.where(ok[..., None], rgb, np.array(BACKGROUND, dtype=np.float64)) + return np.clip(np.round(out), 0, 255).astype(np.uint8) + + +def globe(img, lat0, lon0, size): + """Orthographic view centred on (lat0, lon0), gentle limb darkening.""" + v = (np.arange(size) + 0.5) / size * 2 - 1 + X, Y = np.meshgrid(v, -v) + lat, lon, ok = orthographic_inverse(X, Y, lat0, lon0) + rgb = _sample_rgb(img, lat, lon) * (0.72 + 0.28 * np.sqrt(np.clip(1 - X**2 - Y**2, 0, 1)))[..., None] + out = np.where(ok[..., None], rgb, np.array(BACKGROUND, dtype=np.float64)) + return np.clip(np.round(out), 0, 255).astype(np.uint8) + + +def graticule(arr, proj, step=30): + """Draw a light lat/lon grid (every `step`°) onto an equal-area world map produced by `reproject`.""" + fwd, xmax, ymax = {"equal_earth": (equal_earth_forward, EE_XMAX, EE_YMAX), + "mollweide": (mollweide_forward, MW_XMAX, MW_YMAX)}[proj] + h, w = arr.shape[:2] + im = Image.fromarray(arr) + draw = ImageDraw.Draw(im) + to_px = lambda x, y: list(zip(((x / xmax + 1) / 2 * w).tolist(), ((1 - y / ymax) / 2 * h).tolist())) + t = np.linspace(-90, 90, 181) + for lon in range(-180, 181, step): + draw.line(to_px(*fwd(t, np.full_like(t, float(lon)))), fill=(200, 210, 225), width=1) + s = np.linspace(-180, 180, 361) + for lat in range(-90 + step, 90, step): + draw.line(to_px(*fwd(np.full_like(s, float(lat)), s)), fill=(200, 210, 225), width=1) + return np.array(im) + + +def continent_centres(g, ocean, min_share=0.03): + """Centroids (lat, lon) of land bodies holding ≥ min_share of all land, largest first.""" + from .graph import components + land = ~np.asarray(ocean) + lab = components(g, land) + area = np.bincount(lab[land], weights=g.area_km2[land]) + out = [] + for k in np.argsort(-area): + if area[k] < min_share * area.sum(): + break + m = lab == k + c = (g.xyz[m] * g.area_km2[m, None]).sum(axis=0) + la, lo = xyz_to_latlon(c[None]) + out.append((float(la[0]), float(lo[0]))) + return out + + +def write_all(relief, pdir, g, ocean, width, globe_size=1024): + """Write proj_equal_earth.png, proj_mollweide.png, globe_*.png and globes_sheet.png into pdir.""" + for proj in ("equal_earth", "mollweide"): + Image.fromarray(graticule(reproject(relief, proj, width), proj)).save(pdir / f"proj_{proj}.png") + views = [("north_pole", 90.0, 0.0), ("south_pole", -90.0, 0.0)] + views = [(f"continent_{i + 1}", la, lo) for i, (la, lo) in enumerate(continent_centres(g, ocean))] + views + tiles = [] + for name, la, lo in views: + im = Image.fromarray(globe(relief, la, lo, globe_size)) + im.save(pdir / f"globe_{name}.png") + tiles.append((f"{name} ({la:.0f}°, {lo:.0f}°)", im)) + cols, t = 4, globe_size // 2 + rows = (len(tiles) + cols - 1) // cols + sheet = Image.new("RGB", (cols * t, rows * (t + 16)), BACKGROUND) + draw = ImageDraw.Draw(sheet) + for k, (label, im) in enumerate(tiles): + x, y = (k % cols) * t, (k // cols) * (t + 16) + sheet.paste(im.resize((t, t)), (x, y + 16)) + draw.text((x + 4, y + 2), label, fill=(230, 230, 230)) + sheet.save(pdir / "globes_sheet.png") diff --git a/mapgen/render.py b/mapgen/render.py new file mode 100644 index 0000000..2a408e6 --- /dev/null +++ b/mapgen/render.py @@ -0,0 +1,526 @@ +"""Stage `render`: equirectangular rasters, cells.npz, metadata, previews, contact sheet.""" +from __future__ import annotations + +import colorsys +import json + +import numpy as np +from PIL import Image, ImageDraw +from scipy.spatial import cKDTree + +from . import geo, projections, viewer_export +from . import plateaus as PL +from .crust import AGE_NAMES +from .environment import LEGENDS, REGIONS, ZONES +from .ice import ICE_NAMES +from .noise import fbm +from .seabed import MINERAL_NAMES, SEABED_NAMES +from .sphere import east_north, latlon_to_xyz + +CONTINUOUS = { + "elevation": ("z_surface_m", 0.5, -12000.0, "m"), + "T_mean": ("T_mean", 0.01, -100.0, "degC"), + "T_range": ("T_range", 0.01, 0.0, "degC"), + "P_ann": ("P_ann", 0.5, 0.0, "mm/yr"), + "P_jun": ("P_jun", 0.5, 0.0, "mm/yr"), + "P_dec": ("P_dec", 0.5, 0.0, "mm/yr"), + "po2": ("po2_bar", 2e-4, 0.0, "bar"), + "gravity": ("gravity_g", 1e-4, 0.0, "g"), + "pressure": ("pressure_bar", 2e-4, 0.0, "bar"), + "o2_fraction": ("o2_fraction", 2e-5, 0.0, "fraction"), + "fire": ("fire_reactivity", 1e-4, 0.0, "x"), + "vent_potential": ("vent_potential", 1e-4, 0.0, "0-1"), + "bottom_temp": ("bottom_temp_c", 0.01, -10.0, "degC"), + "sediment": ("sediment_m", 0.5, 0.0, "m"), + "plant_height": ("plant_height_x", 2e-3, 0.0, "x"), + "sst": ("sst", 0.01, -100.0, "degC"), + "productivity": ("productivity", 1e-4, 0.0, "0-1"), + "current_speed": ("current_speed", 1e-4, 0.0, "m/s"), + "upwelling": ("upwelling", 0.05, -1600.0, "m/yr"), +} +CATEGORICAL = {"plates": "plate", "age_class": "age_class", "holdridge": "holdridge", + "seasonality": "seasonality", "landform": "landform", "lithology": "lithology", + "ground": "ground", "ice": "ice", + "seabed_type": "seabed_type", "seabed_mineral": "seabed_mineral", "deposits": "deposit_main"} +DETAIL_M = np.array([60, 40, 150, 400, 80, 150, 300, 300, 300, 60, 30, 120], dtype=np.float64) # by landform +CHUNK = 128 +COAST_DETAIL_M = 250.0 + + +def _by_rows(fn, v, dtype, tail=()): + """fn applied to CHUNK-row slices of v (fn elementwise): the same values as fn(v), without fn's full-size + float64 scratch (a raster ramp builds five H×W×3 float64 temporaries: ~4 GB at 8192 px).""" + v = np.asarray(v) + if v.ndim < 2 or v.shape[0] <= CHUNK: + return fn(v) + out = np.empty(v.shape + tuple(tail), dtype) + for y0 in range(0, v.shape[0], CHUNK): + out[y0:y0 + CHUNK] = fn(v[y0:y0 + CHUNK]) + return out + + +def encode(v, scale, offset): + return _by_rows(lambda x: np.clip(np.round((np.asarray(x, dtype=np.float64) - offset) / scale), 0, 65535) + .astype(np.uint16), v, np.uint16) + + +def decode(raw, scale, offset): + return np.asarray(raw, dtype=np.float64) * scale + offset + + +def _rows_xyz(rows, W, H): + lat = 90.0 - (rows + 0.5) / H * 180.0 + lon = (np.arange(W) + 0.5) / W * 360.0 - 180.0 + LA, LO = np.meshgrid(lat, lon, indexing="ij") + return latlon_to_xyz(LA.ravel(), LO.ravel()) + + +def _make_render_jit(): + """Compiled ramp and sampling loops (numba optional; WORLDGEN_NO_JIT=1 turns it off): per pixel the same float + steps in the same order as the numpy statements they replace, so the same bytes, without the temporaries.""" + import os + if os.environ.get("WORLDGEN_NO_JIT"): + return None + try: + import numba + except ImportError: + return None + + @numba.njit(cache=True, nogil=True) + def ramp(x, lo, span, stops, out): # x: float64 (n,), out: uint8 (n, c) + m = stops.shape[0] + for j in range(x.shape[0]): + t = (x[j] - lo) / span + t = 0.0 if t < 0.0 else (1.0 if t > 1.0 else t) + t = t * (m - 1) + i = min(np.int64(t), m - 2) + f = t - i + for c in range(stops.shape[1]): + out[j, c] = np.uint8(np.int64(stops[i, c] * (1 - f) + stops[i + 1, c] * f)) + + @numba.njit(cache=True, nogil=True) + def sample(v, idx, wts, out): # out[p] = Σ_k v[idx[p,k]] * w[p,k], k left to right (as np.sum, k < 8) + for p in range(idx.shape[0]): + s = v[idx[p, 0]] * np.float64(wts[p, 0]) + for k in range(1, idx.shape[1]): + s += v[idx[p, k]] * np.float64(wts[p, k]) + out[p] = s + return ramp, sample + + +_render_jit = _make_render_jit() + + +def pixel_neighbours(g, W, H, k=3): + from .graph import workers + tree = cKDTree(g.xyz) + idx = np.empty((H, W, k), np.int32) + wts = np.empty((H, W, k), np.float32) + for y0 in range(0, H, CHUNK): + rows = np.arange(y0, min(H, y0 + CHUNK)) + d, i = tree.query(_rows_xyz(rows, W, H), k=k, workers=workers()) # per point: the same answer + w = 1.0 / np.maximum(d, 1e-9) ** 2 + w /= w.sum(axis=1, keepdims=True) + idx[rows] = i.reshape(len(rows), W, k) + wts[rows] = w.reshape(len(rows), W, k) + return idx, wts + + +def sample_cont(v, idx, wts): + out = np.empty(idx.shape[:2]) + if (_render_jit is not None and v.dtype == np.float64 and v.ndim == 1 and wts.dtype == np.float32 + and 1 <= idx.shape[-1] < 8 and idx.shape == wts.shape): + k = idx.shape[-1] + _render_jit[1](np.ascontiguousarray(v), np.ascontiguousarray(idx).reshape(-1, k), + np.ascontiguousarray(wts).reshape(-1, k), out.reshape(-1)) + return out + for y0 in range(0, idx.shape[0], CHUNK): + s = slice(y0, y0 + CHUNK) + out[s] = np.sum(v[idx[s]] * wts[s], axis=-1) + return out + + +def sample_cat(v, idx): + return v[idx[..., 0]] + + +def pixel_land(ocean_k, z): + """Pixel land mask: unanimous neighbour cells decide; at the coast (mixed) the sub-cell height does.""" + all_sea = ocean_k.all(axis=-1) + all_land = ~ocean_k.any(axis=-1) + return all_land | (~all_sea & ~all_land & (z > 0)) + + +def hillshade(z, radius_km, az=315.0, alt=45.0, exag=4.0, low_memory=False): + """low_memory: the same values, CHUNK rows at a time (one halo row each side keeps the central differences).""" + if low_memory: + H = z.shape[0] + out = np.empty(z.shape) + for y0 in range(0, H, CHUNK): + a, b = max(y0 - 1, 0), min(y0 + CHUNK + 1, H) + part = _hillshade_rows(z[a:b], a, H, radius_km, az, alt, exag) + out[y0:min(y0 + CHUNK, H)] = part[y0 - a: y0 - a + min(CHUNK, H - y0)] + return out + return _hillshade_rows(z, 0, z.shape[0], radius_km, az, alt, exag) + + +def _hillshade_rows(z, row0, H, radius_km, az, alt, exag): + """Hillshade of rows row0.. of an H-row raster (z holds those rows; edges of z use one-sided differences).""" + W = z.shape[1] + lat = 90.0 - (np.arange(row0, row0 + z.shape[0]) + 0.5) / H * 180.0 + dy = np.pi * radius_km * 1000.0 / H + dx = 2 * np.pi * radius_km * 1000.0 * np.maximum(np.cos(np.radians(lat)), 0.01) / W + gy, gx = np.gradient(z) + dzdx = gx / dx[:, None] * exag + dzdn = -gy / dy * exag + norm = np.sqrt(dzdx**2 + dzdn**2 + 1.0) + a, b = np.radians(az), np.radians(alt) + L = (np.sin(a) * np.cos(b), np.cos(a) * np.cos(b), np.sin(b)) + return np.clip((-dzdx * L[0] - dzdn * L[1] + L[2]) / norm, 0.0, 1.0) + + +def _ramp(v, lo, hi, stops): + stops = np.asarray(stops, dtype=np.float64) + + def part(x): + if (_render_jit is not None and isinstance(x, np.ndarray) and x.dtype == np.float64 and stops.ndim == 2 + and len(stops) >= 2 and not np.isnan(x).any()): + out = np.empty(x.shape + (stops.shape[1],), np.uint8) + _render_jit[0](np.ascontiguousarray(x).reshape(-1), float(lo), float(hi - lo), stops, + out.reshape(-1, stops.shape[1])) + return out + t = np.clip((x - lo) / (hi - lo), 0, 1) * (len(stops) - 1) + i = np.minimum(t.astype(np.int64), len(stops) - 2) + f = (t - i)[..., None] + return (stops[i] * (1 - f) + stops[i + 1] * f).astype(np.uint8) + return _by_rows(part, v, np.uint8, (stops.shape[-1],)) + + +def holdridge_palette(): + ramp = np.array([[216, 200, 160], [208, 196, 140], [200, 200, 120], [152, 168, 96], [106, 150, 80], + [70, 125, 68], [50, 105, 62], [37, 90, 56]], dtype=np.float64) + tint = {"polar": ((232, 236, 239), 0.9), "subpolar": ((170, 176, 160), 0.55), "boreal": ((70, 100, 80), 0.35), + "tropical": ((20, 90, 40), 0.15)} + cols = [] + for r in REGIONS: + k = len(ZONES[r]) + for j in range(k): + c = ramp[int(round((j / max(k - 1, 1)) * (len(ramp) - 1)))] + if r in tint: + c = c * (1 - tint[r][1]) + np.array(tint[r][0]) * tint[r][1] + cols.append(c) + return np.array(cols, dtype=np.uint8) + + +def category_palette(n): + return np.array([[int(255 * c) for c in colorsys.hsv_to_rgb((i * 0.618034) % 1.0, 0.55, 0.9)] + for i in range(max(n, 1))], dtype=np.uint8) + + +def _region_palette(cols): + """38 zone colours from per-region (dry, wet) colour pairs.""" + out = [] + for r in REGIONS: + dry, wet = (np.array(c, dtype=np.float64) for c in cols[r]) + k = len(ZONES[r]) + out += [dry + (wet - dry) * (j / max(k - 1, 1)) for j in range(k)] + return np.array(out, dtype=np.uint8) + + +STYLES = { # alien palettes (config [render] style); None = Earth-like + "tidal-lock": { + "zones": {"polar": ((34, 38, 46), (34, 38, 46)), "subpolar": ((70, 72, 78), (96, 104, 112)), + "boreal": ((120, 96, 64), (150, 110, 60)), "cool temperate": ((170, 120, 60), (190, 140, 70)), + "warm temperate": ((110, 70, 44), (130, 86, 50)), "subtropical": ((46, 36, 34), (60, 44, 38)), + "tropical": ((22, 20, 22), (30, 26, 28))}, + "ocean": [[4, 12, 16], [10, 34, 40], [40, 80, 84]], "lake": (60, 96, 104), "river": (70, 110, 118), + "ground": {6: (92, 86, 80)}, "sheet": (200, 222, 240), "sea_ice": (170, 200, 226), "sea_ice_seasonal": (120, 150, 170)}, + "salt-mirror": { + "zones": {"polar": ((236, 240, 244), (236, 240, 244)), "subpolar": ((228, 228, 234), (214, 214, 228)), + "boreal": ((238, 236, 228), (208, 208, 224)), "cool temperate": ((242, 238, 226), (200, 204, 222)), + "warm temperate": ((240, 230, 204), (196, 204, 220)), "subtropical": ((238, 224, 186), (192, 206, 218)), + "tropical": ((234, 214, 166), (186, 206, 216))}, + "ocean": [[70, 120, 130], [150, 195, 200], [215, 235, 235]], "lake": (150, 196, 204), "river": (150, 196, 204), + "ground": {6: (252, 250, 244)}, "sheet": (248, 250, 252), "sea_ice": (236, 244, 246), "sea_ice_seasonal": (220, 236, 238)}, +} + + +def style_palette(style): + return None if style is None else STYLES[style] + + +def relief_rgb(z, hs, zone, ground, ice, lake, land=None, vary=None, style=None, low_memory=False): + """vary: optional (brightness factor, tint) per pixel for open land (not lakes or ice): tint > 0 drier/yellower, + < 0 lusher (deeper green). style: a STYLES key (alien palette) or None. low_memory: the same pixels, made + CHUNK rows at a time (no full-size float64 scratch).""" + if low_memory: + out = np.empty(np.shape(z) + (3,), np.uint8) + for y0 in range(0, np.shape(z)[0], CHUNK): + s = slice(y0, y0 + CHUNK) + out[s] = relief_rgb(z[s], hs[s], zone[s], ground[s], ice[s], lake[s], None if land is None else land[s], + None if vary is None else tuple(np.asarray(v)[s] for v in vary), style) + return out + land = z > 0 if land is None else land + st = style_palette(style) + pal = holdridge_palette() if st is None else _region_palette(st["zones"]) + rgb = pal[np.clip(zone, 0, 37)].astype(np.float64) + ocean = _ramp(z, -6500.0 if st is None else -1500.0, 0.0, + [[11, 43, 90], [30, 90, 150], [143, 198, 224]] if st is None else st["ocean"]).astype(np.float64) + rgb = np.where(land[..., None], rgb, ocean) + grounds = ((1, (111, 143, 106)), (2, (120, 128, 100)), (5, (63, 111, 74)), (6, (239, 233, 220))) if st is None \ + else tuple(st["ground"].items()) + for code, col in grounds: + rgb = np.where((land & (ground == code))[..., None], col, rgb) + if vary is not None: + bright, tint = (np.asarray(v, dtype=np.float64)[..., None] for v in vary) + shift = np.where(tint > 0, tint * np.array([1.0, 0.6, -0.8]), -tint * np.array([-0.9, -0.2, -0.6])) + if st is not None: + shift = np.abs(tint) * np.array([0.4, 0.4, 0.4]) * np.sign(tint) + rgb = np.where(land[..., None], rgb * bright + shift, rgb) + rgb = np.where((lake & land)[..., None], (79, 143, 192) if st is None else st["lake"], rgb) + rgb = np.where(np.isin(ice, [1, 2])[..., None], (244, 248, 251) if st is None else st["sheet"], rgb) + rgb = np.where((ice == 4)[..., None], (225, 235, 242) if st is None else st["sea_ice"], rgb) + rgb = np.where((ice == 3)[..., None], 0.5 * rgb + 0.5 * np.array([220, 232, 240] if st is None else st["sea_ice_seasonal"]), rgb) + shade = np.where(land, 0.55 + 0.45 * hs, 0.85 + 0.15 * hs) + return np.clip(rgb * shade[..., None], 0, 255).astype(np.uint8) + + +RIVER_RGB = (58, 112, 176) + + +def draw_rivers(rgb, g, recv, river, strahler, min_order=2, colour=RIVER_RGB): + """Draw river segments (cell centre → receiver) of order ≥ min_order; width grows with order.""" + H, W = rgb.shape[:2] + im = Image.fromarray(rgb) + draw = ImageDraw.Draw(im) + x = (np.asarray(g.lon) + 180.0) / 360.0 * W - 0.5 + y = (90.0 - np.asarray(g.lat)) / 180.0 * H - 0.5 + cells = np.flatnonzero(river & (recv != np.arange(g.n)) & (strahler >= min_order)) + for i in cells[np.argsort(strahler[cells])]: + j = recv[i] + if abs(x[i] - x[j]) > W / 2: + continue + width = max(1, int(round(int(strahler[i]) * W / 8192))) + draw.line([(x[i], y[i]), (x[j], y[j])], fill=tuple(colour), width=width) + return np.array(im) + + +def draw_currents(rgb, g, current, ocean, per_row=60, colour=(255, 255, 255)): + """Arrows along the surface current on a lattice ≈ `per_row` across the map; length ∝ speed (1 m/s ≈ one + lattice step), sea only; currents under 2 cm/s get none.""" + H, W = rgb.shape[:2] + step = W / per_row + current = np.asarray(current, dtype=np.float64) + ocean = np.asarray(ocean, bool) + e, n = east_north(g.xyz) + tree = cKDTree(g.xyz) + ys, xs = np.meshgrid(np.arange(step / 2, H, step), np.arange(step / 2, W, step), indexing="ij") + lat, lon = 90.0 - (ys.ravel() + 0.5) / H * 180.0, (xs.ravel() + 0.5) / W * 360.0 - 180.0 + _, cell = tree.query(latlon_to_xyz(lat, lon)) + ue, un = np.sum(current[cell] * e[cell], axis=1), np.sum(current[cell] * n[cell], axis=1) + sp = np.hypot(ue, un) + im = Image.fromarray(rgb) + draw = ImageDraw.Draw(im) + for x, y, a, b, s, ok in zip(xs.ravel(), ys.ravel(), ue, un, sp, ocean[cell]): + if not ok or s < 0.02: + continue + L = min(s, 1.5) * 0.9 * step + dx, dy = a / s * L, -b / s * L + x0, y0, x1, y1 = x - dx / 2, y - dy / 2, x + dx / 2, y + dy / 2 + draw.line([(x0, y0), (x1, y1)], fill=tuple(colour), width=1) + for ang in (2.6, -2.6): # arrowhead: two barbs ±150° + c, sn = np.cos(ang), np.sin(ang) + draw.line([(x1, y1), (x1 + 0.35 * (c * dx - sn * dy), y1 + 0.35 * (sn * dx + c * dy))], + fill=tuple(colour), width=1) + return np.array(im) + + +RAMPS = { # key: (unit, lo, hi, colour stops) — shared by previews, viewer textures and legends + "elevation": ("m", -6000.0, 6000.0, [[8, 30, 70], [30, 90, 150], [140, 200, 225], [60, 120, 60], + [150, 160, 90], [140, 110, 80], [235, 235, 235]]), + "T_mean": ("degC", -40.0, 40.0, [[40, 60, 160], [240, 240, 240], [180, 30, 30]]), + "T_range": ("degC", 0.0, 60.0, [[68, 1, 84], [33, 145, 140], [253, 231, 37]]), + "P": ("mm/yr", 0.0, 4000.0, [[150, 110, 60], [230, 220, 150], [60, 150, 70], [30, 70, 170]]), + "po2": ("bar", 0.08, 0.32, [[60, 40, 90], [240, 240, 240], [200, 90, 20]]), + "gravity": ("g", 0.3, 1.8, [[20, 120, 180], [240, 240, 240], [120, 40, 40]]), + "pressure": ("bar", 0.3, 2.5, [[40, 60, 120], [240, 240, 240], [150, 60, 30]]), + "o2_fraction": ("fraction", 0.10, 0.40, [[60, 40, 90], [240, 240, 240], [200, 90, 20]]), + "fire": ("x", 0.0, 1.5, [[40, 80, 160], [240, 240, 240], [220, 120, 20]]), + "vent_potential": ("0-1", 0.0, 1.0, [[10, 20, 50], [60, 90, 160], [250, 190, 60], [255, 250, 220]]), + "bottom_temp": ("degC", -2.0, 30.0, [[30, 40, 120], [80, 160, 200], [240, 200, 120]]), + "sediment": ("m", 0.0, 3000.0, [[40, 30, 60], [140, 110, 80], [240, 220, 170]]), + "plant_height": ("x", 0.5, 3.5, [[150, 120, 60], [240, 240, 240], [30, 110, 50]]), + "sst": ("degC", -2.0, 32.0, [[30, 40, 120], [60, 150, 200], [240, 240, 200], [220, 90, 40]]), + "productivity": ("0-1", 0.0, 1.0, [[10, 20, 60], [20, 110, 120], [120, 200, 90], [240, 240, 120]]), + "current_speed": ("m/s", 0.0, 1.0, [[10, 20, 50], [40, 90, 170], [120, 200, 230], [255, 255, 255]]), + "upwelling": ("m/yr", -200.0, 200.0, [[40, 60, 160], [240, 240, 240], [30, 140, 80]]), +} +VIEWER_LAYERS = [ # id, name, source (continuous raster name | categorical name | "relief") + ("relief", "Relief", "relief"), ("biomes", "Biomes (Holdridge)", "holdridge"), + ("elevation", "Elevation", "elevation"), ("temperature", "Mean temperature", "T_mean"), + ("rainfall", "Rainfall", "P_ann"), ("seasonality", "Seasonality", "seasonality"), + ("landform", "Landform", "landform"), ("ground", "Ground", "ground"), ("ice", "Ice", "ice"), ("deposits", "Mineral deposits", "deposits"), + ("plates", "Plates", "plates"), ("o2", "O₂ partial pressure", "po2"), ("gravity", "Gravity", "gravity"), + ("pressure", "Air pressure", "pressure"), ("fire", "Fire reactivity", "fire"), + ("seabed", "Sea-floor type", "seabed_type"), ("minerals", "Sea-floor minerals", "seabed_mineral"), + ("bottom_temp", "Bottom temperature", "bottom_temp"), ("sediment", "Sediment", "sediment"), + ("vent_potential", "Vent potential", "vent_potential"), + ("currents", "Ocean currents", "current_speed"), ("sst", "Sea-surface temperature", "sst"), + ("productivity", "Sea productivity", "productivity"), +] + + +def _ramp_key(name): + return "P" if name.startswith("P_") else name + + +def _colorize(name, v): + _, lo, hi, stops = RAMPS[_ramp_key(name)] + return _ramp(v, lo, hi, stops) + + +def _preview(img: Image.Image, width: int, nearest: bool) -> Image.Image: + return img.resize((width, width // 2), Image.NEAREST if nearest else Image.LANCZOS) + + +def contact_sheet(tiles, tile_w): + cols = 3 + th = tile_w // 2 + 14 + rows = (len(tiles) + cols - 1) // cols + sheet = Image.new("RGB", (cols * tile_w, rows * th), (20, 20, 24)) + draw = ImageDraw.Draw(sheet) + for k, (name, im) in enumerate(tiles): + x, y = (k % cols) * tile_w, (k // cols) * th + sheet.paste(im.resize((tile_w, tile_w // 2)), (x, y + 14)) + draw.text((x + 4, y + 1), name, fill=(235, 235, 235)) + return sheet + + +class _Writer: + """Saves files on a background thread (PNG/zlib encoding releases the GIL) while the next layer is computed; + at most `depth` waiting, so memory stays bounded. The files are the same bytes as saved in line.""" + def __init__(self, threads: int = 2, depth: int = 4): + from concurrent.futures import ThreadPoolExecutor + self.ex, self.pending, self.depth = ThreadPoolExecutor(threads), [], depth + + def __call__(self, fn, *args, **kw): + while len(self.pending) >= self.depth: + self.pending.pop(0).result() + self.pending.append(self.ex.submit(fn, *args, **kw)) + + def close(self): + try: + for f in self.pending: + f.result() + finally: + self.ex.shutdown(wait=True, cancel_futures=True) + + +def run(ctx) -> dict: + writer = _Writer() + try: + return _run(ctx, writer) + finally: + writer.close() + + +def _run(ctx, save) -> dict: + g, cfg, d = ctx.grid, ctx.cfg, ctx.data + W = int(cfg["build"]["raster_width"]) + if ctx.res < int(cfg["build"]["res_final"]): + W = int(cfg["build"].get("dev_raster_width", W)) + H = W // 2 + pw = int(cfg["build"]["preview_width"]) + out = ctx.out_dir or ctx.root / "out" / f"r{ctx.res}" + rdir = out / "raster" + pdir = ctx.preview_dir or ctx.root / "previews" / f"r{ctx.res}" + rdir.mkdir(parents=True, exist_ok=True) + pdir.mkdir(parents=True, exist_ok=True) + idx, wts = pixel_neighbours(g, W, H) + legends = {**LEGENDS, "age_class": AGE_NAMES, "ice": ICE_NAMES, "plates": [p["id"] for p in ctx.tect["plate"]], + "seabed_type": SEABED_NAMES, "seabed_mineral": MINERAL_NAMES} + meta = {"width": W, "height": H, "projection": "equirectangular", "radius_km": g.radius_km, + "units": cfg.get("units", {}), "continuous": {}, "categorical": {}} + tiles = [] + cats = {name: sample_cat(np.asarray(d[key]), idx) for name, key in CATEGORICAL.items()} + + z = sample_cont(np.asarray(d["z_surface_m"], dtype=np.float64), idx, wts) + amp = DETAIL_M[np.clip(cats["landform"], 0, len(DETAIL_M) - 1)] + for y0 in range(0, H, CHUNK): + rows = np.arange(y0, min(H, y0 + CHUNK)) + z[rows] += amp[rows] * fbm(_rows_xyz(rows, W, H), ctx.seed + 61, 5, 64.0).reshape(len(rows), W) + lake = sample_cat(np.asarray(d["lake"]), idx) + ocean_cells = np.asarray(d["ocean"]) if "ocean" in d else np.asarray(d["z_surface_m"]) <= 0 + ocean_k = ocean_cells[idx] + coast = ocean_k.any(axis=-1) & ~ocean_k.all(axis=-1) + for y0 in range(0, H, CHUNK): # extra fine detail where the coastline runs + rows = np.arange(y0, min(H, y0 + CHUNK)) + cz = COAST_DETAIL_M * fbm(_rows_xyz(rows, W, H), ctx.seed + 67, 5, 256.0).reshape(len(rows), W) + z[rows] += np.where(coast[rows], cz, 0.0) + land = pixel_land(ocean_k, z) + if "lake_level_m" in d: # the water surface for drawing: lake cells at their level (beds stay in `elevation`) + lev = np.asarray(d["lake_level_m"], dtype=np.float64) + dz = sample_cont(np.where(np.isfinite(lev), lev - np.asarray(d["z_surface_m"], dtype=np.float64), 0.0), idx, wts) + save(Image.fromarray(encode(z + dz, *CONTINUOUS["elevation"][1:3])).save, rdir / "surface.png") + meta["continuous"]["surface"] = {"file": "surface.png", "scale": CONTINUOUS["elevation"][1], + "offset": CONTINUOUS["elevation"][2], "unit": "m"} + hs = hillshade(z, g.radius_km, low_memory=ctx.low_memory) + style = cfg.get("render", {}).get("style") + rgb = relief_rgb(z, hs, cats["holdridge"], cats["ground"], cats["ice"], lake, land, style=style, + low_memory=ctx.low_memory) + rgb = draw_rivers(rgb, g, np.asarray(d["recv"]), np.asarray(d["river"]), np.asarray(d["strahler"]), + colour=RIVER_RGB if style is None else STYLES[style]["river"]) + relief = Image.fromarray(rgb) + projections.write_all(rgb, pdir, g, ocean_cells, pw, globe_size=max(pw // 2, 64)) + save(relief.copy().save, rdir / "relief.png") + tiles.append(("relief", _preview(relief, pw, False))) + + vdir = out / "viewer" + by_src = {} # viewer layers are written as soon as their colours exist + for vid, vname, src in VIEWER_LAYERS: + by_src.setdefault(src, []).append((vid, vname)) + ventries = {} + + def emit(src, rgb_, legend): + for vid, vname in by_src.get(src, []): + ventries[vid] = viewer_export.write_layer(vdir, {"id": vid, "name": vname, "rgb": rgb_, "legend": legend}) + + emit("relief", rgb, None) + for name, (key, scale, offset, unit) in CONTINUOUS.items(): + v = z if name == "elevation" else sample_cont(np.asarray(d[key], dtype=np.float64), idx, wts) + save(Image.fromarray(encode(v, scale, offset)).save, rdir / f"{name}.png") + meta["continuous"][name] = {"file": f"{name}.png", "scale": scale, "offset": offset, "unit": unit} + crgb = _colorize(name, v) + if name == "current_speed" and "current" in d: + crgb = draw_currents(crgb, g, d["current"], ocean_cells) + if name != "elevation": + tiles.append((name, _preview(Image.fromarray(crgb), pw, False))) + if name in by_src: + unit_, lo, hi, stops = RAMPS[_ramp_key(name)] + emit(name, crgb, viewer_export.continuous_legend(unit_, lo, hi, stops)) + del crgb + for name, v in cats.items(): + save(Image.fromarray(v.astype(np.uint8), "L").save, rdir / f"{name}.png") + meta["categorical"][name] = {"file": f"{name}.png", "legend": legends[name]} + pal = holdridge_palette() if name == "holdridge" else category_palette(len(legends[name])) + crgb = pal[np.clip(v, 0, len(pal) - 1)] + tiles.append((name, _preview(Image.fromarray(crgb), pw, True))) + emit(name, crgb, viewer_export.categorical_legend(legends[name], pal)) + del crgb + (out / "fields.json").write_text(json.dumps(meta, indent=1)) + + per_cell = {k: v for k, v in d.items() if isinstance(v, np.ndarray) and v.shape[:1] == (g.n,)} + save(lambda: np.savez_compressed(out / "cells.npz", **per_cell)) + (out / "cells_meta.json").write_text(json.dumps( + {"res": ctx.res, "radius_km": g.radius_km, "n_cells": g.n, "units": cfg.get("units", {}), + "planet": {**cfg["planet"], **({"sun_lock": cfg["climate"].get("lock_at", [0.0, 0.0])} + if cfg.get("climate", {}).get("lock") else {})}, "name": cfg.get("render", {}).get("name", "World"), "style": style, + "legends": legends, "fields": sorted(per_cell), + "plateaus": [PL.features(p, ctx.seed, g.radius_km) for p in ctx.tect.get("plateau", [])]}, indent=1)) + + for name, im in tiles: + im.save(pdir / f"{name}.png") + contact_sheet(tiles, max(pw // 3, 64)).save(pdir / "contact_sheet.png") + viewer_export.write_index(vdir, [ventries[vid] for vid, _, _ in VIEWER_LAYERS]) + geo.write_all(ctx, out, land, lake) + return {} diff --git a/mapgen/seabed.py b/mapgen/seabed.py new file mode 100644 index 0000000..7715c44 --- /dev/null +++ b/mapgen/seabed.py @@ -0,0 +1,95 @@ +"""Stage `seabed`: sea-floor data for every ocean cell — vent potential, +sea-floor type, minerals, bottom temperature, sediment. Land cells get 0 / "—" (the fields describe the sea floor).""" +from __future__ import annotations + +import numpy as np +from scipy.spatial import cKDTree + +from . import plateaus as PL +from .config import params +from .elevation import hotspot_track +from .graph import distance_to +from .noise import fbm +from .sphere import latlon_to_xyz + +(SB_NONE, SB_TERRIGENOUS, SB_CARBONATE, SB_CLAY, SB_BASALT, SB_VOLCANIC, SB_CONTINENTAL, SB_VENTS, + SB_TRENCH) = range(9) +SEABED_NAMES = ["—", "terrigenous sediment", "carbonate ooze", "pelagic clay", "young basalt", "volcanic", + "continental / plateau", "vent field", "trench"] +MI_NONE, MI_SULFIDES, MI_NODULES, MI_COBALT, MI_PHOSPHORITE = range(5) +MINERAL_NAMES = ["none", "polymetallic sulfides", "manganese nodules", "cobalt crusts", "phosphorite"] +DEFAULTS = {"ridge_km": 100.0, "back_arc_km": 300.0, "back_arc_width_km": 150.0, "back_arc": 0.6, + "hotspot_km": 150.0, "hotspot_spacing_km": 150.0, "terrigenous_km": 300.0, "carbonate_depth_m": 4500.0, + "carbonate_t_c": 10.0, "young_myr": 10.0, "trench_km": 80.0, "trench_depth_m": 5000.0, + "vent_field": 0.8, "sulfides": 0.7, "volcanic": 0.3, "thermocline_m": 800.0} + + +def _cells_near(tree, radius_km, p, reach_km): + return np.asarray(tree.query_ball_point(p, 2.0 * np.sin(min(np.pi, reach_km / radius_km) / 2.0)), dtype=np.int64) + + +def _bump(v, g, tree, p, r_km, amp): + i = _cells_near(tree, g.radius_km, p, 3.0 * r_km) + if len(i): + d = g.radius_km * np.arccos(np.clip(g.xyz[i] @ p, -1.0, 1.0)) + v[i] = np.maximum(v[i], amp * np.exp(-(d / r_km) ** 2)) + + +def volcanic_potential(g, tect: dict, vel, seed: int, P: dict): + """0–1 per cell from hotspot chains (strongest at the active end) and plateau volcanic fields.""" + tree, v = cKDTree(g.xyz), np.zeros(g.n) + for h in tect.get("hotspot", []): + for p, s in hotspot_track(g, h, vel, P["hotspot_spacing_km"]): + _bump(v, g, tree, p, P["hotspot_km"], 1.0 - s / h["length_km"]) + for pl in tect.get("plateau", []): + f, vent = PL.features(pl, seed, g.radius_km), float(pl.get("vent", 1.0)) + _bump(v, g, tree, latlon_to_xyz(*f["hotspot"]), 200.0, vent) + for c in f["cones"]: + _bump(v, g, tree, latlon_to_xyz(c["lat"], c["lon"]), c["radius_km"] + 30.0, 0.8 * vent) + for c in f["calderas"]: + _bump(v, g, tree, latlon_to_xyz(c["lat"], c["lon"]), c["radius_km"] + 30.0, vent) + return v + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "seabed", DEFAULTS) + z, ocean, cont, age, d_div, d_over, d_sub, t_mean, vel = ctx.need( + "elevation_eroded_m", "ocean", "continental", "ocean_age_myr", "d_div_km", "d_over_km", "d_sub_km", "T_mean", + "vel") + z, ocean, cont = np.asarray(z, dtype=np.float64), np.asarray(ocean, bool), np.asarray(cont, bool) + age, t_mean = np.asarray(age, dtype=np.float64), np.asarray(t_mean, dtype=np.float64) + plateau_id = np.asarray(ctx.data.get("plateau_id", np.full(g.n, -1))) + depth = np.maximum(-z, 0.0) + volc = volcanic_potential(g, ctx.tect, np.asarray(vel), ctx.seed, P) + ridge = np.exp(-(np.asarray(d_div) / P["ridge_km"]) ** 2) + back_arc = P["back_arc"] * np.exp(-((np.asarray(d_over) - P["back_arc_km"]) / P["back_arc_width_km"]) ** 2) + cluster = np.clip(0.5 + 1.5 * fbm(g.xyz, ctx.seed + 301, 4, 30.0), 0.0, 1.0) # vents come in fields + vent = np.where(ocean, np.maximum.reduce([ridge, back_arc, volc]) * (0.5 + 0.5 * cluster), 0.0) + d_land = distance_to(g, ~ocean) if (~ocean).any() else np.full(g.n, np.inf) + + t = np.full(g.n, SB_CLAY, np.int8) + t[(depth < P["carbonate_depth_m"]) & (t_mean > P["carbonate_t_c"])] = SB_CARBONATE + t[cont] = SB_CONTINENTAL + t[d_land < P["terrigenous_km"]] = SB_TERRIGENOUS + t[~cont & (age < P["young_myr"])] = SB_BASALT + t[volc > P["volcanic"]] = SB_VOLCANIC + t[(np.asarray(d_sub) < P["trench_km"]) & (depth > P["trench_depth_m"])] = SB_TRENCH + t[vent > P["vent_field"]] = SB_VENTS + t[~ocean] = SB_NONE + + m = np.zeros(g.n, np.int8) + m[(t == SB_CLAY) & (age > 30.0) & (d_land > 1000.0) & (depth > 4000.0)] = MI_NODULES + m[((t == SB_VOLCANIC) | (plateau_id >= 0)) & (depth >= 800.0) & (depth <= 2500.0)] = MI_COBALT + m[cont & (depth < 500.0)] = MI_PHOSPHORITE + m[vent > P["sulfides"]] = MI_SULFIDES + m[~ocean] = MI_NONE + + deep = 1.0 + 3.0 * np.cos(np.radians(g.lat)) ** 2 # ≈ 1 °C polar … 4 °C tropical + w = np.clip((P["thermocline_m"] - depth) / P["thermocline_m"], 0.0, 1.0) # shallow floors: towards the surface + bottom = np.where(ocean, deep + w * (np.maximum(t_mean, -1.8) - deep), 0.0) + terr = 2500.0 * np.exp(-d_land / 250.0) # aprons off the continents + sed = np.where(cont, 300.0 + terr, np.minimum(5.0 * age, 800.0) + terr) # pelagic rain ≈ 5 m per Myr + sed = np.where(ocean, sed, 0.0) + return {"vent_potential": vent.astype(np.float32), "seabed_type": t, "seabed_mineral": m, + "bottom_temp_c": bottom.astype(np.float32), "sediment_m": sed.astype(np.float32)} diff --git a/mapgen/sketch.py b/mapgen/sketch.py new file mode 100644 index 0000000..a0249bf --- /dev/null +++ b/mapgen/sketch.py @@ -0,0 +1,121 @@ +"""Stage `sketch`: sample the sketch + user masks onto cells.""" +from __future__ import annotations + +from pathlib import Path + +import numpy as np +from PIL import Image + +from . import zones as ZN +from .config import params +from .noise import fbm +from .pipeline import StageError +from .sphere import latlon_to_xyz, xyz_to_latlon + +W, H = 2000, 1000 +SKETCH = ("land", "mountains", "desert", "rainforest", "trench") +MASKS = ("land_hint", "mountain_hint", "o2_zones", "gravity_zones", "lock") +DEFAULTS = {"warp_km": 1000.0, "warp_freq": 1.2, "detail_warp_km": 200.0, "detail_warp_freq": 5.0, "moves": []} + + +def load_png01(path: Path) -> np.ndarray: + return np.array(Image.open(path).convert("L"), dtype=np.float64) / 255.0 + + +def load_mask(path: Path) -> np.ndarray: + """Grayscale override mask → [−1, 1]; mid-grey (128 / 32768), transparency and missing detail = neutral 0.""" + try: + im = Image.open(path) + im.load() + except (OSError, ValueError) as e: + raise StageError(f"cannot read mask {path.name}: {e}") from e + if im.mode in ("I", "I;16", "I;16B", "I;16L"): + v = (np.array(im).astype(np.float64) - 32768.0) / 32768.0 + dead = 1.0 / 32768.0 + else: + la = np.array(im.convert("LA")).astype(np.float64) + v = np.where(la[..., 1] > 0, (la[..., 0] - 127.5) / 127.5, 0.0) + dead = 1.0 / 255.0 + return np.clip(np.where(np.abs(v) <= dead, 0.0, v), -1.0, 1.0) + + +def sample_equirect(img: np.ndarray, lat, lon) -> np.ndarray: + h, w = img.shape + x = (np.asarray(lon, dtype=np.float64) + 180.0) / 360.0 * w - 0.5 + y = (90.0 - np.asarray(lat, dtype=np.float64)) / 180.0 * h - 0.5 + x0 = np.floor(x).astype(np.int64) + y0 = np.floor(y).astype(np.int64) + fx, fy = x - x0, y - y0 + xa, xb = x0 % w, (x0 + 1) % w + ya, yb = np.clip(y0, 0, h - 1), np.clip(y0 + 1, 0, h - 1) + top = img[ya, xa] * (1 - fx) + img[ya, xb] * fx + bot = img[yb, xa] * (1 - fx) + img[yb, xb] * fx + return top * (1 - fy) + bot * fy + + +def rotation(a, b): + """3×3 rotation carrying [lat, lon] a onto b along the great circle (Rodrigues).""" + pa, pb = latlon_to_xyz(*np.array(a, dtype=np.float64)), latlon_to_xyz(*np.array(b, dtype=np.float64)) + k = np.cross(pa, pb) + s, c = np.linalg.norm(k), float(np.dot(pa, pb)) + if s < 1e-12: + return np.eye(3) + k /= s + K = np.array([[0, -k[2], k[1]], [k[2], 0, -k[0]], [-k[1], k[0], 0]]) + return np.eye(3) + s * K + (1 - c) * K @ K + + +def continent_mask(land, lat, lon): + """Pixels of the sketch landmass containing (lat, lon), 2 px coastal fringe included; lon wraps.""" + from scipy import ndimage + lab, _ = ndimage.label(land) + for a, b in zip(lab[:, 0], lab[:, -1]): + if a and b and a != b: + lab[lab == b] = a + h, w = land.shape + y, x = min(max(int((90.0 - lat) / 180.0 * h), 0), h - 1), int((lon + 180.0) / 360.0 * w) % w + if not lab[y, x]: + raise StageError(f"sketch move: no sketch land at [{lat}, {lon}]") + return ndimage.binary_dilation(lab == lab[y, x], iterations=2) + + +def warped_latlon(g, P, seed): + """Sample positions displaced by a two-scale domain warp (RMS ≈ warp_km, detail_warp_km).""" + p = g.xyz.copy() + for amp, freq, s in ((P["warp_km"], P["warp_freq"], 201), (P["detail_warp_km"], P["detail_warp_freq"], 211)): + if amp <= 0: + continue + disp = np.stack([fbm(g.xyz, seed + s + k, 4, freq) for k in range(3)], axis=1) + p = p + (amp / g.radius_km) * disp / max(float(disp.std()), 1e-12) + return xyz_to_latlon(p) + + +def run(ctx) -> dict: + g = ctx.grid + P = params(ctx.cfg, "sketch", DEFAULTS) + sk = ctx.root / "sketch" + missing = [n for n in SKETCH if not (sk / f"{n}.png").exists()] + if missing: + raise StageError(f"sketch files missing {missing}: draw them or run `mapgen.py new-world` first") + lat, lon = warped_latlon(g, P, ctx.seed) # the drawing is a loose guide: warp it (user masks are not) + imgs = {n: load_png01(sk / f"{n}.png") for n in SKETCH} + moved = [] # [[sketch.moves]]: rotate whole landmasses across the sphere + for mv in P["moves"]: + comp = continent_mask(imgs["land"] > 0.5, *mv["at"]) + moved.append((rotation(mv["at"], mv["to"]), {n: np.where(comp, im, 0.0) for n, im in imgs.items()})) + imgs = {n: np.where(comp, 0.0, im) for n, im in imgs.items()} + p = latlon_to_xyz(lat, lon) + out = {} + for n in SKETCH: + v = sample_equirect(imgs[n], lat, lon) + for rot, parts in moved: + qlat, qlon = xyz_to_latlon(p @ rot) # rot.T applied to row vectors + v = np.maximum(v, sample_equirect(parts[n], qlat, qlon)) + out[f"sk_{n}"] = v + for m in MASKS: + p = ctx.root / "masks" / f"{m}.png" + out[f"m_{m}"] = sample_equirect(load_mask(p), g.lat, g.lon) if p.exists() else np.zeros(g.n) + zo2, zg = ZN.masks(g.xyz, ctx.tect.get("zone", []), ctx.seed, g.radius_km) # config zones on top of the masks + out["m_o2_zones"] = np.clip(out["m_o2_zones"] + zo2, -1.0, 1.0) + out["m_gravity_zones"] = np.clip(out["m_gravity_zones"] + zg, -1.0, 1.0) + return out diff --git a/mapgen/sphere.py b/mapgen/sphere.py new file mode 100644 index 0000000..394b9ef --- /dev/null +++ b/mapgen/sphere.py @@ -0,0 +1,68 @@ +"""Unit-sphere geometry and plate kinematics. Vectors are (..., 3); lat/lon in degrees.""" +from __future__ import annotations + +import numpy as np + + +def latlon_to_xyz(lat, lon): + la, lo = np.radians(lat), np.radians(lon) + return np.stack([np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)], axis=-1) + + +def xyz_to_latlon(p): + p = p / np.linalg.norm(p, axis=-1, keepdims=True) + return np.degrees(np.arcsin(np.clip(p[..., 2], -1, 1))), np.degrees(np.arctan2(p[..., 1], p[..., 0])) + + +def east_north(p): + """Local unit east/north vectors. At the poles east is taken as +y (finite, arbitrary).""" + p = np.asarray(p, dtype=np.float64) + e = np.stack([-p[..., 1], p[..., 0], np.zeros_like(p[..., 0])], axis=-1) # z × p + ne = np.linalg.norm(e, axis=-1, keepdims=True) + polar = ne[..., 0] < 1e-12 + e = np.where(polar[..., None], np.array([0.0, 1.0, 0.0]), e / np.maximum(ne, 1e-300)) + n = np.cross(p, e) + return e, n + + +def tangent_dir(p, q): + """Unit tangent at p pointing along the great circle toward q.""" + t = q - np.sum(p * q, axis=-1, keepdims=True) * p + return t / np.maximum(np.linalg.norm(t, axis=-1, keepdims=True), 1e-15) + + +def gc_dist_km(p, q, radius_km): + return radius_km * np.arccos(np.clip(np.sum(p * q, axis=-1), -1.0, 1.0)) + + +def motion_to_omega(lat, lon, azimuth_deg, speed_cm_yr, radius_km): + """Angular velocity (rad/yr) that moves the point (lat, lon) along azimuth at speed.""" + p = latlon_to_xyz(np.array([lat]), np.array([lon])) + e, n = east_north(p) + az = np.radians(azimuth_deg) + v = (np.sin(az) * e[0] + np.cos(az) * n[0]) * (speed_cm_yr / 100.0) # m/yr + return np.cross(p[0], v) / (radius_km * 1000.0) + + +def velocity(p, omega, radius_km): + """Surface velocity (m/yr) of points p on plates with angular velocity omega (broadcast).""" + return np.cross(omega, p) * (radius_km * 1000.0) + + +def rotate_about(p, v, ang): + """Rotate tangent vectors v about normals p by ang radians (counter-clockwise seen from outside).""" + ang = np.asarray(ang, dtype=np.float64)[..., None] + return v * np.cos(ang) + np.cross(p, v) * np.sin(ang) + + +def azimuth_deg(center, p): + """Bearing (deg clockwise from north, [0, 360)) from center (3,) to points p (..., 3).""" + e, n = east_north(center[None]) + t = tangent_dir(np.broadcast_to(center, p.shape), p) + return np.degrees(np.arctan2(t @ e[0], t @ n[0])) % 360.0 + + +def great_circle_point(p0, t0, dist_km, radius_km): + """Point reached from p0 moving dist_km along unit tangent t0.""" + a = dist_km / radius_km + return np.cos(a) * p0 + np.sin(a) * t0 diff --git a/mapgen/testdata/world.toml b/mapgen/testdata/world.toml new file mode 100644 index 0000000..e650e39 --- /dev/null +++ b/mapgen/testdata/world.toml @@ -0,0 +1,31 @@ +# Test fixture world (the unit tests' planet); not an example — see example/config/world.toml. +[planet] +radius_km = 12742.0 +gravity_g = 1.05 +day_hours = 31.149 +year_days = 216 +tilt_deg = 20.0 +sea_level_pressure_bar = 1.0 +scale_height_km = 8.0 +o2_fraction = 0.21 + +[units] +span_m = 0.331013 +moment_s = 2.403473 +league_km = 15.4437 + +[build] +seed = 1296 +res_dev = 4 +res_final = 5 +raster_width = 8192 +dev_raster_width = 2048 +preview_width = 2048 +land_fraction = 0.19 + +[climate] +hadley_edge_deg = 20.0 +ferrel_edge_deg = 55.0 + +[sketch] +moves = [] diff --git a/mapgen/testing.py b/mapgen/testing.py new file mode 100644 index 0000000..b70e5fe --- /dev/null +++ b/mapgen/testing.py @@ -0,0 +1,59 @@ +"""Test fixtures shared with tools built on mapgen (e.g. a viewer's tests): a tiny two-continent world at H3 res 2.""" +from __future__ import annotations + +import atexit +import functools +import re +import shutil +import tempfile +from pathlib import Path + +import numpy as np +from PIL import Image + +FIXTURE_TOML = Path(__file__).resolve().parent / "testdata" / "world.toml" + +TECT = """ +[[plate]] +id = "a" +seed = [10.0, -30.0] +kind = "continental" +motion = [90.0, 4.0] +[[plate]] +id = "b" +seed = [10.0, 30.0] +kind = "continental" +motion = [270.0, 4.0] +[[plate]] +id = "c" +seed = [0.0, 150.0] +kind = "oceanic" +motion = [0.0, 3.0] +""" + + +def small_world(tmp: Path) -> None: + (tmp / "config").mkdir() + w = FIXTURE_TOML.read_text() + w = w.replace("raster_width = 8192", "raster_width = 256").replace("preview_width = 2048", "preview_width = 128") + w = w.replace("dev_raster_width = 2048", "dev_raster_width = 256") + w = re.sub(r"moves = \[.*?\n\]", "moves = []", w, flags=re.S) + (tmp / "config" / "world.toml").write_text(w) + (tmp / "config" / "tectonics.toml").write_text(TECT) + (tmp / "sketch").mkdir() + lat = 90 - (np.arange(100) + 0.5) * 1.8 + lon = (np.arange(200) + 0.5) * 1.8 - 180 + LA, LO = np.meshgrid(lat, lon, indexing="ij") + land = ((np.abs(LA - 10) < 35) & (np.abs(LO) < 70)).astype(np.uint8) * 255 + for n in ("land", "mountains", "desert", "rainforest", "trench"): + Image.fromarray(land if n == "land" else np.zeros_like(land), "L").save(tmp / "sketch" / f"{n}.png") + + +@functools.lru_cache(maxsize=1) +def built_world() -> Path: + from mapgen import pipeline as P + tmp = Path(tempfile.mkdtemp(prefix="worldgen-fixture-")) + atexit.register(shutil.rmtree, tmp, True) # 5 MB per build: never leave them in /tmp + small_world(tmp) + P.build(tmp, 2, log=lambda m: None) + return tmp diff --git a/mapgen/viewer_export.py b/mapgen/viewer_export.py new file mode 100644 index 0000000..e6f51b4 --- /dev/null +++ b/mapgen/viewer_export.py @@ -0,0 +1,39 @@ +"""Viewer textures: colourised equirectangular JPEG layers + layers.json.""" +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np +from PIL import Image + + +def hex_color(c) -> str: + return "#%02x%02x%02x" % tuple(int(v) for v in list(c)[:3]) + + +def categorical_legend(names, palette) -> dict: + return {"type": "categorical", + "items": [{"name": n, "color": hex_color(palette[i % len(palette)])} for i, n in enumerate(names)]} + + +def continuous_legend(unit, lo, hi, stops) -> dict: + vals = np.linspace(lo, hi, len(stops)) + return {"type": "continuous", "unit": unit, + "stops": [[round(float(v), 4), hex_color(s)] for v, s in zip(vals, stops)]} + + +def write_layer(vdir: Path, L: dict) -> dict: + """Write one layer's JPEG; returns its layers.json entry.""" + vdir.mkdir(parents=True, exist_ok=True) + Image.fromarray(np.ascontiguousarray(L["rgb"], dtype=np.uint8)).save(vdir / f"{L['id']}.jpg", quality=90) + return {"id": L["id"], "name": L["name"], "file": f"{L['id']}.jpg", "legend": L["legend"]} + + +def write_index(vdir: Path, entries: list[dict]) -> None: + vdir.mkdir(parents=True, exist_ok=True) + (vdir / "layers.json").write_text(json.dumps(entries, indent=1)) + + +def write_all(vdir: Path, layers: list[dict]) -> None: + write_index(vdir, [write_layer(vdir, L) for L in layers]) diff --git a/mapgen/zones.py b/mapgen/zones.py new file mode 100644 index 0000000..46386f7 --- /dev/null +++ b/mapgen/zones.py @@ -0,0 +1,38 @@ +"""Config O₂ / low-gravity zones: smooth, noise-warped discs in mask units +(O₂ ×(1 + 0.5 v), gravity ×(1 + 0.7 v)), added to the painted masks in the `sketch` stage. Very gradual: at most +≈ 1.9 · |effect| / radius per km (≤ 0.23 O₂ points and ≤ 0.014 g per 30 km for the configured zones).""" +from __future__ import annotations + +import numpy as np + +from .noise import fbm, name_seed +from .sphere import gc_dist_km, latlon_to_xyz + +WARP = 0.2 # outline wobble: distances stretched or shrunk by up to 20 % (low-frequency noise) +WARP_FREQ = 1.0 + + +def smootherstep(t): + t = np.clip(t, 0.0, 1.0) + return t * t * t * (t * (6.0 * t - 15.0) + 10.0) + + +def contribution(xyz, zone: dict, seed: int, radius_km: float): + """v · (1 − smootherstep(d′ / r)), d′ = d · (1 + WARP · warp): v at the centre, 0 beyond r / (1 − WARP).""" + xyz = np.asarray(xyz, dtype=np.float64) + d = gc_dist_km(xyz, latlon_to_xyz(*zone["center"]), radius_km) + r = float(zone["radius_km"]) + out = np.zeros(len(xyz)) + near = d < r / (1.0 - WARP) + if near.any(): + w = np.clip(2.0 * fbm(xyz[near], name_seed(seed, zone["name"]), 3, WARP_FREQ), -1.0, 1.0) + out[near] = zone["v"] * (1.0 - smootherstep(d[near] * (1.0 + WARP * w) / r)) + return out + + +def masks(xyz, zones: list, seed: int, radius_km: float): + """(m_o2, m_gravity): the config zones summed (the caller adds the painted masks and clips to [−1, 1]).""" + m = {"o2": np.zeros(len(xyz)), "gravity": np.zeros(len(xyz))} + for z in zones: + m[z["field"]] += contribution(xyz, z, seed, radius_km) + return m["o2"], m["gravity"] |