jpmschweitzerandClaude Opus 5.5
76e3f3ec96
D-258 rulings 5 and 6: the pre-seed biome input is class ids only, on two maps, because compute_biome was already computing them apart and throwing one away. It gave every land pixel its Whittaker (climate) class, then painted terrain over it (ocean bands, ice, altitude snow, mountain rock, lava and ash, the dry, lunar and ferric ground). planet_simulation: - compute_biome_layers returns (biome, terrain). Each rule declares the layer it writes, by RULE not by class: ice from the cold climate box is biome, altitude snow is terrain with the lowland biome kept beneath. Rules that place life (thermophiles, mats, crust, the oasis rings) write the biome and clear the terrain under them. - 255 means "nothing on this layer", not 0, because class 0 is ocean_deep, itself terrain. Ids stay literal. - compute_biome is now compose_layers() of the pair, so the renderers are untouched. Proven byte-identical, dtype included, on 65 real bodies: every 8th standard-atmosphere body plus every thin, thick, reducing and trace one, captured before the edit and compared after it. - simulate() carries both maps as biome_layer / terrain_layer. reach atlas planet bake-biome [--body --limit --dry-run --check]: - Writes biomemap.png + terrainmap.png (8-bit, 1024x512, tEXt: layer, none, classes, decision) beside each heightmap. Scope is the heightmap bake set with an atmosphere (257 bodies; airless are deferred, per D-258). - Pre-seed by construction: simulate() is keyed on the body frontmatter's seed, and no world seed is accepted anywhere. - --check re-simulates and compares DECODED PIXELS, never bytes, since PNG encoding drifts across Pillow/zlib versions (T-1291) and a check that trips on that gets muted. Mutation-proven: one flipped pixel exits 1 and names the file; the restored file exits 0. - import_heightmaps' body lookup is extracted as body_def_for and shared, not copied. Tests (make test-tooling): rock keeps its biome beneath, thin-atmosphere ground is terrain with only scattered life, water is terrain only, no pixel is empty on both maps, and stacking equals the rendered grid. Routing rock through the biome map fails two of them by name. The router drift test covers the new verb. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
1012 lines
45 KiB
Python
1012 lines
45 KiB
Python
"""
|
||
planet_simulation.py
|
||
--------------------
|
||
Terrain simulation stack for the Settled Reach planet generator.
|
||
|
||
Consumes a body_definition dict (output of body_definition_parser.py)
|
||
and produces a terrain dict consumed by planet_renderer.render_globe().
|
||
|
||
Output terrain dict:
|
||
{
|
||
"elevation": float32 (H, W) [0, 1] normalised elevation
|
||
"temperature": float32 (H, W) [0, 1] 0=coldest, 1=hottest
|
||
"moisture": float32 (H, W) [0, 1] 0=driest, 1=wettest
|
||
"biome": int8 (H, W) biome class index (composed)
|
||
"biome_layer": uint8 (H, W) climate/life class, 255 = none
|
||
"terrain_layer": uint8 (H, W) landform/substrate class, 255 = none
|
||
"surface_water": bool (H, W) ocean/lake mask
|
||
"hillshade": float32 (H, W) [0, 1] lighting from slope+aspect
|
||
"sea_level": float elevation threshold
|
||
}
|
||
|
||
Pipeline:
|
||
1. Elevation - continent mask + domain-warped FBM + tectonic ridges + erosion
|
||
2. Temperature - analytical formula: star + latitude + altitude
|
||
3. Moisture - Hadley cells + ocean proximity + rain shadow
|
||
4. Hillshade - surface normals from elevation gradient
|
||
5. Rivers - downhill carving from moisture-seeded sources
|
||
6. Biome - extended Whittaker lookup + modifier stack
|
||
|
||
Grid: 512 x 256 (longitude x latitude), equirectangular.
|
||
Row 0 = north pole, row 255 = south pole.
|
||
Col 0 = 180W, col 511 = 180E.
|
||
"""
|
||
|
||
import logging
|
||
import math
|
||
import numpy as np
|
||
from scipy.ndimage import gaussian_filter
|
||
|
||
from tooling.domains.atlas.planet.biome_config import (
|
||
WHITTAKER_TABLE, CLASS_T_BAND, EXOTIC_CLASSES, CRATER_SCALING,
|
||
)
|
||
from tooling.core import console
|
||
|
||
log = logging.getLogger(__name__)
|
||
# Canonical heightmap grid (D-202 amended, #963): bumped to 1024×512 so the
|
||
# stored elevation has real mid-scale detail for the lower cascade layers. The
|
||
# elevation noise is resolution-independent (normalized coords + absolute
|
||
# frequencies), so a higher grid samples the SAME continuous terrain at finer
|
||
# density — features keep their physical size and generation stays deterministic.
|
||
# Pixel-unit operations (gaussian sigma, crater radii, filter windows) scale by
|
||
# `GRID_W / 512` so smoothing/morphology behave identically at any resolution.
|
||
GRID_W = 1024
|
||
GRID_H = 512
|
||
|
||
|
||
def oasis_ring_iterations(grid_w: int) -> tuple[int, int]:
|
||
"""Pixel-unit `binary_dilation` iteration counts for the oasis vegetation
|
||
ring (T-964, PR #210 review — hoshe), scaled by `grid_w / 512` per the
|
||
module's GRID_W pixel-unit convention (see the comment above). Base counts
|
||
(2, 4) are the values authored at the original 512-wide grid; `round(...)`
|
||
is used (not floor/ceil) so the scaling is symmetric around exact
|
||
doubling/halving points, and `max(1, ...)` guards against a degenerate
|
||
zero-iteration dilation at very small `grid_w`.
|
||
|
||
Parameterized on `grid_w` (not reading the module-global `GRID_W`
|
||
directly) so it is independently unit-testable at both the historical
|
||
512 baseline and the current production 1024 value without needing to
|
||
monkeypatch the module global.
|
||
"""
|
||
scale = grid_w / 512.0
|
||
ring1_iters = max(1, round(2 * scale))
|
||
ring2_iters = max(1, round(4 * scale))
|
||
return ring1_iters, ring2_iters
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Seeded RNG
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _rng(seed: int, salt: int = 0) -> np.random.Generator:
|
||
return np.random.default_rng(seed ^ (salt * 2654435761))
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Noise primitives
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _hash2(x: np.ndarray, y: np.ndarray, seed: int) -> np.ndarray:
|
||
s = np.int64(seed & 0xFFFF)
|
||
h = (x.astype(np.int64) * np.int64(1619) +
|
||
y.astype(np.int64) * np.int64(31337) +
|
||
s * np.int64(6971)) & np.int64(0xFFFFFFFF)
|
||
h = ((h >> 16) ^ h) * np.int64(0x45d9f3b) & np.int64(0xFFFFFFFF)
|
||
h = ((h >> 16) ^ h) * np.int64(0x45d9f3b) & np.int64(0xFFFFFFFF)
|
||
return (h & np.int64(0xFFFF)).astype(np.float32) / 65535.0
|
||
|
||
|
||
def _vnoise(u, v, freq, seed):
|
||
"""Standard 2D value noise — NOT seamless. Use _vnoise_s for longitude axis."""
|
||
uf = u * freq; vf = v * freq
|
||
x0 = np.floor(uf).astype(np.int32); y0 = np.floor(vf).astype(np.int32)
|
||
x1 = x0 + 1; y1 = y0 + 1
|
||
tx = uf - x0; ty = vf - y0
|
||
tx = tx * tx * (3.0 - 2.0 * tx)
|
||
ty = ty * ty * (3.0 - 2.0 * ty)
|
||
v00 = _hash2(x0, y0, seed); v10 = _hash2(x1, y0, seed)
|
||
v01 = _hash2(x0, y1, seed); v11 = _hash2(x1, y1, seed)
|
||
return (v00*(1-tx)*(1-ty) + v10*tx*(1-ty) +
|
||
v01*(1-tx)*ty + v11*tx*ty).astype(np.float32)
|
||
|
||
|
||
def _hash3(x, y, z, seed):
|
||
"""Hash for 3D integer coords."""
|
||
s = np.int64(seed & 0xFFFF)
|
||
h = (x.astype(np.int64) * np.int64(1619) +
|
||
y.astype(np.int64) * np.int64(31337) +
|
||
z.astype(np.int64) * np.int64(49979) +
|
||
s * np.int64(6971)) & np.int64(0xFFFFFFFF)
|
||
h = ((h >> 16) ^ h) * np.int64(0x45d9f3b) & np.int64(0xFFFFFFFF)
|
||
h = ((h >> 16) ^ h) * np.int64(0x45d9f3b) & np.int64(0xFFFFFFFF)
|
||
return (h & np.int64(0xFFFF)).astype(np.float32) / 65535.0
|
||
|
||
|
||
def _vnoise_seamless(u, v, freq, seed):
|
||
"""
|
||
Seamless value noise in the U (longitude) axis only.
|
||
Maps u -> (cos(u*2π), sin(u*2π)) before hashing, so the noise
|
||
field is periodic in U with period 1 — no seam at the date line.
|
||
V (latitude) is not periodic — poles are endpoints, not a loop.
|
||
"""
|
||
# Project U onto a circle: (cx, cy)
|
||
# Divide circle radius by 2π so one full revolution spans the same
|
||
# distance as freq units on the flat V axis — corrects aspect ratio.
|
||
angle = u * (2.0 * math.pi)
|
||
r = freq / (2.0 * math.pi)
|
||
cx = np.cos(angle) * r
|
||
cy = np.sin(angle) * r
|
||
vf = v * freq
|
||
|
||
# Integer lattice in 3D (cx, cy, vf)
|
||
x0 = np.floor(cx).astype(np.int32); x1 = x0 + 1
|
||
y0 = np.floor(cy).astype(np.int32); y1 = y0 + 1
|
||
z0 = np.floor(vf).astype(np.int32); z1 = z0 + 1
|
||
|
||
# Smoothstep weights
|
||
tx = cx - x0; tx = tx * tx * (3.0 - 2.0 * tx)
|
||
ty = cy - y0; ty = ty * ty * (3.0 - 2.0 * ty)
|
||
tz = vf - z0; tz = tz * tz * (3.0 - 2.0 * tz)
|
||
|
||
# Trilinear interpolation over 8 corners
|
||
v000 = _hash3(x0, y0, z0, seed); v100 = _hash3(x1, y0, z0, seed)
|
||
v010 = _hash3(x0, y1, z0, seed); v110 = _hash3(x1, y1, z0, seed)
|
||
v001 = _hash3(x0, y0, z1, seed); v101 = _hash3(x1, y0, z1, seed)
|
||
v011 = _hash3(x0, y1, z1, seed); v111 = _hash3(x1, y1, z1, seed)
|
||
|
||
return (v000*(1-tx)*(1-ty)*(1-tz) + v100*tx*(1-ty)*(1-tz) +
|
||
v010*(1-tx)*ty*(1-tz) + v110*tx*ty*(1-tz) +
|
||
v001*(1-tx)*(1-ty)*tz + v101*tx*(1-ty)*tz +
|
||
v011*(1-tx)*ty*tz + v111*tx*ty*tz).astype(np.float32)
|
||
|
||
|
||
def _fbm(u, v, seed, octaves=6, lacunarity=2.0, gain=0.50, base_freq=2.0):
|
||
"""FBM using seamless noise in U — no longitude seam."""
|
||
result = np.zeros_like(u, dtype=np.float32)
|
||
amp = 1.0; freq = base_freq; total = 0.0
|
||
rng = np.random.default_rng(seed)
|
||
for _ in range(octaves):
|
||
oct_seed = int(rng.integers(0, 0x7FFFFFFF))
|
||
result += amp * _vnoise_seamless(u, v, freq, oct_seed)
|
||
total += amp
|
||
amp *= gain; freq *= lacunarity
|
||
return result / (total + 1e-9)
|
||
|
||
|
||
def _domain_warp(u, v, seed, strength=0.35):
|
||
"""Domain warp using seamless FBM — preserves no-seam property."""
|
||
wu = _fbm(u + 1.7, v + 9.2, seed + 1, octaves=4) * 2.0 - 1.0
|
||
wv = _fbm(u + 8.3, v + 2.8, seed + 2, octaves=4) * 2.0 - 1.0
|
||
# Only warp u periodically — keep v warp non-periodic (poles stay poles)
|
||
return (u + wu * strength) % 1.0, np.clip(v + wv * strength * 0.5, 0.0, 1.0)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Coordinate grids
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _make_grids():
|
||
u_1d = np.linspace(0, 1, GRID_W, dtype=np.float32)
|
||
v_1d = np.linspace(0, 1, GRID_H, dtype=np.float32)
|
||
u, v = np.meshgrid(u_1d, v_1d)
|
||
lat_frac = -(v - 0.5) * 2.0 # +1 = north, -1 = south
|
||
lon_frac = (u - 0.5) * 2.0
|
||
lat_rad = lat_frac * (math.pi / 2.0)
|
||
return u, v, lat_frac, lon_frac, lat_rad
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 1. Elevation
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _continent_mask(u, v, seed, land_fraction):
|
||
def _norm(a):
|
||
lo, hi = a.min(), a.max()
|
||
return (a - lo) / (hi - lo + 1e-9)
|
||
|
||
def _contrast(a, strength=3.0):
|
||
"""
|
||
S-curve contrast: pushes highs toward 1 and lows toward 0
|
||
regardless of the field mean. More reliable than power curves
|
||
which behave differently depending on the field's distribution.
|
||
strength controls steepness — higher = sharper separation.
|
||
"""
|
||
# Sigmoid centred at 0.5: f(x) = 1/(1+exp(-k*(x-0.5)))
|
||
k = strength * 8.0
|
||
return 1.0 / (1.0 + np.exp(-k * (a - 0.5)))
|
||
|
||
# Primary: large continental plates
|
||
wu1, wv1 = _domain_warp(u, v, seed, strength=0.45)
|
||
primary = _norm(_fbm(wu1, wv1, seed + 10, octaves=5, gain=0.58, base_freq=1.2))
|
||
|
||
# Secondary: independent medium-scale field.
|
||
# S-curve contrast gives reliable highs and lows regardless of seed.
|
||
wu2, wv2 = _domain_warp(u, v, seed + 11, strength=0.40)
|
||
sec_raw = _norm(_fbm(wu2, wv2, seed + 20, octaves=5, gain=0.55, base_freq=1.8))
|
||
secondary = _contrast(sec_raw, strength=2.5)
|
||
|
||
# Rift: anisotropic thin elongated features
|
||
wu3, wv3 = _domain_warp(u, v, seed + 17, strength=0.30)
|
||
rift = _norm(_fbm(wu3, wv3 * 0.35, seed + 30, octaves=4, gain=0.52, base_freq=3.5))
|
||
|
||
# Multiplicative gate: secondary zeroes kill primary → ocean channels
|
||
separated = primary * (0.4 + secondary * 0.6)
|
||
combined = separated * 0.82 + (rift - 0.5) * 0.18
|
||
|
||
return _norm(combined).astype(np.float32)
|
||
|
||
|
||
def _tectonic_ridges(u, v, seed, n_plates=8):
|
||
rng = _rng(seed, 99)
|
||
px = rng.uniform(0, 1, n_plates).astype(np.float32)
|
||
py = rng.uniform(0, 1, n_plates).astype(np.float32)
|
||
H, W = u.shape
|
||
|
||
# Domain-warp coords before Voronoi — bends ridge positions into curves
|
||
wu1 = _fbm(u * 1.5 + 3.1, v * 1.5 + 7.4, seed + 201, octaves=3,
|
||
gain=0.55, base_freq=1.8) * 2.0 - 1.0
|
||
wv1 = _fbm(u * 1.5 + 8.6, v * 1.5 + 2.2, seed + 202, octaves=3,
|
||
gain=0.55, base_freq=1.8) * 2.0 - 1.0
|
||
wu2 = _fbm(u * 4.0 + 1.3, v * 4.0 + 5.7, seed + 203, octaves=2,
|
||
gain=0.50, base_freq=3.5) * 2.0 - 1.0
|
||
wv2 = _fbm(u * 4.0 + 6.1, v * 4.0 + 0.9, seed + 204, octaves=2,
|
||
gain=0.50, base_freq=3.5) * 2.0 - 1.0
|
||
|
||
uw = (u + wu1 * 0.22 + wu2 * 0.08) % 1.0
|
||
vw = np.clip(v + wv1 * 0.18 + wv2 * 0.06, 0.0, 1.0)
|
||
|
||
dist1 = np.full((H, W), np.inf, dtype=np.float32)
|
||
dist2 = np.full((H, W), np.inf, dtype=np.float32)
|
||
for i in range(n_plates):
|
||
du = np.minimum(np.abs(uw - px[i]), 1.0 - np.abs(uw - px[i]))
|
||
dv = np.abs(vw - py[i])
|
||
d = np.sqrt(du**2 + dv**2)
|
||
mask = d < dist1
|
||
dist2 = np.where(mask, dist1, np.minimum(dist2, d))
|
||
dist1 = np.where(mask, d, dist1)
|
||
|
||
# Two ridge widths: broad ranges + sharp collision zones
|
||
broad = np.exp(-((dist2 - dist1) / 0.06) ** 2) * 0.5
|
||
sharp = np.exp(-((dist2 - dist1) / 0.025) ** 2) * 1.0
|
||
ridge_raw = np.clip(broad + sharp, 0, 1)
|
||
|
||
# Amplitude variation along ridge
|
||
ridge_noise = _fbm(u, v, seed + 50, octaves=4, gain=0.55, base_freq=4.0)
|
||
|
||
# Fracture zones — cross-cutting features (transform faults, rift valleys)
|
||
# Anisotropic: stretch u relative to v for elongated cross features
|
||
fracture = _fbm(u * 0.4, v, seed + 77, octaves=3, gain=0.6, base_freq=6.0)
|
||
fracture = np.clip(fracture - 0.55, 0, 1) * 2.0
|
||
|
||
return np.clip(ridge_raw * (0.35 + 0.65 * ridge_noise)
|
||
+ fracture * 0.20, 0, 1).astype(np.float32)
|
||
|
||
|
||
def _erode(terrain, passes, seed):
|
||
result = terrain.copy()
|
||
# Resolution scale: smoothing radii and the per-pixel slope (which halves as
|
||
# the grid doubles, since np.gradient is in pixel units) scale with width so
|
||
# erosion behaves identically at any GRID size.
|
||
scale = result.shape[1] / 512.0
|
||
for _ in range(passes):
|
||
gy, gx = np.gradient(result)
|
||
slope = np.sqrt(gx**2 + gy**2)
|
||
smooth = gaussian_filter(result, sigma=1.2 * scale)
|
||
weight = np.clip(slope * 6.0 * scale, 0.0, 1.0)
|
||
result = result * (1.0 - weight * 0.35) + smooth * (weight * 0.35)
|
||
gy, gx = np.gradient(result)
|
||
slope = np.sqrt(gx**2 + gy**2)
|
||
flow = gaussian_filter(slope, sigma=3.0 * scale)
|
||
flow = (flow - flow.min()) / (flow.max() - flow.min() + 1e-9)
|
||
result = result - flow * 0.06
|
||
return np.clip(result, 0.0, 1.0)
|
||
|
||
|
||
def compute_elevation(body_def, u, v, lat_frac):
|
||
seed = body_def["seed"]
|
||
planet_class = body_def["planet_class"].replace("_ringed", "")
|
||
land_frac = body_def["terrain"]["land_fraction"]
|
||
tectonics = body_def["terrain"].get("tectonics", "active")
|
||
|
||
plate_map = {"extreme": 12, "active": 8, "low": 5, "none": 3}
|
||
erosion_map = {"extreme": 1, "active": 3, "low": 4, "none": 2}
|
||
n_plates = plate_map.get(tectonics, 8)
|
||
erosion_p = erosion_map.get(tectonics, 3)
|
||
|
||
ocean_pct = (1.0 - land_frac) * 100.0
|
||
detail = _fbm(u, v, seed + 300, octaves=5, gain=0.45, base_freq=4.0)
|
||
|
||
if tectonics == "none":
|
||
# No tectonic activity: gentle base terrain, no ridges, no continents.
|
||
# Craters dominate on these worlds.
|
||
base = _fbm(u, v, seed + 100, octaves=4, gain=0.50, base_freq=1.5)
|
||
elev = base * 0.60 + detail * 0.40
|
||
else:
|
||
# Tectonic worlds: continent mask + ridges scaled by activity level.
|
||
cont = _continent_mask(u, v, seed, land_frac)
|
||
ridges = _tectonic_ridges(u, v, seed, n_plates=n_plates)
|
||
|
||
sea_level_est = float(np.percentile(cont, ocean_pct))
|
||
land_mask = cont >= sea_level_est
|
||
|
||
# Ridge prominence scales with tectonic activity
|
||
ridge_weight = {"low": 0.12, "active": 0.25, "extreme": 0.38}
|
||
rw = ridge_weight.get(tectonics, 0.25)
|
||
|
||
elev = (cont * (0.80 - rw)
|
||
+ ridges * rw * land_mask
|
||
+ detail * 0.20)
|
||
|
||
# Craters happen everywhere. Atmosphere controls how many impactors
|
||
# survive entry; tectonics controls how many craters get resurfaced.
|
||
# Both reduce density, neither toggles craters off entirely.
|
||
crater_factor = (CRATER_SCALING["atmosphere"].get(body_def["physical"]["atmosphere"], 0.25)
|
||
* CRATER_SCALING["tectonics"].get(tectonics, 0.3))
|
||
|
||
if crater_factor > 0.02:
|
||
rng = _rng(seed, 77)
|
||
base_count = CRATER_SCALING["base_count"]
|
||
n_craters = max(5, int(base_count * crater_factor))
|
||
cy_c = rng.uniform(0, GRID_H, n_craters).astype(np.float32)
|
||
cx_c = rng.uniform(0, GRID_W, n_craters).astype(np.float32)
|
||
# Power-law: most craters are small, a few large. Radii are in pixels,
|
||
# so scale with resolution to keep craters the same physical size.
|
||
raw_sizes = rng.power(0.4, n_craters) # skewed toward 0
|
||
sizes = ((2 + raw_sizes * 28) * (GRID_W / 512.0)).astype(np.float32)
|
||
# Depth scales with crater factor — eroded worlds have shallower craters
|
||
depth_scale = 0.5 + 0.5 * crater_factor
|
||
depths = ((0.05 + raw_sizes * 0.20) * depth_scale).astype(np.float32)
|
||
|
||
rows = np.arange(GRID_H, dtype=np.float32)
|
||
cols = np.arange(GRID_W, dtype=np.float32)
|
||
rr, cc = np.meshgrid(rows, cols, indexing='ij')
|
||
# Latitude correction: scale longitude distance by cos(lat) so
|
||
# craters are circular on the sphere, not stretched at the poles.
|
||
lat_rad = (0.5 - rr / GRID_H) * math.pi # +pi/2 at north, -pi/2 at south
|
||
cos_lat = np.cos(lat_rad)
|
||
cos_lat = np.clip(cos_lat, 0.1, 1.0) # avoid division issues at poles
|
||
craters = np.zeros_like(elev)
|
||
for i in range(n_craters):
|
||
dy = rr - cy_c[i]
|
||
dx = cc - cx_c[i]
|
||
# Wrap longitude for craters near the date line
|
||
dx = np.minimum(np.abs(dx), GRID_W - np.abs(dx))
|
||
# Scale dx by cos(lat) at the crater center
|
||
center_lat = (0.5 - cy_c[i] / GRID_H) * math.pi
|
||
dx_scaled = dx / max(math.cos(center_lat), 0.1)
|
||
d = np.sqrt(dy**2 + dx_scaled**2)
|
||
r = sizes[i]
|
||
dep = depths[i]
|
||
# Crater profile: flat floor inside 0.6r, raised rim at 0.9-1.1r,
|
||
# smooth falloff outside. More realistic than gaussian dimple.
|
||
floor = np.clip(1.0 - d / (r * 0.6), 0, 1)
|
||
rim = np.exp(-((d - r) / (r * 0.25))**2)
|
||
craters -= dep * floor * 0.8 # excavate floor
|
||
craters += dep * rim * 0.3 # raise rim
|
||
elev = elev + craters
|
||
elev = np.clip(elev, 0.0, None) # floor at 0
|
||
elif planet_class == "frozen":
|
||
elev = gaussian_filter(elev, sigma=1.5 * (GRID_W / 512.0)).astype(np.float32)
|
||
elif planet_class == "volcanic":
|
||
erosion_p = max(1, erosion_p - 1)
|
||
|
||
elev = _erode(elev, passes=erosion_p, seed=seed)
|
||
|
||
lo, hi = elev.min(), elev.max()
|
||
elev = (elev - lo) / (hi - lo + 1e-9)
|
||
sea_level = float(np.percentile(elev, ocean_pct))
|
||
|
||
# Polar ice — smooth land elevation toward a low plateau at high latitudes.
|
||
# Only applies when there's an atmosphere to deliver precipitation/ice.
|
||
# Airless bodies have no polar caps — cold rock stays rock.
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "ocean")
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
has_polar_ice = atmo not in ("none",) and hydro not in ("none", "subsurface")
|
||
|
||
if has_polar_ice:
|
||
ice_lat = body_def["terrain"].get("polar_ice_lat", 0.80)
|
||
lat_abs = np.abs(lat_frac)
|
||
ice_blend = np.clip((lat_abs - ice_lat) / (1.0 - ice_lat + 0.01), 0, 1)
|
||
if planet_class == "frozen":
|
||
ice_blend = np.clip(ice_blend * 2.0, 0, 1)
|
||
land_mask = elev >= sea_level
|
||
ice_target = sea_level + 0.05
|
||
elev = np.where(
|
||
land_mask,
|
||
elev * (1.0 - ice_blend * 0.6) + ice_target * (ice_blend * 0.6),
|
||
elev)
|
||
elev = np.clip(elev, 0.0, 1.0).astype(np.float32)
|
||
|
||
sea_level = float(np.percentile(elev, ocean_pct))
|
||
# Dry worlds (no/subsurface hydrosphere): low elevation is dry basin, not ocean.
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "ocean")
|
||
if hydro in ("none", "subsurface"):
|
||
surf_water = np.zeros_like(elev, dtype=bool)
|
||
else:
|
||
surf_water = elev < sea_level
|
||
return elev, sea_level, surf_water
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 2. Temperature
|
||
# ---------------------------------------------------------------------------
|
||
|
||
# Stellar luminosity relative to Sol (approximate midpoint per spectral type)
|
||
STAR_LUMINOSITY = {
|
||
"O": 100000.0, "B": 1000.0, "A": 10.0,
|
||
"F": 2.5, "G": 1.0, "K": 0.4, "M": 0.04,
|
||
}
|
||
|
||
|
||
# CLASS_T_BAND loaded from biomes.toml via biome_config
|
||
|
||
def compute_temperature(body_def, elevation, sea_level, lat_frac):
|
||
star_type = body_def["star"]["type"]
|
||
distance_au = body_def["orbit"]["distance_au"]
|
||
axial_tilt = body_def["orbit"]["axial_tilt_deg"]
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
planet_class = body_def["planet_class"].replace("_ringed", "")
|
||
geothermal = body_def.get("environment", {}).get("geothermal_flux", "low")
|
||
|
||
# Equilibrium temperature — descriptor-anchored.
|
||
#
|
||
# We compute the raw stellar physics (Stefan-Boltzmann) to get a
|
||
# physically grounded value, then clamp it to the temperature band
|
||
# appropriate for the planet_class. This ensures the wiki's descriptors
|
||
# (temperate, frozen, arid…) are always honoured even when orbital
|
||
# parameters were set with "close enough" precision.
|
||
#
|
||
# Within the clamped band, the raw value still drives relative warmth:
|
||
# a close-in temperate world sits at the warm end of the temperate band,
|
||
# a far-out one at the cool end. The fiction wins; physics sets the gradient.
|
||
lum = body_def.get("star", {}).get("luminosity_solar",
|
||
STAR_LUMINOSITY.get(star_type, 1.0))
|
||
t_raw = 278.5 * (lum ** 0.25) / math.sqrt(max(distance_au, 0.01))
|
||
|
||
greenhouse = {"none": 0, "thin": 8, "standard": 33, "thick": 80}
|
||
t_raw += greenhouse.get(atmo, 0)
|
||
|
||
|
||
temperature_clamped = False
|
||
temperature_raw_K = float(t_raw)
|
||
|
||
if planet_class in CLASS_T_BAND:
|
||
t_lo, t_hi = CLASS_T_BAND[planet_class]
|
||
t_base = float(np.clip(t_raw, t_lo, t_hi))
|
||
if t_raw < t_lo or t_raw > t_hi:
|
||
temperature_clamped = True
|
||
log.debug(f" T_raw={t_raw:.0f}K clamped to [{t_lo},{t_hi}] "
|
||
f"for {planet_class} ({body_def.get('id','')})")
|
||
else:
|
||
t_base = t_raw
|
||
|
||
tilt_factor = 1.0 - (axial_tilt / 90.0) * 0.5
|
||
# Atmosphere controls heat redistribution — thicker atmo = smaller
|
||
# equator-pole gradient. Thin/no atmo = extreme day/night but we
|
||
# still want the planet class to read correctly at the poles.
|
||
atmo_gradient_scale = {"none": 0.6, "thin": 0.7, "standard": 1.0, "thick": 1.2}
|
||
lat_gradient = 60.0 * tilt_factor * atmo_gradient_scale.get(atmo, 1.0)
|
||
t_lat = t_base - lat_gradient * np.abs(lat_frac)
|
||
|
||
max_relief_km = body_def.get("terrain", {}).get("max_elevation_km", 10.0)
|
||
elev_land = np.where(elevation >= sea_level,
|
||
(elevation - sea_level) / (1.0 - sea_level + 1e-9), 0.0)
|
||
elev_km = elev_land * max_relief_km
|
||
lapse = 6.5 if atmo != "none" else 2.0
|
||
t_final = t_lat - lapse * elev_km
|
||
|
||
class_offset = {"frozen": -30, "volcanic": 20, "arid": 10}
|
||
t_final += class_offset.get(planet_class, 0)
|
||
|
||
geo_boost = {"low": 0, "moderate": 5, "high": 15, "extreme": 35}
|
||
t_final += geo_boost.get(geothermal, 0)
|
||
|
||
# Soft floor: prevent planet class from being contradicted at the poles.
|
||
# An arid world shouldn't have ice caps; a volcanic world shouldn't freeze.
|
||
# Clamp the minimum temperature to the class band's lower bound.
|
||
if planet_class in CLASS_T_BAND:
|
||
t_floor = CLASS_T_BAND[planet_class][0]
|
||
t_final = np.maximum(t_final, t_floor)
|
||
|
||
# Return absolute Kelvin grid plus audit metadata.
|
||
# Biome lookup needs absolute values; renderer normalises for display.
|
||
return t_final.astype(np.float32), temperature_clamped, temperature_raw_K
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 3. Moisture
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def compute_moisture(body_def, elevation, sea_level, temperature,
|
||
lat_frac, lon_frac):
|
||
# Normalise temperature locally for moisture computation
|
||
t_norm = np.clip((temperature - temperature.min()) /
|
||
(temperature.max() - temperature.min() + 1e-9), 0, 1)
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
planet_class = body_def["planet_class"].replace("_ringed", "")
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "none")
|
||
|
||
if atmo == "none":
|
||
return np.zeros((GRID_H, GRID_W), dtype=np.float32)
|
||
|
||
lat_abs = np.abs(lat_frac)
|
||
|
||
# Hadley cell bands
|
||
itcz = np.clip(1.0 - (lat_abs / 0.33), 0, 1)
|
||
subtr = np.clip(1.0 - np.abs(lat_abs - 0.50) / 0.17, 0, 1)
|
||
polar = np.clip((lat_abs - 0.67) / 0.33, 0, 1)
|
||
hadley = np.clip(itcz * 0.85 + subtr * 0.10 + polar * 0.40, 0, 1)
|
||
|
||
# Ocean proximity
|
||
surf_water = elevation < sea_level
|
||
if surf_water.any():
|
||
from scipy.ndimage import distance_transform_edt
|
||
dist = distance_transform_edt(~surf_water).astype(np.float32)
|
||
ocean_prox = 1.0 - np.clip(dist / (dist.max() * 0.5 + 1e-9), 0, 1)
|
||
else:
|
||
ocean_prox = np.zeros((GRID_H, GRID_W), dtype=np.float32)
|
||
|
||
# Rain shadow — westerly winds: windward (west face) is wet
|
||
shift = max(1, GRID_W // 80)
|
||
elev_above = np.clip(elevation - sea_level, 0, None)
|
||
elev_sh = np.clip(np.roll(elevation, shift, axis=1) - sea_level, 0, None)
|
||
shadow_raw = np.clip(elev_sh - elev_above * 0.5, 0, None)
|
||
shadow_raw = shadow_raw / (shadow_raw.max() + 1e-9)
|
||
rain_shadow = 1.0 - shadow_raw * 0.70
|
||
|
||
moisture = (hadley * 0.40
|
||
+ ocean_prox * 0.45
|
||
+ t_norm * 0.15) * rain_shadow
|
||
|
||
class_scale = {
|
||
"arid": 0.25, "oceanic": 1.30, "forest": 1.30,
|
||
"frozen": 0.55, "volcanic": 0.40, "barren": 0.05,
|
||
}
|
||
moisture *= class_scale.get(planet_class, 1.0)
|
||
|
||
hydro_scale = {
|
||
"ocean": 1.2, "liquid_water": 1.2,
|
||
"subsurface": 0.1, "none": 0.05,
|
||
}
|
||
moisture *= hydro_scale.get(hydro, 1.0)
|
||
|
||
moisture = gaussian_filter(moisture.astype(np.float32), sigma=2.0 * (GRID_W / 512.0))
|
||
# Only normalize if the raw range is substantial — otherwise the
|
||
# normalization re-inflates near-zero moisture on dry worlds back to [0,1].
|
||
m_min, m_max = moisture.min(), moisture.max()
|
||
if m_max > 0.05:
|
||
moisture = ((moisture - m_min) / (m_max - m_min + 1e-9)).astype(np.float32)
|
||
else:
|
||
# Effectively dry — clamp to near-zero
|
||
moisture = np.clip(moisture / 0.05, 0, 1).astype(np.float32)
|
||
return moisture
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 4. Hillshade
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def compute_hillshade(elevation,
|
||
sun_azimuth_deg=315.0,
|
||
sun_altitude_deg=45.0):
|
||
scale = GRID_W / 8.0
|
||
gy, gx = np.gradient(elevation * scale)
|
||
mag = np.sqrt(gx**2 + gy**2 + 1.0)
|
||
nx = -gx / mag; ny = -gy / mag; nz = 1.0 / mag
|
||
|
||
az = math.radians(sun_azimuth_deg)
|
||
alt = math.radians(sun_altitude_deg)
|
||
lx = math.cos(alt) * math.cos(az)
|
||
ly = math.cos(alt) * math.sin(az)
|
||
lz = math.sin(alt)
|
||
|
||
diffuse = np.clip(nx * lx + ny * ly + nz * lz, 0.0, 1.0)
|
||
return (0.25 + 0.75 * diffuse).astype(np.float32)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Rivers are NOT computed here (D-208, #963): river networks are derived by the
|
||
# Rust cascade's D8 drainage from the heightmap — the single source of river
|
||
# truth, with mouths that reach the sea by construction. The old heuristic
|
||
# `compute_rivers` was removed to avoid implying the Python sim owns rivers.
|
||
# ---------------------------------------------------------------------------
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 6. Biome
|
||
# ---------------------------------------------------------------------------
|
||
|
||
# WHITTAKER_TABLE, EXOTIC_CLASSES loaded from biomes.toml via biome_config
|
||
|
||
|
||
# "Nothing on this layer" in the biome and terrain maps. Not 0: class 0 is
|
||
# ocean_deep, itself a terrain class, and every id stays literal (D-258).
|
||
LAYER_NONE = 255
|
||
|
||
|
||
def compose_layers(biome, terrain):
|
||
"""The single class grid the renderers draw: terrain where there is any,
|
||
the biome beneath it everywhere else."""
|
||
return np.where(terrain != LAYER_NONE, terrain, biome).astype(np.int8)
|
||
|
||
|
||
def compute_biome(body_def, elevation, sea_level, surface_water,
|
||
temperature, moisture):
|
||
"""The composed class grid — what the reliefmap and globe renders colour."""
|
||
return compose_layers(*compute_biome_layers(
|
||
body_def, elevation, sea_level, surface_water, temperature, moisture))
|
||
|
||
|
||
def compute_biome_layers(body_def, elevation, sea_level, surface_water,
|
||
temperature, moisture):
|
||
"""Biome and terrain as two maps (D-258, 2026-09-25 rulings 5 and 6).
|
||
|
||
The biome map is the climate and life class; the terrain map is the
|
||
landform or substrate a terrain rule laid over it. Each rule below writes
|
||
the layer it belongs to, by RULE rather than by class: ice from the cold
|
||
climate box is biome, snow on a tropical summit is terrain with the tropical
|
||
biome kept beneath it. A biome rule that lands on terrain (a crust of life
|
||
on bare rock, a shrub patch on a geothermal spot) replaces the terrain there,
|
||
exactly as it replaced it in the single grid, so compose_layers() of the
|
||
result is byte-identical to what compute_biome returned before the split.
|
||
|
||
Both maps use LAYER_NONE for "nothing on this layer": water and
|
||
thin-atmosphere ground carry no biome, most land carries no terrain.
|
||
"""
|
||
H, W = elevation.shape
|
||
biome = np.full((H, W), LAYER_NONE, dtype=np.uint8)
|
||
terrain = np.full((H, W), LAYER_NONE, dtype=np.uint8)
|
||
land = ~surface_water
|
||
|
||
def paint_terrain(mask, cls):
|
||
terrain[mask] = cls
|
||
|
||
def paint_biome(mask, cls):
|
||
biome[mask] = cls
|
||
terrain[mask] = LAYER_NONE
|
||
|
||
def composed():
|
||
return np.where(terrain != LAYER_NONE, terrain, biome)
|
||
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
|
||
# --- Atmosphere gate ---
|
||
# Worlds with no or thin atmosphere can't support vegetation.
|
||
# Skip the Whittaker table entirely — classify by elevation and
|
||
# temperature only, using rock/dust/ice classes.
|
||
if atmo in ("none", "thin"):
|
||
# Dry terrain classes: 27=dust plain, 28=rocky highland,
|
||
# 29=warm dust, 30=cold rock. No vegetation possible.
|
||
# No ice on airless worlds — cold rock stays rock.
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "none")
|
||
has_ice_source = hydro not in ("none", "subsurface") or atmo == "thin"
|
||
|
||
elev_norm = np.where(land,
|
||
(elevation - sea_level) / (1.0 - sea_level + 1e-9),
|
||
0.0)
|
||
cf = np.full(land.sum(), 28, dtype=np.int8) # default: rocky highland
|
||
tf = temperature[land].ravel()
|
||
en = elev_norm[land].ravel()
|
||
|
||
# Moon vs planet: moons use grey lunar palette, planets use warm rock
|
||
is_lunar = body_def.get("body_type") == "moon"
|
||
|
||
if is_lunar:
|
||
# Lunar classes: 31=highland, 32=mare (dark basin), 33=midland
|
||
cf[:] = 33 # default: midland grey
|
||
cf[en > 0.50] = 31 # highland
|
||
cf[en < 0.20] = 32 # mare (dark basin floor)
|
||
if has_ice_source:
|
||
cf[tf < 200] = 17 # ice (only if water source)
|
||
else:
|
||
# Temperature-based classification using dry terrain classes
|
||
if has_ice_source:
|
||
cf[tf < 200] = 17 # ice/snow (only if water source)
|
||
else:
|
||
cf[tf < 200] = 30 # cold rock (no water = no ice)
|
||
cf[(tf >= 200) & (tf < 260)] = 30 # cold rock
|
||
cf[(tf >= 260) & (tf < 310)] = 28 # rocky highland
|
||
cf[(tf >= 310) & (tf < 340)] = 29 # warm dust
|
||
cf[tf >= 340] = 15 # hot desert (scorched)
|
||
|
||
# Elevation variation
|
||
if has_ice_source:
|
||
cf[(en > 0.70) & (tf < 273)] = 17 # high + cold = ice cap
|
||
cf[(en < 0.20) & (tf >= 260)] = 27 # low elevation = dust plain
|
||
|
||
# Bare ground, no life: terrain, with no biome beneath it.
|
||
paint_terrain(land, cf)
|
||
else:
|
||
# --- Standard Whittaker lookup for breathable/toxic atmospheres ---
|
||
tf = temperature[land].ravel()
|
||
mf = moisture[land].ravel()
|
||
cf = np.full(tf.shape, 17, dtype=np.int8) # default: ice
|
||
|
||
# Temperature fed to biome is absolute Kelvin — compare directly
|
||
for (tlo, thi, mlo, mhi, cls) in WHITTAKER_TABLE:
|
||
mask = (tf >= tlo) & (tf <= thi) & (mf >= mlo) & (mf <= mhi)
|
||
cf[mask] = cls
|
||
|
||
biome[land] = cf # the climate class — the biome map's whole reason to exist
|
||
|
||
# Ocean depth bands
|
||
if surface_water.any():
|
||
depth = np.clip((sea_level - elevation) / (sea_level + 1e-9), 0, 1)
|
||
paint_terrain(surface_water & (depth < 0.15), 2)
|
||
paint_terrain(surface_water & (depth >= 0.15) & (depth < 0.50), 1)
|
||
paint_terrain(surface_water & (depth >= 0.50), 0)
|
||
|
||
# Frozen ocean — override ocean biome with ice shelf (class 26).
|
||
# Distinct from land ice (17) — slightly different appearance,
|
||
# blue tint suggests ocean beneath.
|
||
# Add noise to the freeze threshold so the boundary isn't a straight
|
||
# latitude line — ice edges are irregular in reality.
|
||
seed = body_def["seed"]
|
||
u_grid, v_grid, _, _, _ = _make_grids()
|
||
ice_noise = _fbm(u_grid, v_grid, seed + 900, octaves=4,
|
||
gain=0.5, base_freq=3.0) * 2.0 - 1.0
|
||
freeze_threshold = 271.0 + ice_noise * 8.0 # ±8K variation
|
||
frozen_ocean = surface_water & (temperature < freeze_threshold)
|
||
paint_terrain(frozen_ocean, 26)
|
||
|
||
# Very cold override — only on worlds with atmosphere (ice needs deposition).
|
||
# Biome, not terrain: it is the climate's own verdict, the same class the
|
||
# Whittaker ice box gives, so it carries that climate.
|
||
if atmo not in ("none",):
|
||
paint_biome((temperature < 243.0) & land, 17) # below -30C → ice
|
||
|
||
# Elevation overrides — mountain rock and permanent snow.
|
||
# Only apply snow on worlds with atmosphere (ice needs deposition).
|
||
elev_norm = np.where(land,
|
||
(elevation - sea_level) / (1.0 - sea_level + 1e-9),
|
||
0.0)
|
||
hydro_here = body_def.get("environment", {}).get("hydrosphere", "none")
|
||
has_ice_deposition = atmo not in ("none",) and hydro_here not in ("none", "subsurface")
|
||
if has_ice_deposition:
|
||
# Altitude snow is terrain: a glacier with the lowland biome beneath.
|
||
paint_terrain(land & (elev_norm > 0.85), 17)
|
||
# `temperature` is absolute kelvin here, but the two rules below were
|
||
# written against a 0-1 scale and compared kelvin to 0.35 from the first
|
||
# commit: mountain rock could never fire, and the sulfuric test was always
|
||
# true (T-1295). 0.35 only means something as a fraction of this world's own
|
||
# range — the same normalisation simulate() applies for the renderer — so
|
||
# that is what they now compare against. A world with no range reads 0.
|
||
t_span = float(temperature.max() - temperature.min())
|
||
t_norm = (temperature - temperature.min()) / t_span if t_span > 0 else np.zeros_like(temperature)
|
||
paint_terrain(land & (elev_norm > 0.65) & (t_norm < 0.35), 18)
|
||
|
||
# ── Modifier stack ─────────────────────────────────────────────────────
|
||
env = body_def.get("environment", {})
|
||
geothermal = env.get("geothermal_flux", "low")
|
||
chemosyn = env.get("chemosynthetic", False)
|
||
uv_index = env.get("uv_index", "moderate")
|
||
substrate = env.get("substrate", "silicate")
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
planet_class = body_def["planet_class"].replace("_ringed", "")
|
||
|
||
# Geothermal: volcanic worlds get lava/ash at high elevations
|
||
if geothermal in ("extreme", "high") and planet_class == "volcanic":
|
||
paint_terrain(land & (elev_norm > 0.75), EXOTIC_CLASSES["lava_field"])
|
||
paint_terrain(land & (elev_norm > 0.45) & (elev_norm <= 0.75), EXOTIC_CLASSES["ash_field"])
|
||
|
||
# The four rules below place LIFE (thermophiles, mats, crusts, scrub), so
|
||
# they write the biome map and replace any terrain under them.
|
||
|
||
# Thermophilic fields near heat vents on any high-geothermal world
|
||
if geothermal in ("extreme", "high") and not chemosyn:
|
||
hot = (temperature > 303.0) & land & (elev_norm < 0.45)
|
||
paint_biome(hot, EXOTIC_CLASSES["thermophilic_field"])
|
||
|
||
# Chemosynthetic worlds (Europa-type): cold surface, geothermal warmth
|
||
if chemosyn:
|
||
geo_warm = (temperature > 263.0) & (temperature < 293.0) & land
|
||
paint_biome(geo_warm, EXOTIC_CLASSES["chemosynthetic_mat"])
|
||
|
||
# UV radiation: cryptobiotic crust on exposed terrain with thin/no atmo
|
||
if uv_index in ("extreme", "high") and atmo in ("none", "thin"):
|
||
shown = composed()
|
||
exposed = (land & (elev_norm > 0.15) & (elev_norm < 0.65)
|
||
& (moisture < 0.30)
|
||
& (shown != 17) & (shown != 18) & (shown != 19))
|
||
paint_biome(exposed, EXOTIC_CLASSES["cryptobiotic_crust"])
|
||
|
||
# Sulfuric substrate: scrub on volcanic mid-elevations. Unreachable as
|
||
# written — it converts rock (18) only below elev_norm 0.65, and rock exists
|
||
# only above it (T-1295; left for the balancing pass).
|
||
if substrate == "sulfuric":
|
||
scrub = land & (elev_norm > 0.25) & (elev_norm < 0.65) & (t_norm > 0.35)
|
||
paint_biome(scrub & (composed() == 18), EXOTIC_CLASSES["sulfuric_scrub"])
|
||
|
||
# ── Anomaly scatter ─────────────────────────────────────────────────
|
||
# Sparse micro-features that break biome uniformity and tell stories.
|
||
# A high-frequency noise field selects ~2-5% of cells for anomaly
|
||
# replacement. The anomaly type depends on the surrounding biome context.
|
||
if atmo not in ("none",):
|
||
seed = body_def["seed"]
|
||
u_grid, v_grid, _, _, _ = _make_grids()
|
||
scatter_noise = _fbm(u_grid, v_grid, seed + 800, octaves=3,
|
||
gain=0.6, base_freq=12.0)
|
||
# High threshold = sparse features (~3% of land)
|
||
scatter_mask = (scatter_noise > 0.72) & land
|
||
|
||
if scatter_mask.any():
|
||
b_local = composed()[scatter_mask]
|
||
t_local = temperature[scatter_mask]
|
||
m_local = moisture[scatter_mask]
|
||
e_local = elev_norm[scatter_mask]
|
||
new_b = b_local.copy()
|
||
|
||
# Temperate/forest → volcanic vent (lava at high elevation)
|
||
veg_mask = np.isin(b_local, [5, 6, 7, 8, 9, 10, 11])
|
||
new_b[veg_mask & (e_local > 0.50)] = 19 # lava field
|
||
new_b[veg_mask & (e_local > 0.35) & (e_local <= 0.50)] = 25 # ash
|
||
|
||
# Desert/dry → oasis with vegetation ring (only if water exists)
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "none")
|
||
has_water = hydro not in ("none", "subsurface")
|
||
dry_mask = np.isin(b_local, [13, 14, 15, 27, 28, 29])
|
||
if has_water:
|
||
# Very rare lake in desert lowlands
|
||
new_b[dry_mask & (e_local < 0.10) & (m_local > 0.20)] = 2 # shallow water
|
||
# Vegetation around moisture (oasis fringe — works even without
|
||
# standing water, represents subsurface moisture reaching roots)
|
||
new_b[dry_mask & (m_local > 0.15) & (e_local >= 0.10)] = 7 # savanna
|
||
|
||
# Frozen → geothermal hotspot with pioneer vegetation
|
||
cold_mask = np.isin(b_local, [16, 17])
|
||
new_b[cold_mask & (t_local > 260)] = 12 # shrubland (hardy plants)
|
||
|
||
# Volcanic → cooling zone with pioneer life
|
||
lava_mask = np.isin(b_local, [19, 25])
|
||
new_b[lava_mask & (t_local < 310) & (m_local > 0.30)] = 24 # lithic pioneer
|
||
|
||
# Route each repainted cell to its layer: vents, ash and oasis
|
||
# water are ground; savanna, shrubland and pioneer growth are life.
|
||
changed = new_b != b_local
|
||
rows, cols = (a[changed] for a in np.nonzero(scatter_mask))
|
||
vals = new_b[changed]
|
||
ground = np.isin(vals, (19, 25, 2))
|
||
terrain[rows[ground], cols[ground]] = vals[ground]
|
||
biome[rows[~ground], cols[~ground]] = vals[~ground]
|
||
terrain[rows[~ground], cols[~ground]] = LAYER_NONE
|
||
|
||
# Vegetation ring around oasis lakes: dilate water cells from the
|
||
# scatter pass and assign graduated vegetation to the ring.
|
||
# water → coast vegetation → savanna/shrub → original biome
|
||
oasis_water = (composed() == 2) & land # scattered lake cells on land
|
||
if oasis_water.any():
|
||
from scipy.ndimage import binary_dilation
|
||
# Pixel-unit iteration counts scale by GRID_W/512 (module
|
||
# convention, see the GRID_W comment above) so the vegetation
|
||
# ring covers the same physical distance at any resolution — a
|
||
# fixed pixel count would shrink the ring (relative to the
|
||
# planet's real size) as GRID_W grows past the 512 baseline.
|
||
# LIVE BEHAVIOR CHANGE (T-964, PR #210 review — hoshe): GRID_W is
|
||
# ALREADY 1024 in production (not a future value), so this scales
|
||
# ring1/ring2 from the historical 2/4 iterations to 4/8 on EVERY
|
||
# planet generated today — see `oasis_ring_iterations`'s test
|
||
# coverage (`test_planet_simulation.py`) pinning both the 512
|
||
# (historical) and 1024 (current production) values explicitly.
|
||
ring1_iters, ring2_iters = oasis_ring_iterations(GRID_W)
|
||
ring1 = binary_dilation(oasis_water, iterations=ring1_iters) & ~oasis_water & land
|
||
ring2 = (
|
||
binary_dilation(oasis_water, iterations=ring2_iters)
|
||
& ~oasis_water
|
||
& ~ring1
|
||
& land
|
||
)
|
||
# Inner ring: lush vegetation (coast/lowland green)
|
||
paint_biome(ring1, 4) # lowland
|
||
# Outer ring: transitional (savanna/shrub)
|
||
paint_biome(ring2, 12) # shrubland
|
||
|
||
return biome, terrain
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Top-level simulate()
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def simulate(body_def: dict) -> dict:
|
||
"""
|
||
Run the full simulation stack for one body.
|
||
|
||
Parameters
|
||
----------
|
||
body_def : dict — from body_definition_parser.parse_system()
|
||
|
||
Returns
|
||
-------
|
||
dict terrain dict consumed by planet_renderer.render_globe()
|
||
Empty dict for gas giants (renderer handles those procedurally).
|
||
"""
|
||
planet_class = body_def.get("planet_class", "barren").replace("_ringed", "")
|
||
if planet_class == "gas_giant":
|
||
return {}
|
||
|
||
u, v, lat_frac, lon_frac, lat_rad = _make_grids()
|
||
|
||
elevation, sea_level, surface_water = compute_elevation(
|
||
body_def, u, v, lat_frac)
|
||
|
||
temperature, temp_clamped, temp_raw_K = compute_temperature(
|
||
body_def, elevation, sea_level, lat_frac)
|
||
|
||
moisture = compute_moisture(
|
||
body_def, elevation, sea_level, temperature, lat_frac, lon_frac)
|
||
|
||
hillshade = compute_hillshade(elevation)
|
||
|
||
biome_layer, terrain_layer = compute_biome_layers(
|
||
body_def, elevation, sea_level, surface_water, temperature, moisture)
|
||
biome = compose_layers(biome_layer, terrain_layer)
|
||
|
||
# Normalise temperature to [0,1] for renderer display — biome already computed
|
||
t_min, t_max = temperature.min(), temperature.max()
|
||
temperature_norm = ((temperature - t_min) / (t_max - t_min + 1e-9)).astype(np.float32)
|
||
|
||
return {
|
||
"elevation": elevation,
|
||
"temperature": temperature_norm, # normalised [0,1] for renderer
|
||
"moisture": moisture,
|
||
"biome": biome, # composed — what the renders colour
|
||
"biome_layer": biome_layer, # D-258: the baked biomemap (LAYER_NONE = none)
|
||
"terrain_layer": terrain_layer, # D-258: the baked terrainmap
|
||
"surface_water": surface_water,
|
||
"hillshade": hillshade,
|
||
"sea_level": sea_level,
|
||
"_grid_w": GRID_W,
|
||
"_grid_h": GRID_H,
|
||
# Audit trail
|
||
"temperature_clamped": temp_clamped,
|
||
"temperature_raw_K": round(temp_raw_K, 1),
|
||
"temperature_band_K": list(CLASS_T_BAND.get(
|
||
body_def.get("planet_class","").replace("_ringed",""), [None,None])),
|
||
}
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# CLI
|
||
# ---------------------------------------------------------------------------
|
||
|
||
if __name__ == "__main__":
|
||
import sys
|
||
import json
|
||
import time
|
||
import os
|
||
from PIL import Image
|
||
|
||
if len(sys.argv) < 2:
|
||
console.event("Usage: python -m tooling.domains.atlas.planet.planet_simulation body_def.json [--save-grids]")
|
||
sys.exit(1)
|
||
|
||
with open(sys.argv[1]) as f:
|
||
bd = json.load(f)
|
||
|
||
save_grids = "--save-grids" in sys.argv
|
||
|
||
console.event(f"Simulating: {bd['id']} ({bd['planet_class']})")
|
||
t0 = time.time()
|
||
terrain = simulate(bd)
|
||
|
||
if not terrain:
|
||
console.event("Gas giant — no terrain simulation.")
|
||
sys.exit(0)
|
||
|
||
dt = time.time() - t0
|
||
console.event(f"Done in {dt:.1f}s")
|
||
console.event(f" sea_level: {terrain['sea_level']:.3f}")
|
||
console.event(f" land cells: {(~terrain['surface_water']).sum()}")
|
||
|
||
ids, counts = np.unique(terrain['biome'], return_counts=True)
|
||
console.event(f" biomes: {list(zip(ids.tolist(), counts.tolist()))}")
|
||
|
||
if save_grids:
|
||
out = f"/tmp/{bd['id']}_grids"
|
||
os.makedirs(out, exist_ok=True)
|
||
for name in ("elevation", "temperature", "moisture", "hillshade"):
|
||
arr = terrain[name]
|
||
Image.fromarray((arr * 255).astype("uint8"), "L").save(
|
||
f"{out}/{name}.png")
|
||
console.event(f"Grids saved → {out}/")
|