feat(assets): add Sol system handcrafted terrain pipeline

Custom pipeline for GJ-0 (Sol) that imports real NASA/USGS
planetary data instead of procedural generation. Produces the
same output format (heightmap.png, globe.png, markers.json).

Real data bodies:
- Earth: ETOPO2022 elevation + WorldClim climate + 14 rivers
- Mars: MOLA DEM + ferric biome classes + terraformed water
- Luna: LOLA DEM + lunar biome palette

Procedural fallback for Mercury, Venus, Phobos, Deimos.
Synthetic elevation from albedo for Io, Europa, Ganymede,
Callisto, Titan, Enceladus. Gas giants use existing renderer.

New biome classes 34-36 (ferric_dust/highland/lowland) for
Mars iron oxide surface. Earth features: 50 cities (smart
scatter by continent), 15 named rivers, oceans, mountains.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

.gitignore

@@ -39,6 +39,8 @@ spikes/**/*.npz
 
 # Planet generator intermediates
 tooling/planet-gen/__pycache__/
+tooling/planet-gen/sol_data/.cache/
+tooling/planet-gen/sol_data/__pycache__/
 *.tmp.npz
 
 # Generated terrain grids (large, regenerated from pipeline)

tooling/planet-gen/biomes.toml

@@ -232,6 +232,11 @@ boreal                = [245, 278]   # cold forest / taiga
 29 = { name = "warm_dust",            cartographic = [210, 175, 120], photographic = [188, 155, 105] }
 30 = { name = "cold_rock",            cartographic = [160, 140, 115], photographic = [135, 118,  95] }
 
+# Ferric terrain (iron oxide — Mars, arid iron-rich worlds)
+34 = { name = "ferric_dust",          cartographic = [185, 110,  65], photographic = [158,  88,  48] }
+35 = { name = "ferric_highland",      cartographic = [165, 100,  60], photographic = [138,  78,  42] }
+36 = { name = "ferric_lowland",       cartographic = [200, 130,  75], photographic = [172, 105,  58] }
+
 # Lunar terrain (grey rock)
 31 = { name = "lunar_highland",       cartographic = [165, 165, 162], photographic = [138, 138, 135] }
 32 = { name = "lunar_mare",           cartographic = [120, 120, 118], photographic = [100, 100,  98] }

tooling/planet-gen/sol_data/__init__.py

@@ -0,0 +1 @@
+# Sol system real-world data importers

tooling/planet-gen/sol_data/download.py

@@ -0,0 +1,90 @@
+"""
+Caching downloader for planetary science datasets.
+
+Downloads are stored in sol_data/.cache/ and reused on subsequent runs.
+Supports resume for large files and optional SHA-256 verification.
+"""
+
+import hashlib
+import os
+import sys
+import urllib.request
+from pathlib import Path
+
+CACHE_DIR = Path(__file__).resolve().parent / ".cache"
+
+
+def _progress_hook(block_num, block_size, total_size):
+    """Print download progress."""
+    downloaded = block_num * block_size
+    if total_size > 0:
+        pct = min(100.0, downloaded * 100.0 / total_size)
+        mb = downloaded / (1024 * 1024)
+        total_mb = total_size / (1024 * 1024)
+        sys.stdout.write(f"\r    downloading: {mb:.1f}/{total_mb:.1f} MB ({pct:.0f}%)")
+    else:
+        mb = downloaded / (1024 * 1024)
+        sys.stdout.write(f"\r    downloading: {mb:.1f} MB")
+    sys.stdout.flush()
+
+
+def ensure_cached(url: str, filename: str, sha256: str = None) -> Path:
+    """
+    Download a file if not already cached. Returns path to cached file.
+
+    Parameters
+    ----------
+    url      : download URL
+    filename : local filename within the cache directory
+    sha256   : optional hex digest for verification
+    """
+    CACHE_DIR.mkdir(parents=True, exist_ok=True)
+    local_path = CACHE_DIR / filename
+
+    if local_path.exists():
+        if sha256:
+            actual = _sha256(local_path)
+            if actual != sha256:
+                print(f"    WARNING: checksum mismatch for {filename}, re-downloading")
+                local_path.unlink()
+            else:
+                return local_path
+        else:
+            return local_path
+
+    print(f"    fetching {filename} from {url[:80]}...")
+    tmp_path = local_path.with_suffix(".tmp")
+
+    try:
+        # Many government data servers (USGS, NOAA) require a User-Agent
+        opener = urllib.request.build_opener()
+        opener.addheaders = [
+            ("User-Agent", "SettledReach-PlanetGen/1.0 (terrain pipeline)"),
+        ]
+        urllib.request.install_opener(opener)
+        urllib.request.urlretrieve(url, str(tmp_path), reporthook=_progress_hook)
+        print()  # newline after progress
+    except Exception as e:
+        if tmp_path.exists():
+            tmp_path.unlink()
+        raise RuntimeError(f"Download failed for {filename}: {e}") from e
+
+    if sha256:
+        actual = _sha256(tmp_path)
+        if actual != sha256:
+            tmp_path.unlink()
+            raise RuntimeError(
+                f"Checksum mismatch for {filename}: "
+                f"expected {sha256[:16]}..., got {actual[:16]}..."
+            )
+
+    tmp_path.rename(local_path)
+    return local_path
+
+
+def _sha256(path: Path) -> str:
+    h = hashlib.sha256()
+    with open(path, "rb") as f:
+        for chunk in iter(lambda: f.read(8192), b""):
+            h.update(chunk)
+    return h.hexdigest()

tooling/planet-gen/sol_data/earth.py

@@ -0,0 +1,395 @@
+"""
+Earth (GJ0d) terrain builder.
+
+Data sources:
+  - Elevation: ETOPO 2022 60 arc-second (NOAA) — GeoTIFF
+  - Temperature: WorldClim v2.1 annual mean (10 arc-min) — GeoTIFF
+  - Precipitation: WorldClim v2.1 annual total (10 arc-min) — GeoTIFF
+  - Rivers: Natural Earth 10m rivers — GeoJSON
+
+All sources are equirectangular with col 0 = 180°W. ETOPO and WorldClim
+use col 0 = 180°W natively. Natural Earth uses -180 to 180 longitude.
+"""
+
+import json
+import math
+import os
+import struct
+import zipfile
+import numpy as np
+from pathlib import Path
+from scipy.ndimage import gaussian_filter
+
+from sol_data.download import ensure_cached
+from sol_data.shared import (
+    GRID_W, GRID_H,
+    load_tiff_as_array, load_raw_binary, load_image_as_elevation,
+    resample_to_grid, normalize_01, compute_sea_level,
+    greenwich_to_dateline, compute_hillshade, assemble_terrain,
+)
+
+# ─── Data source URLs ───────────────────────────────────────────────────────
+
+# ETOPO 2022 60 arc-second — surface elevation (ice surface, not bedrock)
+# ~130 MB GeoTIFF, 21600 x 10800, int16 metres
+ETOPO_URL = "https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO2022/data/60s/60s_surface_elev_gtif/ETOPO_2022_v1_60s_N90W180_surface.tif"
+ETOPO_FILE = "ETOPO_2022_v1_60s_N90W180_surface.tif"
+
+# WorldClim v2.1 — 10 arc-minute resolution (migrated to geodata.ucdavis.edu)
+# Temperature: mean annual, °C × 10 (int16), in a zip
+WCLIM_TEMP_URL = "https://geodata.ucdavis.edu/climate/worldclim/2_1/base/wc2.1_10m_tavg.zip"
+WCLIM_TEMP_FILE = "wc2.1_10m_tavg.zip"
+
+# Precipitation: annual total mm (int16), in a zip
+WCLIM_PREC_URL = "https://geodata.ucdavis.edu/climate/worldclim/2_1/base/wc2.1_10m_prec.zip"
+WCLIM_PREC_FILE = "wc2.1_10m_prec.zip"
+
+# Natural Earth 10m rivers — GeoJSON from GitHub
+RIVERS_URL = "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_10m_rivers_lake_centerlines.geojson"
+RIVERS_FILE = "ne_10m_rivers_lake_centerlines.geojson"
+
+# Earth physical constants
+EARTH_OCEAN_FRACTION = 0.71
+EARTH_MIN_ELEV_M = -10994.0  # Mariana Trench
+EARTH_MAX_ELEV_M = 8849.0    # Everest
+
+
+# ─── River filtering ────────────────────────────────────────────────────────
+
+# Rivers to include (smart scatter: 1-2 per continent + Rhine)
+INCLUDED_RIVERS = {
+    # Europe
+    "Danube", "Volga", "Rhine",
+    # North America
+    "Mississippi", "St. Lawrence",
+    # South America
+    "Amazon", "Paraná",
+    # Africa
+    "Nile", "Congo",
+    # West Asia
+    "Tigris",
+    # East/South Asia
+    "Yangtze", "Ganges", "Mekong",
+    # Australia
+    "Murray",
+}
+
+# Fuzzy matching — some NE names differ slightly
+RIVER_NAME_ALIASES = {
+    "Parana": "Paraná",
+    "Chang Jiang": "Yangtze",
+    "Huang He": "Yellow",
+    "Ganga": "Ganges",
+    "Nil": "Nile",
+    "Danau": "Danube",
+    "Donau": "Danube",
+    "Rhin": "Rhine",
+    "Rhein": "Rhine",
+    "Saint Lawrence": "St. Lawrence",
+    "St Lawrence": "St. Lawrence",
+    "Río Paraná": "Paraná",
+    "Rio Parana": "Paraná",
+}
+
+
+def _match_river_name(feature_name: str) -> str:
+    """Check if a Natural Earth river name matches our included set."""
+    if not feature_name:
+        return None
+    name = feature_name.strip()
+    # Direct match
+    if name in INCLUDED_RIVERS:
+        return name
+    # Alias match
+    if name in RIVER_NAME_ALIASES:
+        alias = RIVER_NAME_ALIASES[name]
+        if alias in INCLUDED_RIVERS:
+            return alias
+    # Substring match (e.g. "Mississippi River" contains "Mississippi")
+    for included in INCLUDED_RIVERS:
+        if included.lower() in name.lower() or name.lower() in included.lower():
+            return included
+    return None
+
+
+# ─── Data loaders ───────────────────────────────────────────────────────────
+
+def _load_etopo() -> np.ndarray:
+    """Load ETOPO 2022 elevation data, return raw metres array."""
+    path = ensure_cached(ETOPO_URL, ETOPO_FILE)
+    print(f"    loading ETOPO: {path}")
+    try:
+        arr = load_tiff_as_array(str(path))
+    except Exception as e:
+        raise RuntimeError(
+            f"Failed to load ETOPO GeoTIFF: {e}\n"
+            f"If PIL can't read this TIFF, install Pillow with TIFF support "
+            f"or convert to raw binary."
+        ) from e
+    print(f"    ETOPO shape: {arr.shape}, range: [{arr.min():.0f}, {arr.max():.0f}] m")
+    return arr
+
+
+def _load_worldclim_temperature() -> np.ndarray:
+    """
+    Load WorldClim v2.1 annual mean temperature.
+    Returns temperature in Kelvin at native resolution.
+    """
+    zip_path = ensure_cached(WCLIM_TEMP_URL, WCLIM_TEMP_FILE)
+    print(f"    loading WorldClim temperature: {zip_path}")
+
+    # The zip contains monthly TIFFs (tavg_01.tif to tavg_12.tif).
+    # Compute annual mean from all 12 months.
+    cache_dir = zip_path.parent
+    monthly_sum = None
+    count = 0
+
+    with zipfile.ZipFile(zip_path) as zf:
+        tif_names = sorted([n for n in zf.namelist() if n.endswith(".tif")])
+        for tif_name in tif_names:
+            extracted = cache_dir / Path(tif_name).name
+            if not extracted.exists():
+                zf.extract(tif_name, cache_dir)
+                # Handle nested paths in zip
+                nested = cache_dir / tif_name
+                if nested != extracted and nested.exists():
+                    nested.rename(extracted)
+            try:
+                arr = load_tiff_as_array(str(extracted))
+            except Exception:
+                # Try the nested path
+                nested = cache_dir / tif_name
+                if nested.exists():
+                    arr = load_tiff_as_array(str(nested))
+                else:
+                    continue
+            # Replace nodata with NaN
+            arr[arr < -999] = np.nan
+            if monthly_sum is None:
+                monthly_sum = arr.copy()
+            else:
+                monthly_sum += arr
+            count += 1
+
+    if count == 0:
+        raise RuntimeError("No temperature TIFFs found in WorldClim archive")
+
+    # Annual mean (WorldClim tavg is °C × 10)
+    temp_C = (monthly_sum / count) / 10.0
+    # Convert to Kelvin
+    temp_K = temp_C + 273.15
+    # Replace NaN (ocean/nodata) with a reasonable ocean temperature
+    temp_K = np.nan_to_num(temp_K, nan=288.0)
+
+    print(f"    WorldClim temp shape: {temp_K.shape}, "
+          f"range: [{np.nanmin(temp_K):.0f}, {np.nanmax(temp_K):.0f}] K")
+    return temp_K
+
+
+def _load_worldclim_precipitation() -> np.ndarray:
+    """
+    Load WorldClim v2.1 annual precipitation (sum of 12 months).
+    Returns precipitation in mm/year at native resolution.
+    """
+    zip_path = ensure_cached(WCLIM_PREC_URL, WCLIM_PREC_FILE)
+    print(f"    loading WorldClim precipitation: {zip_path}")
+
+    cache_dir = zip_path.parent
+    annual_sum = None
+
+    with zipfile.ZipFile(zip_path) as zf:
+        tif_names = sorted([n for n in zf.namelist() if n.endswith(".tif")])
+        for tif_name in tif_names:
+            extracted = cache_dir / Path(tif_name).name
+            if not extracted.exists():
+                zf.extract(tif_name, cache_dir)
+                nested = cache_dir / tif_name
+                if nested != extracted and nested.exists():
+                    nested.rename(extracted)
+            try:
+                arr = load_tiff_as_array(str(extracted))
+            except Exception:
+                nested = cache_dir / tif_name
+                if nested.exists():
+                    arr = load_tiff_as_array(str(nested))
+                else:
+                    continue
+            arr[arr < -999] = 0.0
+            if annual_sum is None:
+                annual_sum = arr.copy()
+            else:
+                annual_sum += arr
+
+    if annual_sum is None:
+        raise RuntimeError("No precipitation TIFFs found in WorldClim archive")
+
+    print(f"    WorldClim precip shape: {annual_sum.shape}, "
+          f"range: [{annual_sum.min():.0f}, {annual_sum.max():.0f}] mm/yr")
+    return annual_sum
+
+
+def _load_rivers_geojson() -> list:
+    """
+    Load Natural Earth rivers GeoJSON and extract polylines for included rivers.
+    Returns list of (name, [(row, col), ...]) in grid coordinates.
+    """
+    path = ensure_cached(RIVERS_URL, RIVERS_FILE)
+    print(f"    loading rivers: {path}")
+
+    with open(path) as f:
+        geojson = json.load(f)
+
+    rivers = []
+    for feature in geojson.get("features", []):
+        props = feature.get("properties", {})
+        fname = props.get("name") or props.get("name_en") or ""
+        matched = _match_river_name(fname)
+        if not matched:
+            continue
+
+        geom = feature.get("geometry", {})
+        geom_type = geom.get("type", "")
+        coords_list = []
+
+        if geom_type == "LineString":
+            coords_list = [geom["coordinates"]]
+        elif geom_type == "MultiLineString":
+            coords_list = geom["coordinates"]
+        else:
+            continue
+
+        for coords in coords_list:
+            path_grid = []
+            for lon, lat in coords:
+                # Convert lon/lat to grid coordinates
+                # Grid: row 0 = 90°N, row 255 = 90°S
+                #        col 0 = 180°W, col 511 = 180°E
+                row = int((90.0 - lat) / 180.0 * GRID_H)
+                col = int((lon + 180.0) / 360.0 * GRID_W)
+                row = max(0, min(GRID_H - 1, row))
+                col = max(0, min(GRID_W - 1, col))
+                # Deduplicate: skip if same grid cell as previous point.
+                # Natural Earth has hundreds of lon/lat points per river,
+                # many of which land on the same 512x256 cell. Without
+                # dedup, the renderer sees len(path)=300 and draws width 6.
+                if path_grid and path_grid[-1] == (row, col):
+                    continue
+                path_grid.append((row, col))
+            if len(path_grid) >= 2:
+                rivers.append((matched, path_grid))
+
+    # Deduplicate: keep longest segment per river name
+    by_name = {}
+    for name, path in rivers:
+        if name not in by_name or len(path) > len(by_name[name]):
+            by_name[name] = path
+
+    print(f"    matched {len(by_name)} rivers: {', '.join(sorted(by_name.keys()))}")
+    return [(name, path) for name, path in by_name.items()]
+
+
+# ─── Main builder ───────────────────────────────────────────────────────────
+
+def build_terrain(body_def: dict) -> dict:
+    """
+    Build Earth terrain dict from real-world data.
+
+    Returns the same dict format as planet_simulation.simulate().
+    """
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+    from planet_simulation import compute_biome
+
+    print("  Earth: loading real-world data...")
+
+    # ── 1. Elevation ────────────────────────────────────────────────────
+    etopo_raw = _load_etopo()
+
+    # ETOPO 2022 N90W180 is already col 0 = 180°W — no shift needed
+    # Resample to grid
+    elevation_m = resample_to_grid(etopo_raw, GRID_H, GRID_W, order=1)
+
+    # Normalise to [0, 1]
+    elevation = normalize_01(elevation_m, EARTH_MIN_ELEV_M, EARTH_MAX_ELEV_M)
+
+    # Sea level: Earth's ocean fraction is ~0.71
+    sea_level = compute_sea_level(elevation, EARTH_OCEAN_FRACTION)
+    surface_water = elevation < sea_level
+
+    print(f"    elevation: sea_level={sea_level:.4f}, "
+          f"ocean={surface_water.sum()}/{GRID_H*GRID_W} cells")
+
+    # ── 2. Temperature ──────────────────────────────────────────────────
+    temp_raw_K = _load_worldclim_temperature()
+
+    # WorldClim uses col 0 = 180°W — no shift needed
+    temperature_K = resample_to_grid(temp_raw_K, GRID_H, GRID_W, order=1)
+
+    # Fill ocean areas with latitude-dependent ocean temperature
+    v = np.linspace(0, 1, GRID_H, dtype=np.float32)
+    lat_abs = np.abs(v - 0.5) * 2.0  # 0 at equator, 1 at poles
+    ocean_temp = 301.0 - lat_abs[:, np.newaxis] * 30.0  # ~28°C equator, ~-2°C poles
+    temperature_K = np.where(surface_water, ocean_temp, temperature_K)
+
+    print(f"    temperature: [{temperature_K.min():.0f}, {temperature_K.max():.0f}] K")
+
+    # ── 3. Moisture ─────────────────────────────────────────────────────
+    precip_raw = _load_worldclim_precipitation()
+
+    # WorldClim uses col 0 = 180°W — no shift needed
+    precip = resample_to_grid(precip_raw, GRID_H, GRID_W, order=1)
+
+    # Normalise to [0, 1] — global max is ~10000 mm/yr (tropical rainforest)
+    moisture = normalize_01(precip, 0.0, 6000.0)
+    # Ocean moisture = high (drives adjacent land humidity)
+    moisture = np.where(surface_water, 0.9, moisture)
+
+    print(f"    moisture: [{moisture.min():.2f}, {moisture.max():.2f}]")
+
+    # ── 4. Biome classification ─────────────────────────────────────────
+    # Use the existing Whittaker table with real temperature and moisture
+    biome = compute_biome(body_def, elevation, sea_level, surface_water,
+                          temperature_K, moisture)
+    n_biomes = len(np.unique(biome))
+    print(f"    biomes: {n_biomes} classes present")
+
+    # ── 5. Hillshade ────────────────────────────────────────────────────
+    hillshade = compute_hillshade(elevation)
+
+    # ── 6. Rivers ───────────────────────────────────────────────────────
+    named_rivers = _load_rivers_geojson()
+
+    # Clip rivers: stop each path when it hits surface water.
+    # Rivers like the Amazon/Nile/Rhine otherwise draw through seas.
+    clipped = []
+    for name, path in named_rivers:
+        clipped_path = []
+        for r, c in path:
+            if surface_water[r, c]:
+                break
+            clipped_path.append((r, c))
+        if len(clipped_path) >= 2:
+            clipped.append((name, clipped_path))
+
+    n_orig = len(named_rivers)
+    n_kept = len(clipped)
+    print(f"    rivers: {n_kept}/{n_orig} kept after water clipping")
+    named_rivers = clipped
+    rivers = [path for _, path in named_rivers]
+
+    # ── 7. Assemble ─────────────────────────────────────────────────────
+    terrain = assemble_terrain(
+        elevation=elevation,
+        temperature_K=temperature_K,
+        moisture=moisture,
+        biome=biome,
+        surface_water=surface_water,
+        hillshade=hillshade,
+        rivers=rivers,
+        sea_level=sea_level,
+    )
+
+    # Store river names for the marker overlay
+    terrain["_river_names"] = {i: name for i, (name, _) in enumerate(named_rivers)}
+
+    return terrain

tooling/planet-gen/sol_data/gas_giants.py

@@ -0,0 +1,25 @@
+"""
+Gas giant body definition helpers for Jupiter, Saturn, Uranus, Neptune.
+
+Gas giants have no solid surface — the existing planet_renderer._render_gas_giant()
+handles band patterns procedurally. This module only provides configuration
+validation and body_def enhancement. No terrain dict is produced.
+
+The actual overrides are in sol_overrides.json and applied by the body
+definition parser. This module exists for future enhancement (ring tuning,
+storm placement, etc).
+"""
+
+
+def validate_gas_giant_def(body_def: dict) -> bool:
+    """Check that a gas giant body_def has required fields for rendering."""
+    pc = body_def.get("planet_class", "")
+    if "gas_giant" not in pc and pc not in ("gas_giant",):
+        return False
+
+    gg = body_def.get("gas_giant", {})
+    if not gg.get("band_palette"):
+        print(f"    WARNING: {body_def['id']} missing gas_giant.band_palette")
+        return False
+
+    return True

tooling/planet-gen/sol_data/ice_moons.py

@@ -0,0 +1,143 @@
+"""
+Ice moon terrain builder — Europa, Ganymede, Callisto, Enceladus.
+
+These bodies lack high-quality global DEMs. We use available mosaics
+(albedo/reflectance) to derive synthetic elevation:
+  - Bright = ice ridges/highlands (high)
+  - Dark = mare/chaos terrain/craters (low)
+
+Each moon gets specific temperature and appearance tuning.
+"""
+
+import numpy as np
+from pathlib import Path
+from scipy.ndimage import gaussian_filter
+
+from sol_data.download import ensure_cached
+from sol_data.shared import (
+    GRID_W, GRID_H,
+    load_image_as_elevation, resample_to_grid, normalize_01,
+    compute_hillshade, assemble_terrain,
+    temperature_grid_analytical,
+)
+
+# ─── Per-moon configuration ─────────────────────────────────────────────────
+
+MOON_CONFIG = {
+    "GJ0f-2": {  # Europa
+        "name": "Europa",
+        "mosaic_url": "https://astrogeology.usgs.gov/cache/images/3c79b3867c0dc5ec2ea33e485a079e58_europa_voyager_galileo_ssi_global_mosaic_500m.jpg",
+        "mosaic_file": "europa_galileo_mosaic.jpg",
+        "base_temp_K": 102.0,
+        "lat_gradient_K": 10.0,
+        "sigma": 2.0,         # smooth albedo → elevation
+        "invert_albedo": False,  # bright = ridges (high)
+    },
+    "GJ0f-3": {  # Ganymede
+        "name": "Ganymede",
+        "mosaic_url": "https://astrogeology.usgs.gov/cache/images/f60b3c06c92f59834f2d4cf9b46cb8f7_ganymede_voyager_galileo_global_mosaic_1km.jpg",
+        "mosaic_file": "ganymede_galileo_mosaic.jpg",
+        "base_temp_K": 110.0,
+        "lat_gradient_K": 15.0,
+        "sigma": 3.0,
+        "invert_albedo": False,
+    },
+    "GJ0f-4": {  # Callisto
+        "name": "Callisto",
+        "mosaic_url": "https://astrogeology.usgs.gov/cache/images/26b4e80eeb35d46c53d56cded56deeef_callisto_voyager_galileo_global_mosaic_1km.jpg",
+        "mosaic_file": "callisto_galileo_mosaic.jpg",
+        "base_temp_K": 115.0,
+        "lat_gradient_K": 12.0,
+        "sigma": 4.0,
+        "invert_albedo": False,
+    },
+    "GJ0g-2": {  # Enceladus
+        "name": "Enceladus",
+        "mosaic_url": "https://astrogeology.usgs.gov/cache/images/1e9fede316c8c47fdc0b96f4c09e4915_enceladus_cassini_iss_global_mosaic_100m.jpg",
+        "mosaic_file": "enceladus_cassini_mosaic.jpg",
+        "base_temp_K": 75.0,
+        "lat_gradient_K": 8.0,
+        "sigma": 2.0,
+        "invert_albedo": False,
+    },
+}
+
+
+def _load_mosaic_as_elevation(config: dict) -> np.ndarray:
+    """Load a global mosaic and convert to synthetic elevation."""
+    try:
+        path = ensure_cached(config["mosaic_url"], config["mosaic_file"])
+        print(f"    loading {config['name']} mosaic: {path}")
+        albedo = load_image_as_elevation(str(path),
+                                          invert=config.get("invert_albedo", False))
+        albedo = resample_to_grid(albedo, GRID_H, GRID_W, order=1)
+    except Exception as e:
+        print(f"    WARNING: {config['name']} mosaic unavailable ({e}), synthetic")
+        albedo = _synthetic_ice_terrain(config["name"])
+
+    # Smooth albedo to create plausible topography
+    sigma = config.get("sigma", 3.0)
+    elevation = gaussian_filter(albedo, sigma=sigma)
+    return normalize_01(elevation)
+
+
+def _synthetic_ice_terrain(name: str) -> np.ndarray:
+    """Generate synthetic ice moon terrain if mosaic unavailable."""
+    seed = hash(name) & 0xFFFFFFFF
+    rng = np.random.default_rng(seed)
+    base = rng.random((GRID_H, GRID_W)).astype(np.float32)
+    base = gaussian_filter(base, sigma=6.0)
+    # Add craters
+    for _ in range(20):
+        cy, cx = rng.integers(0, GRID_H), rng.integers(0, GRID_W)
+        r = rng.integers(5, 20)
+        y, x = np.ogrid[-cy:GRID_H-cy, -cx:GRID_W-cx]
+        mask = x*x + y*y <= r*r
+        base[mask] *= 0.5
+    return normalize_01(base)
+
+
+def build_terrain(body_def: dict) -> dict:
+    """Build ice moon terrain dict from mosaic data."""
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+    from planet_simulation import compute_biome
+
+    body_id = body_def["id"]
+    config = MOON_CONFIG.get(body_id)
+
+    if config is None:
+        raise ValueError(f"No ice moon config for {body_id}")
+
+    print(f"  {config['name']}: loading data...")
+
+    # ── 1. Elevation ────────────────────────────────────────────────────
+    elevation = _load_mosaic_as_elevation(config)
+    sea_level = 0.0
+    surface_water = np.zeros((GRID_H, GRID_W), dtype=bool)
+
+    # ── 2. Temperature ──────────────────────────────────────────────────
+    temperature_K = temperature_grid_analytical(
+        base_T_K=config["base_temp_K"],
+        elevation=elevation,
+        lapse_rate_K_per_unit=5.0,
+        lat_gradient_K=config["lat_gradient_K"],
+    )
+    temperature_K = np.maximum(temperature_K, 40.0)
+
+    # ── 3. Moisture ─────────────────────────────────────────────────────
+    moisture = np.zeros((GRID_H, GRID_W), dtype=np.float32)
+
+    # ── 4. Biome ────────────────────────────────────────────────────────
+    biome = compute_biome(body_def, elevation, sea_level, surface_water,
+                          temperature_K, moisture)
+
+    # ── 5. Hillshade ────────────────────────────────────────────────────
+    hillshade = compute_hillshade(elevation)
+
+    return assemble_terrain(
+        elevation=elevation, temperature_K=temperature_K,
+        moisture=moisture, biome=biome,
+        surface_water=surface_water, hillshade=hillshade,
+        rivers=[], sea_level=sea_level,
+    )

tooling/planet-gen/sol_data/io_moon.py

@@ -0,0 +1,117 @@
+"""
+Io (GJ0f-1) terrain builder.
+
+Io is the most volcanically active body in the solar system due to
+tidal heating from Jupiter. Surface is covered in sulfur and volcanic
+deposits. No published global DEM exists at useful resolution — we use
+the Galileo/Voyager global mosaic (albedo) to derive synthetic elevation.
+
+Data source:
+  - Surface: USGS Io Galileo/Voyager global mosaic
+  - Elevation: synthetic from albedo (dark = caldera/lava, bright = sulfur)
+
+Properties:
+  - Surface temp: ~130K background, 400-1800K at volcanic hotspots
+  - planet_class: "volcanic", atmosphere: "none"
+"""
+
+import numpy as np
+from pathlib import Path
+from scipy.ndimage import gaussian_filter
+
+from sol_data.download import ensure_cached
+from sol_data.shared import (
+    GRID_W, GRID_H,
+    load_image_as_elevation, resample_to_grid, normalize_01,
+    compute_hillshade, assemble_terrain,
+    temperature_grid_analytical,
+)
+
+# Io global mosaic (Galileo SSI + Voyager) — JPEG from USGS
+# If direct download isn't available, fall back to procedural
+IO_MOSAIC_URL = "https://astrogeology.usgs.gov/cache/images/bf08a5b6fa0c2ed73117dc1b6c516fa8_io_galileo_voyager_global_mosaic_1km.jpg"
+IO_MOSAIC_FILE = "io_galileo_mosaic.jpg"
+
+IO_BACKGROUND_TEMP_K = 130.0
+IO_HOTSPOT_TEMP_K = 600.0
+
+
+def _load_io_mosaic() -> np.ndarray:
+    """Load Io global mosaic and convert to synthetic elevation."""
+    try:
+        path = ensure_cached(IO_MOSAIC_URL, IO_MOSAIC_FILE)
+        print(f"    loading Io mosaic: {path}")
+        albedo = load_image_as_elevation(str(path), invert=False)
+    except Exception as e:
+        print(f"    WARNING: Io mosaic unavailable ({e}), generating synthetic")
+        return _synthetic_io_terrain()
+
+    # Resample to grid
+    albedo = resample_to_grid(albedo, GRID_H, GRID_W, order=1)
+
+    # Convert albedo to elevation:
+    # Dark regions (low albedo) = calderas/lava flows = low elevation
+    # Bright regions (high albedo) = sulfur deposits = high elevation
+    # Smooth to create plausible topography
+    elevation = gaussian_filter(albedo, sigma=3.0)
+    elevation = normalize_01(elevation)
+
+    return elevation
+
+
+def _synthetic_io_terrain() -> np.ndarray:
+    """Generate synthetic Io-like terrain if mosaic unavailable."""
+    rng = np.random.default_rng(42)
+    base = rng.random((GRID_H, GRID_W)).astype(np.float32)
+    base = gaussian_filter(base, sigma=8.0)
+    # Add volcanic calderas (circular depressions)
+    for _ in range(30):
+        cy, cx = rng.integers(0, GRID_H), rng.integers(0, GRID_W)
+        r = rng.integers(3, 15)
+        y, x = np.ogrid[-cy:GRID_H-cy, -cx:GRID_W-cx]
+        mask = x*x + y*y <= r*r
+        base[mask] *= 0.3
+    return normalize_01(base)
+
+
+def build_terrain(body_def: dict) -> dict:
+    """Build Io terrain dict."""
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+    from planet_simulation import compute_biome
+
+    print("  Io: loading data...")
+
+    # ── 1. Elevation ────────────────────────────────────────────────────
+    elevation = _load_io_mosaic()
+    sea_level = 0.0
+    surface_water = np.zeros((GRID_H, GRID_W), dtype=bool)
+
+    # ── 2. Temperature ──────────────────────────────────────────────────
+    # Background ~130K, volcanic hotspots much hotter
+    temperature_K = temperature_grid_analytical(
+        base_T_K=IO_BACKGROUND_TEMP_K,
+        elevation=elevation,
+        lapse_rate_K_per_unit=-200.0,  # low elevation = hot (lava)
+        lat_gradient_K=10.0,
+    )
+    # Volcanic hotspots: low-elevation areas are hot
+    hotspot_mask = elevation < 0.25
+    temperature_K[hotspot_mask] += 300.0
+
+    # ── 3. Moisture ─────────────────────────────────────────────────────
+    moisture = np.zeros((GRID_H, GRID_W), dtype=np.float32)
+
+    # ── 4. Biome ────────────────────────────────────────────────────────
+    biome = compute_biome(body_def, elevation, sea_level, surface_water,
+                          temperature_K, moisture)
+
+    # ── 5. Hillshade ────────────────────────────────────────────────────
+    hillshade = compute_hillshade(elevation)
+
+    return assemble_terrain(
+        elevation=elevation, temperature_K=temperature_K,
+        moisture=moisture, biome=biome,
+        surface_water=surface_water, hillshade=hillshade,
+        rivers=[], sea_level=sea_level,
+    )

tooling/planet-gen/sol_data/luna.py

@@ -0,0 +1,134 @@
+"""
+Luna (GJ0d-1) terrain builder.
+
+Data source:
+  - Elevation: LOLA (Lunar Orbiter Laser Altimeter) DEM
+    Available at various resolutions from USGS Astrogeology.
+    We use the 4ppd (1440×720) or 16ppd version.
+
+Luna properties:
+  - Min elevation: ~-9100 m (South Pole-Aitken basin)
+  - Max elevation: ~10786 m (near Engel'gardt crater rim)
+  - No atmosphere, no water
+  - body_type: "moon" → uses lunar biome palette (classes 31/32/33)
+"""
+
+import numpy as np
+from pathlib import Path
+
+from sol_data.download import ensure_cached
+from sol_data.shared import (
+    GRID_W, GRID_H,
+    load_raw_binary, load_tiff_as_array, load_image_as_elevation,
+    resample_to_grid, normalize_01,
+    compute_hillshade, assemble_terrain,
+    temperature_grid_analytical,
+)
+
+# LOLA GDR — available as PDS IMG files
+# 4ppd (1440 × 720) — compact version
+LOLA_4PPD_URL = "https://pds-geosciences.wustl.edu/lro/lro-l-lola-3-rdr-v1/lrolol_1xxx/data/lola_gdr/cylindrical/img/ldem_4.img"
+LOLA_4PPD_FILE = "lola_gdr_4ppd.img"
+LOLA_4PPD_W = 1440
+LOLA_4PPD_H = 720
+
+# 16ppd (5760 × 2880) — higher quality
+LOLA_16PPD_URL = "https://pds-geosciences.wustl.edu/lro/lro-l-lola-3-rdr-v1/lrolol_1xxx/data/lola_gdr/cylindrical/img/ldem_16.img"
+LOLA_16PPD_FILE = "lola_gdr_16ppd.img"
+LOLA_16PPD_W = 5760
+LOLA_16PPD_H = 2880
+
+# Luna physical constants
+LUNA_MIN_ELEV_M = -9100.0
+LUNA_MAX_ELEV_M = 10786.0
+LUNA_EQUATORIAL_TEMP_K = 220.0   # mean dayside ~220K
+LUNA_POLAR_TEMP_K = 100.0        # permanently shadowed craters ~40K, average ~100K
+
+
+def _load_lola(use_16ppd: bool = False) -> np.ndarray:
+    """Load LOLA DEM, return elevation in metres."""
+    if use_16ppd:
+        url, filename, w, h = LOLA_16PPD_URL, LOLA_16PPD_FILE, LOLA_16PPD_W, LOLA_16PPD_H
+    else:
+        url, filename, w, h = LOLA_4PPD_URL, LOLA_4PPD_FILE, LOLA_4PPD_W, LOLA_4PPD_H
+
+    path = ensure_cached(url, filename)
+    print(f"    loading LOLA: {path} ({w}x{h})")
+
+    # LOLA GDR: little-endian int16 (LSB_INTEGER per PDS label)
+    # with a scaling factor of 0.5 metres.
+    try:
+        arr = load_raw_binary(str(path), w, h, dtype="<i2", offset=0)
+        # LOLA int16 values are in units of 0.5m (scale factor 0.5)
+        arr = arr * 0.5
+    except ValueError:
+        # If int16 doesn't work, try float32
+        arr = load_raw_binary(str(path), w, h, dtype="<f4", offset=0)
+
+    # Handle nodata
+    arr[arr > 20000] = 0.0
+    arr[arr < -20000] = 0.0
+
+    print(f"    LOLA range: [{arr.min():.0f}, {arr.max():.0f}] m")
+    return arr
+
+
+def build_terrain(body_def: dict) -> dict:
+    """Build Luna terrain dict from LOLA data."""
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+    from planet_simulation import compute_biome
+
+    print("  Luna: loading LOLA data...")
+
+    # ── 1. Elevation ────────────────────────────────────────────────────
+    lola_raw = _load_lola(use_16ppd=False)
+
+    # LOLA cylindrical: col 0 = 0° longitude — shift to 180°W
+    from sol_data.shared import greenwich_to_dateline
+    lola_shifted = greenwich_to_dateline(lola_raw)
+
+    elevation_m = resample_to_grid(lola_shifted, GRID_H, GRID_W, order=1)
+    elevation = normalize_01(elevation_m, LUNA_MIN_ELEV_M, LUNA_MAX_ELEV_M)
+
+    # No liquid — sea level at 0
+    sea_level = 0.0
+    surface_water = np.zeros((GRID_H, GRID_W), dtype=bool)
+
+    print(f"    elevation normalised")
+
+    # ── 2. Temperature ──────────────────────────────────────────────────
+    temperature_K = temperature_grid_analytical(
+        base_T_K=LUNA_EQUATORIAL_TEMP_K,
+        elevation=elevation,
+        lapse_rate_K_per_unit=10.0,
+        lat_gradient_K=120.0,  # huge contrast equator to poles
+    )
+    # Clamp minimum
+    temperature_K = np.maximum(temperature_K, 40.0)
+
+    print(f"    temperature: [{temperature_K.min():.0f}, {temperature_K.max():.0f}] K")
+
+    # ── 3. Moisture ─────────────────────────────────────────────────────
+    moisture = np.zeros((GRID_H, GRID_W), dtype=np.float32)
+
+    # ── 4. Biome ────────────────────────────────────────────────────────
+    # body_type: "moon" + atmosphere: "none" → lunar palette (31/32/33)
+    biome = compute_biome(body_def, elevation, sea_level, surface_water,
+                          temperature_K, moisture)
+    print(f"    biomes: {len(np.unique(biome))} classes")
+
+    # ── 5. Hillshade ────────────────────────────────────────────────────
+    hillshade = compute_hillshade(elevation)
+
+    # ── 6. Assemble ─────────────────────────────────────────────────────
+    return assemble_terrain(
+        elevation=elevation,
+        temperature_K=temperature_K,
+        moisture=moisture,
+        biome=biome,
+        surface_water=surface_water,
+        hillshade=hillshade,
+        rivers=[],
+        sea_level=sea_level,
+    )

tooling/planet-gen/sol_data/mars.py

@@ -0,0 +1,190 @@
+"""
+Mars (GJ0e) terrain builder.
+
+Data source:
+  - Elevation: MOLA MEGDR (Mars Orbiter Laser Altimeter)
+    PDS format, big-endian int16, metres relative to areoid.
+    Available at multiple resolutions. We use 4ppd (1440×720)
+    or 16ppd (5760×2880) — both small enough to download quickly.
+
+Mars properties:
+  - Min elevation: ~-8200 m (Hellas Basin)
+  - Max elevation: ~21229 m (Olympus Mons)
+  - Polar ice caps: CO2 + water ice
+  - Thin atmosphere (6 mbar) — classified as "thin" in body_def
+  - Almost no liquid water (hydrosphere: "ice")
+"""
+
+import numpy as np
+from pathlib import Path
+
+from sol_data.download import ensure_cached
+from sol_data.shared import (
+    GRID_W, GRID_H,
+    load_raw_binary, resample_to_grid, normalize_01, compute_sea_level,
+    compute_hillshade, assemble_terrain,
+    temperature_grid_analytical, temperature_equilibrium_K,
+)
+
+# MOLA MEGDR — 4 pixels per degree (1440 × 720), big-endian int16
+# Each pixel = metres relative to Mars areoid
+# PDS binary with no header (data starts at byte 0 for .img files)
+MOLA_4PPD_URL = "https://pds-geosciences.wustl.edu/mgs/mgs-m-mola-5-megdr-l3-v1/mgsl_300x/meg004/megt90n000cb.img"
+MOLA_4PPD_FILE = "mola_megdr_4ppd.img"
+MOLA_4PPD_W = 1440
+MOLA_4PPD_H = 720
+
+# Alternative: 16ppd (5760 × 2880) for higher quality
+MOLA_16PPD_URL = "https://pds-geosciences.wustl.edu/mgs/mgs-m-mola-5-megdr-l3-v1/mgsl_300x/meg016/megt90n000eb.img"
+MOLA_16PPD_FILE = "mola_megdr_16ppd.img"
+MOLA_16PPD_W = 5760
+MOLA_16PPD_H = 2880
+
+# Mars physical constants
+MARS_MIN_ELEV_M = -8200.0   # Hellas Basin
+MARS_MAX_ELEV_M = 21229.0   # Olympus Mons summit
+# Real Mars temperatures — we don't fudge these. Mars colour comes from
+# ferric biome classes (34/35/36) applied based on iron oxide substrate.
+MARS_EQUATORIAL_TEMP_K = 215.0  # daytime average near equator
+MARS_POLAR_TEMP_K = 150.0
+MARS_OCEAN_FRACTION = 0.0    # no liquid water (ice only)
+
+# Ferric biome class IDs (from biomes.toml)
+FERRIC_DUST = 34
+FERRIC_HIGHLAND = 35
+FERRIC_LOWLAND = 36
+
+
+def _load_mola(use_16ppd: bool = False) -> np.ndarray:
+    """Load MOLA DEM, return elevation in metres."""
+    if use_16ppd:
+        url, filename, w, h = MOLA_16PPD_URL, MOLA_16PPD_FILE, MOLA_16PPD_W, MOLA_16PPD_H
+    else:
+        url, filename, w, h = MOLA_4PPD_URL, MOLA_4PPD_FILE, MOLA_4PPD_W, MOLA_4PPD_H
+
+    path = ensure_cached(url, filename)
+    print(f"    loading MOLA: {path} ({w}x{h})")
+
+    # MOLA MEGDR: big-endian int16, metres, no header
+    arr = load_raw_binary(str(path), w, h, dtype=">i2", offset=0)
+
+    # MOLA nodata is typically 32767 or -32768
+    arr[arr > 30000] = 0.0
+    arr[arr < -30000] = 0.0
+
+    print(f"    MOLA range: [{arr.min():.0f}, {arr.max():.0f}] m")
+    return arr
+
+
+def build_terrain(body_def: dict) -> dict:
+    """Build Mars terrain dict from MOLA data."""
+    print("  Mars: loading MOLA data...")
+
+    # ── 1. Elevation ────────────────────────────────────────────────────
+    mola_raw = _load_mola(use_16ppd=False)
+
+    # MOLA is col 0 = 0° longitude — shift to col 0 = 180°W
+    from sol_data.shared import greenwich_to_dateline
+    mola_shifted = greenwich_to_dateline(mola_raw)
+
+    # Resample to grid
+    elevation_m = resample_to_grid(mola_shifted, GRID_H, GRID_W, order=1)
+
+    # Normalise to [0, 1]
+    elevation = normalize_01(elevation_m, MARS_MIN_ELEV_M, MARS_MAX_ELEV_M)
+
+    print(f"    elevation normalised")
+
+    # ── 2. Temperature ──────────────────────────────────────────────────
+    # Analytical: equatorial ~210K, polar ~150K, elevation lapse
+    temperature_K = temperature_grid_analytical(
+        base_T_K=MARS_EQUATORIAL_TEMP_K,
+        elevation=elevation,
+        lapse_rate_K_per_unit=30.0,
+        lat_gradient_K=60.0,
+    )
+
+    # Polar ice caps: very cold at high latitudes
+    v = np.linspace(0, 1, GRID_H, dtype=np.float32)
+    lat_abs = np.abs(v - 0.5) * 2.0
+    polar_rows = lat_abs > 0.75
+    temperature_K[polar_rows, :] = np.minimum(temperature_K[polar_rows, :], 155.0)
+
+    print(f"    temperature: [{temperature_K.min():.0f}, {temperature_K.max():.0f}] K")
+
+    # ── 3. Moisture ─────────────────────────────────────────────────────
+    # Mars has almost no moisture — thin atmosphere
+    moisture = np.zeros((GRID_H, GRID_W), dtype=np.float32)
+    # Slight moisture near polar caps (water ice)
+    moisture[polar_rows, :] = 0.1
+
+    # ── 4. Terraformed water bodies ─────────────────────────────────────
+    # Lore: 800 years of partial terraforming. Water pools in the deepest
+    # basins (Hellas, Utopia, Isidis). ~2% of surface is now liquid water.
+    from sol_data.shared import compute_sea_level as _compute_sl
+    from scipy.ndimage import binary_dilation
+
+    TERRAFORM_OCEAN_FRAC = 0.02  # 2% water coverage
+    sea_level = _compute_sl(elevation, TERRAFORM_OCEAN_FRAC)
+    surface_water = elevation < sea_level
+
+    # Don't flood polar regions — those stay as ice caps, not lakes
+    surface_water[polar_rows, :] = False
+
+    n_water = int(surface_water.sum())
+    print(f"    terraformed water: {n_water} cells "
+          f"(sea_level={sea_level:.4f})")
+
+    # ── 5. Biome classification ─────────────────────────────────────────
+    # Mars biome is built directly — compute_biome() would classify
+    # everything as ice at these temperatures.
+    biome = np.full((GRID_H, GRID_W), FERRIC_DUST, dtype=np.int8)
+
+    # Elevation-based ferric variation
+    biome[elevation > 0.55] = FERRIC_HIGHLAND  # volcanic highlands
+    biome[elevation < 0.25] = FERRIC_LOWLAND   # basin floors
+
+    # Polar ice caps
+    biome[polar_rows, :] = 17  # ice/snow
+
+    # Terraformed green fringe around water bodies — vegetation band
+    # where the thicker local atmosphere and water access allow plants.
+    # ~5 cell band around each water body.
+    veg_ring = binary_dilation(surface_water, iterations=5) & ~surface_water
+    # Don't put vegetation at poles
+    veg_ring[polar_rows, :] = False
+    biome[veg_ring] = 12  # shrubland (olive green — sparse terraformed vegetation)
+
+    # Inner vegetation ring (closer to water = lusher)
+    inner_ring = binary_dilation(surface_water, iterations=2) & ~surface_water
+    inner_ring[polar_rows, :] = False
+    biome[inner_ring] = 8  # temperate grassland (greener)
+
+    # Ocean depth bands for water bodies
+    if surface_water.any():
+        depth = np.clip((sea_level - elevation) / (sea_level + 1e-9), 0, 1)
+        biome[surface_water & (depth < 0.15)] = 2   # shallow
+        biome[surface_water & (depth >= 0.15) & (depth < 0.50)] = 1  # mid
+        biome[surface_water & (depth >= 0.50)] = 0   # deep
+
+    n_ice = int((biome == 17).sum())
+    n_ferric = int(((biome >= 34) & (biome <= 36)).sum())
+    n_veg = int(((biome == 8) | (biome == 12)).sum())
+    n_ocean = int(((biome >= 0) & (biome <= 2)).sum())
+    print(f"    biomes: {len(np.unique(biome))} classes "
+          f"(ferric={n_ferric}, ice={n_ice}, veg={n_veg}, water={n_ocean})")
+
+    # ── 5. Hillshade ────────────────────────────────────────────────────
+    hillshade = compute_hillshade(elevation)
+
+    # ── 6. Assemble ─────────────────────────────────────────────────────
+    return assemble_terrain(
+        elevation=elevation,
+        temperature_K=temperature_K,
+        moisture=moisture,
+        biome=biome,
+        surface_water=surface_water,
+        hillshade=hillshade,
+        rivers=[],  # no rivers on Mars
+        sea_level=sea_level,
+    )

tooling/planet-gen/sol_data/mercury.py

@@ -0,0 +1,120 @@
+"""
+Mercury (GJ0b) terrain builder.
+
+Data source:
+  - Elevation: MESSENGER DEM from USGS Astrogeology
+    665m/px global DEM, GeoTIFF.
+
+Mercury properties:
+  - Min elevation: ~-5380 m
+  - Max elevation: ~4480 m
+  - No atmosphere, no water
+  - Extreme temperature range: ~100K (night) to ~700K (day)
+  - body_type: "planet", planet_class: "barren", atmosphere: "none"
+"""
+
+import numpy as np
+from pathlib import Path
+
+from sol_data.download import ensure_cached
+from sol_data.shared import (
+    GRID_W, GRID_H,
+    load_tiff_as_array, load_raw_binary,
+    resample_to_grid, normalize_01,
+    compute_hillshade, assemble_terrain,
+    temperature_grid_analytical,
+)
+
+# MESSENGER DEM — try PDS binary first (compact), fall back to USGS GeoTIFF
+MESSENGER_PDS_URL = "https://pds-geosciences.wustl.edu/messenger/mess-h-mdis_mla-6-dem-elevation-v1/messdmdem_1001/data/global_dem_16ppd.img"
+MESSENGER_PDS_FILE = "messenger_dem_16ppd.img"
+MESSENGER_PDS_W = 5760
+MESSENGER_PDS_H = 2880
+
+# USGS GeoTIFF fallback (~506 MB, but PIL-loadable)
+MESSENGER_TIFF_URL = "https://planetarymaps.usgs.gov/mosaic/Mercury_Messenger_USGS_DEM_Global_665m_v2.tif"
+MESSENGER_TIFF_FILE = "Mercury_Messenger_USGS_DEM_Global_665m_v2.tif"
+
+MERCURY_MIN_ELEV_M = -5380.0
+MERCURY_MAX_ELEV_M = 4480.0
+MERCURY_EQUATORIAL_TEMP_K = 440.0  # mean dayside
+MERCURY_POLAR_TEMP_K = 200.0
+
+
+def _load_messenger() -> np.ndarray:
+    """Load MESSENGER DEM, return elevation in metres."""
+    # Try PDS binary first (compact ~33 MB)
+    try:
+        path = ensure_cached(MESSENGER_PDS_URL, MESSENGER_PDS_FILE)
+        print(f"    loading MESSENGER PDS: {path}")
+        arr = load_raw_binary(str(path), MESSENGER_PDS_W, MESSENGER_PDS_H,
+                              dtype=">i2", offset=0)
+        arr[arr > 20000] = 0.0
+        arr[arr < -20000] = 0.0
+        print(f"    MESSENGER range: [{arr.min():.0f}, {arr.max():.0f}] m")
+        return arr
+    except Exception as e:
+        print(f"    PDS load failed ({e}), trying USGS GeoTIFF...")
+
+    # Fallback: USGS GeoTIFF (~506 MB)
+    try:
+        path = ensure_cached(MESSENGER_TIFF_URL, MESSENGER_TIFF_FILE)
+        print(f"    loading MESSENGER GeoTIFF: {path}")
+        arr = load_tiff_as_array(str(path))
+        arr[arr < -20000] = 0.0
+        print(f"    MESSENGER shape: {arr.shape}, range: [{arr.min():.0f}, {arr.max():.0f}] m")
+        return arr
+    except Exception as e2:
+        print(f"    GeoTIFF also failed ({e2}), using procedural")
+        return None
+
+
+def build_terrain(body_def: dict) -> dict:
+    """Build Mercury terrain dict from MESSENGER data."""
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+    from planet_simulation import compute_biome
+
+    print("  Mercury: loading MESSENGER data...")
+
+    # ── 1. Elevation ────────────────────────────────────────────────────
+    raw = _load_messenger()
+    if raw is None:
+        import sys
+        sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+        from planet_simulation import simulate
+        return simulate(body_def)
+
+    from sol_data.shared import greenwich_to_dateline
+    shifted = greenwich_to_dateline(raw)
+    elevation_m = resample_to_grid(shifted, GRID_H, GRID_W, order=1)
+    elevation = normalize_01(elevation_m, MERCURY_MIN_ELEV_M, MERCURY_MAX_ELEV_M)
+
+    sea_level = 0.0
+    surface_water = np.zeros((GRID_H, GRID_W), dtype=bool)
+
+    # ── 2. Temperature ──────────────────────────────────────────────────
+    temperature_K = temperature_grid_analytical(
+        base_T_K=MERCURY_EQUATORIAL_TEMP_K,
+        elevation=elevation,
+        lapse_rate_K_per_unit=20.0,
+        lat_gradient_K=240.0,
+    )
+    temperature_K = np.maximum(temperature_K, 100.0)
+
+    # ── 3. Moisture ─────────────────────────────────────────────────────
+    moisture = np.zeros((GRID_H, GRID_W), dtype=np.float32)
+
+    # ── 4. Biome ────────────────────────────────────────────────────────
+    biome = compute_biome(body_def, elevation, sea_level, surface_water,
+                          temperature_K, moisture)
+
+    # ── 5. Hillshade ────────────────────────────────────────────────────
+    hillshade = compute_hillshade(elevation)
+
+    return assemble_terrain(
+        elevation=elevation, temperature_K=temperature_K,
+        moisture=moisture, biome=biome,
+        surface_water=surface_water, hillshade=hillshade,
+        rivers=[], sea_level=sea_level,
+    )

tooling/planet-gen/sol_data/shared.py

@@ -0,0 +1,249 @@
+"""
+Shared utilities for loading and processing real-world planetary data.
+
+All loaders produce arrays compatible with the planet_simulation terrain dict:
+  - Grid size: GRID_H x GRID_W (256 x 512)
+  - Elevation: float32 [0, 1] normalised
+  - Temperature: float32 in absolute Kelvin (normalised to [0,1] later)
+  - Moisture: float32 [0, 1]
+  - Sea level: float elevation threshold
+"""
+
+import math
+import struct
+import numpy as np
+from scipy.ndimage import zoom
+from PIL import Image
+
+# Planetary DEMs can exceed PIL's default decompression bomb limit
+Image.MAX_IMAGE_PIXELS = None
+
+# Match planet_simulation grid
+GRID_W = 512
+GRID_H = 256
+
+
+# ─── Loading ────────────────────────────────────────────────────────────────
+
+def load_tiff_as_array(path: str) -> np.ndarray:
+    """
+    Load a GeoTIFF/TIFF as a numpy array via PIL.
+
+    PIL handles uncompressed and LZW-compressed TIFFs with 8/16/32-bit
+    integer or float samples. For multi-band, returns (H, W, bands).
+    For single-band, returns (H, W).
+    """
+    img = Image.open(path)
+    arr = np.array(img, dtype=np.float32)
+    return arr
+
+
+def load_raw_binary(path: str, width: int, height: int,
+                    dtype: str = ">i2", offset: int = 0) -> np.ndarray:
+    """
+    Load a raw binary raster (PDS IMG, .bin, etc).
+
+    Parameters
+    ----------
+    path   : file path
+    width  : number of columns
+    height : number of rows
+    dtype  : numpy dtype string (e.g. ">i2" for big-endian int16)
+    offset : byte offset to skip (header size)
+    """
+    dt = np.dtype(dtype)
+    expected_bytes = width * height * dt.itemsize
+    with open(path, "rb") as f:
+        f.seek(offset)
+        raw = f.read(expected_bytes)
+    if len(raw) < expected_bytes:
+        raise ValueError(
+            f"Expected {expected_bytes} bytes, got {len(raw)}. "
+            f"Check width/height/dtype/offset."
+        )
+    arr = np.frombuffer(raw, dtype=dt).reshape(height, width).astype(np.float32)
+    return arr
+
+
+def load_image_as_elevation(path: str, invert: bool = False) -> np.ndarray:
+    """
+    Load a greyscale or RGB image and convert to float32 elevation.
+    For RGB, uses luminance. For greyscale, uses the single channel.
+    """
+    img = Image.open(path).convert("L")
+    arr = np.array(img, dtype=np.float32) / 255.0
+    if invert:
+        arr = 1.0 - arr
+    return arr
+
+
+# ─── Resampling ─────────────────────────────────────────────────────────────
+
+def resample_to_grid(arr: np.ndarray, target_h: int = GRID_H,
+                     target_w: int = GRID_W,
+                     order: int = 1) -> np.ndarray:
+    """
+    Resample a 2D array to target grid size.
+
+    order: 0=nearest, 1=bilinear, 3=cubic
+    """
+    if arr.shape == (target_h, target_w):
+        return arr.astype(np.float32)
+    zoom_y = target_h / arr.shape[0]
+    zoom_x = target_w / arr.shape[1]
+    return zoom(arr, (zoom_y, zoom_x), order=order).astype(np.float32)
+
+
+# ─── Normalisation ──────────────────────────────────────────────────────────
+
+def normalize_01(arr: np.ndarray, lo: float = None, hi: float = None) -> np.ndarray:
+    """Normalise array to [0, 1]."""
+    if lo is None:
+        lo = float(arr.min())
+    if hi is None:
+        hi = float(arr.max())
+    if hi - lo < 1e-9:
+        return np.zeros_like(arr, dtype=np.float32)
+    return np.clip((arr - lo) / (hi - lo), 0.0, 1.0).astype(np.float32)
+
+
+def compute_sea_level(elevation: np.ndarray, ocean_fraction: float) -> float:
+    """
+    Compute sea_level threshold such that ocean_fraction of cells are below it.
+    """
+    if ocean_fraction <= 0.0:
+        return 0.0
+    if ocean_fraction >= 1.0:
+        return 1.0
+    return float(np.percentile(elevation, ocean_fraction * 100.0))
+
+
+# ─── Longitude shift ────────────────────────────────────────────────────────
+
+def shift_longitude(arr: np.ndarray, shift_cols: int) -> np.ndarray:
+    """
+    Roll array along the longitude (column) axis.
+
+    The pipeline uses col 0 = 180°W. If source data uses col 0 = 0° (Greenwich),
+    shift by half the width to align.
+    """
+    return np.roll(arr, shift_cols, axis=1)
+
+
+def greenwich_to_dateline(arr: np.ndarray) -> np.ndarray:
+    """
+    Shift from col 0 = 0° (Greenwich) to col 0 = 180°W (dateline).
+    Standard for most NASA/NOAA global datasets → pipeline convention.
+    """
+    return shift_longitude(arr, arr.shape[1] // 2)
+
+
+# ─── Hillshade ──────────────────────────────────────────────────────────────
+
+def compute_hillshade(elevation: np.ndarray,
+                      sun_azimuth_deg: float = 315.0,
+                      sun_altitude_deg: float = 45.0) -> np.ndarray:
+    """
+    Compute hillshade from elevation grid. Matches planet_simulation.compute_hillshade().
+    """
+    scale = elevation.shape[1] / 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.sin(az)
+    ly = -math.cos(alt) * math.cos(az)
+    lz = math.sin(alt)
+
+    shade = np.clip(nx * lx + ny * ly + nz * lz, 0.0, 1.0)
+    return shade.astype(np.float32)
+
+
+# ─── Analytical temperature models ─────────────────────────────────────────
+
+def temperature_equilibrium_K(luminosity_solar: float, distance_au: float,
+                              albedo: float = 0.3) -> float:
+    """
+    Stefan-Boltzmann equilibrium temperature in Kelvin.
+    """
+    L_sun = 3.828e26  # watts
+    sigma = 5.670e-8
+    d_m = distance_au * 1.496e11
+    T_eq = ((luminosity_solar * L_sun * (1 - albedo)) /
+            (16 * math.pi * sigma * d_m**2)) ** 0.25
+    return T_eq
+
+
+def temperature_grid_analytical(
+    base_T_K: float,
+    grid_h: int = GRID_H,
+    grid_w: int = GRID_W,
+    elevation: np.ndarray = None,
+    lapse_rate_K_per_unit: float = 40.0,
+    lat_gradient_K: float = 60.0,
+) -> np.ndarray:
+    """
+    Analytical temperature grid: equator-to-pole gradient + elevation lapse.
+
+    Parameters
+    ----------
+    base_T_K             : equatorial temperature in Kelvin
+    elevation            : normalised [0,1] elevation grid (optional)
+    lapse_rate_K_per_unit: temperature drop per unit elevation
+    lat_gradient_K       : total temperature drop from equator to pole
+    """
+    v = np.linspace(0, 1, grid_h, dtype=np.float32)
+    lat_frac = np.abs(v - 0.5) * 2.0  # 0 at equator, 1 at poles
+    lat_temp = lat_frac[:, np.newaxis] * lat_gradient_K  # broadcast to (H, W)
+    temp = np.full((grid_h, grid_w), base_T_K, dtype=np.float32)
+    temp -= lat_temp
+    if elevation is not None:
+        temp -= elevation * lapse_rate_K_per_unit
+    return temp
+
+
+# ─── Terrain dict assembly ──────────────────────────────────────────────────
+
+def assemble_terrain(
+    elevation: np.ndarray,
+    temperature_K: np.ndarray,
+    moisture: np.ndarray,
+    biome: np.ndarray,
+    surface_water: np.ndarray,
+    hillshade: np.ndarray,
+    rivers: list,
+    sea_level: float,
+) -> dict:
+    """
+    Assemble the terrain dict in the format expected by render_heightmap
+    and render_globe. Temperature is normalised to [0,1] for the output
+    (matching planet_simulation.simulate() lines 891-893).
+    """
+    H, W = elevation.shape
+    river_grid = np.zeros((H, W), dtype=bool)
+    for path in rivers:
+        for r, c in path:
+            if 0 <= r < H and 0 <= c < W:
+                river_grid[r, c] = True
+
+    # Normalise temperature to [0,1] for renderer display
+    t_min, t_max = temperature_K.min(), temperature_K.max()
+    temp_norm = ((temperature_K - t_min) / (t_max - t_min + 1e-9)).astype(np.float32)
+
+    return {
+        "elevation":     elevation.astype(np.float32),
+        "temperature":   temp_norm,
+        "moisture":      moisture.astype(np.float32),
+        "biome":         biome.astype(np.int8),
+        "surface_water": surface_water.astype(bool),
+        "hillshade":     hillshade.astype(np.float32),
+        "river_grid":    river_grid,
+        "rivers":        rivers,
+        "sea_level":     float(sea_level),
+        "_grid_w":       W,
+        "_grid_h":       H,
+    }

tooling/planet-gen/sol_data/titan.py

@@ -0,0 +1,175 @@
+"""
+Titan (GJ0g-1) terrain builder.
+
+Titan is unique: dense nitrogen atmosphere, methane rain cycle,
+methane/ethane lakes and rivers. Surface temperature ~94K uniform.
+
+Data source:
+  - Surface: Cassini ISS global mosaic (4km resolution)
+  - Topography: very sparse Cassini radar altimetry (gap-filled)
+
+Since Cassini topographic data is extremely sparse, we use the ISS
+mosaic albedo to derive synthetic elevation (similar to ice moons)
+with special handling for known methane lake regions.
+
+Properties:
+  - planet_class: "frozen", atmosphere: "dense", hydrosphere: "rivers"
+  - Methane lakes primarily near the north pole (Kraken Mare, Ligeia Mare)
+"""
+
+import numpy as np
+from pathlib import Path
+from scipy.ndimage import gaussian_filter
+
+from sol_data.download import ensure_cached
+from sol_data.shared import (
+    GRID_W, GRID_H,
+    load_image_as_elevation, resample_to_grid, normalize_01,
+    compute_hillshade, assemble_terrain,
+    temperature_grid_analytical,
+)
+
+# Cassini ISS global mosaic
+TITAN_MOSAIC_URL = "https://astrogeology.usgs.gov/cache/images/5e5ba96a58d3b38ee6e7b1e94b8c44e6_titan_iss_p19658_mosaic_global_4km.jpg"
+TITAN_MOSAIC_FILE = "titan_cassini_iss_mosaic.jpg"
+
+TITAN_SURFACE_TEMP_K = 94.0  # nearly uniform
+TITAN_METHANE_LAKE_FRACTION = 0.02  # ~2% of surface is liquid methane
+
+
+def _load_titan_mosaic() -> np.ndarray:
+    """Load Titan mosaic and convert to synthetic elevation."""
+    try:
+        path = ensure_cached(TITAN_MOSAIC_URL, TITAN_MOSAIC_FILE)
+        print(f"    loading Titan mosaic: {path}")
+        albedo = load_image_as_elevation(str(path), invert=False)
+        albedo = resample_to_grid(albedo, GRID_H, GRID_W, order=1)
+    except Exception as e:
+        print(f"    WARNING: Titan mosaic unavailable ({e}), synthetic")
+        albedo = _synthetic_titan_terrain()
+
+    # Dark regions = low (lakes/flat), bright = dunes/highlands
+    elevation = gaussian_filter(albedo, sigma=3.0)
+    return normalize_01(elevation)
+
+
+def _synthetic_titan_terrain() -> np.ndarray:
+    """Generate synthetic Titan terrain."""
+    rng = np.random.default_rng(94)
+    base = rng.random((GRID_H, GRID_W)).astype(np.float32)
+    base = gaussian_filter(base, sigma=6.0)
+    # Titan has equatorial dune fields (higher terrain)
+    v = np.linspace(0, 1, GRID_H, dtype=np.float32)
+    equatorial = np.exp(-((v - 0.5) ** 2) / 0.02)
+    base += equatorial[:, np.newaxis] * 0.3
+    return normalize_01(base)
+
+
+def build_terrain(body_def: dict) -> dict:
+    """Build Titan terrain dict."""
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+    from planet_simulation import compute_biome
+
+    print("  Titan: loading data...")
+
+    # ── 1. Elevation ────────────────────────────────────────────────────
+    elevation = _load_titan_mosaic()
+
+    # Titan has methane lakes — set sea level to create them
+    # Lakes are concentrated at north polar regions
+    # Use a low sea level so that only the darkest (lowest) areas become liquid
+    from sol_data.shared import compute_sea_level
+    sea_level = compute_sea_level(elevation, TITAN_METHANE_LAKE_FRACTION)
+    surface_water = elevation < sea_level
+
+    # Concentrate lakes near north pole (real Titan has lakes mostly 60-90°N)
+    v = np.linspace(0, 1, GRID_H, dtype=np.float32)
+    lat_abs = np.abs(v - 0.5) * 2.0  # 0=equator, 1=poles
+    north_mask = v < 0.2  # north polar region (top 20% of grid = 72-90°N)
+    # Allow lakes only in polar regions — mask out equatorial/southern "seas"
+    equatorial_mask = (lat_abs < 0.6)[:, np.newaxis] * np.ones(GRID_W, dtype=bool)
+    surface_water = surface_water & ~equatorial_mask
+
+    print(f"    methane lakes: {surface_water.sum()} cells")
+
+    # ── 2. Temperature ──────────────────────────────────────────────────
+    # Titan has nearly uniform surface temp due to dense atmosphere + distance
+    temperature_K = np.full((GRID_H, GRID_W), TITAN_SURFACE_TEMP_K, dtype=np.float32)
+    # Very slight pole-equator gradient (~2K)
+    lat_temp = lat_abs[:, np.newaxis] * 2.0
+    temperature_K -= lat_temp
+
+    # ── 3. Moisture ─────────────────────────────────────────────────────
+    # Titan has a methane humidity cycle — higher moisture near poles
+    moisture = np.zeros((GRID_H, GRID_W), dtype=np.float32)
+    # Polar moisture (methane humidity)
+    polar_humid = np.clip(lat_abs[:, np.newaxis] - 0.5, 0, 1) * 0.6
+    moisture += polar_humid
+    # Some equatorial humidity (methane drizzle)
+    equatorial_humid = np.exp(-((v[:, np.newaxis] - 0.5) ** 2) / 0.05) * 0.2
+    moisture += equatorial_humid
+
+    # ── 4. Biome ────────────────────────────────────────────────────────
+    # Titan at 94K with dense atmosphere goes through Whittaker table
+    # Everything will classify as ice/snow (class 17) — which is correct
+    biome = compute_biome(body_def, elevation, sea_level, surface_water,
+                          temperature_K, moisture)
+
+    # Override: methane lakes should be ocean classes, not ice
+    # (The biome function sets ocean depth bands for surface_water, which is
+    #  what we want — methane lakes rendered like ocean)
+    print(f"    biomes: {len(np.unique(biome))} classes")
+
+    # ── 5. Hillshade ────────────────────────────────────────────────────
+    hillshade = compute_hillshade(elevation)
+
+    # ── 6. Rivers ───────────────────────────────────────────────────────
+    # Titan has methane drainage channels — add synthetic ones near poles
+    rivers = _titan_rivers(elevation, surface_water)
+
+    return assemble_terrain(
+        elevation=elevation, temperature_K=temperature_K,
+        moisture=moisture, biome=biome,
+        surface_water=surface_water, hillshade=hillshade,
+        rivers=rivers, sea_level=sea_level,
+    )
+
+
+def _titan_rivers(elevation: np.ndarray, surface_water: np.ndarray) -> list:
+    """
+    Generate synthetic methane drainage channels for Titan.
+    Simple downhill tracing from high-latitude sources to lakes.
+    """
+    rivers = []
+    rng = np.random.default_rng(94)
+
+    # Start from a few points in the north polar region
+    for _ in range(5):
+        r = int(rng.integers(10, 50))  # north polar zone
+        c = int(rng.integers(0, GRID_W))
+        path = [(r, c)]
+        visited = {(r, c)}
+
+        for _ in range(200):
+            if surface_water[r, c]:
+                break
+            best_r, best_c = r, c
+            best_elev = elevation[r, c]
+            for dr, dc in [(-1, 0), (1, 0), (0, -1), (0, 1),
+                           (-1, -1), (-1, 1), (1, -1), (1, 1)]:
+                nr, nc = r + dr, (c + dc) % GRID_W
+                if 0 <= nr < GRID_H and (nr, nc) not in visited:
+                    if elevation[nr, nc] < best_elev:
+                        best_elev = elevation[nr, nc]
+                        best_r, best_c = nr, nc
+            if (best_r, best_c) == (r, c):
+                break
+            r, c = int(best_r), int(best_c)
+            path.append((r, c))
+            visited.add((r, c))
+
+        if len(path) >= 5:
+            rivers.append(path)
+
+    return rivers

tooling/planet-gen/sol_data/venus.py

@@ -0,0 +1,126 @@
+"""
+Venus (GJ0c) terrain builder.
+
+Data source:
+  - Elevation: Magellan radar altimetry from USGS Astrogeology
+    Global topography at ~4.6 km/px, PDS format.
+
+Venus properties:
+  - Surface: volcanic, extremely hot (~735K), dense CO2 atmosphere
+  - No liquid water, thick clouds
+  - Min elevation: ~-2000 m (lowlands)
+  - Max elevation: ~11000 m (Maxwell Montes on Ishtar Terra)
+  - planet_class: "volcanic", atmosphere: "toxic"
+"""
+
+import numpy as np
+from pathlib import Path
+
+from sol_data.download import ensure_cached
+from sol_data.shared import (
+    GRID_W, GRID_H,
+    load_tiff_as_array, load_raw_binary, resample_to_grid, normalize_01,
+    compute_hillshade, assemble_terrain,
+    temperature_grid_analytical,
+)
+
+# Magellan topography — USGS GeoTIFF (reliable, PIL-loadable)
+MAGELLAN_TIFF_URL = "https://planetarymaps.usgs.gov/mosaic/Venus_Magellan_Topography_Global_4641m_v02.tif"
+MAGELLAN_TIFF_FILE = "Venus_Magellan_Topography_Global_4641m_v02.tif"
+
+# PDS fallback (raw binary, dimensions may vary)
+MAGELLAN_PDS_URL = "https://pds-geosciences.wustl.edu/mgn/mgn-v-rdrs-5-dim-v1/mg_3002/gedr/gtdr/gtdr_shtplt.img"
+MAGELLAN_PDS_FILE = "venus_magellan_gtdr.img"
+
+VENUS_MIN_ELEV_M = -2000.0
+VENUS_MAX_ELEV_M = 11000.0
+VENUS_SURFACE_TEMP_K = 735.0  # nearly uniform due to dense atmosphere
+
+
+def _load_magellan() -> np.ndarray:
+    """Load Magellan topography data."""
+    # Try USGS GeoTIFF first (reliable, well-defined format)
+    try:
+        path = ensure_cached(MAGELLAN_TIFF_URL, MAGELLAN_TIFF_FILE)
+        print(f"    loading Magellan GeoTIFF: {path}")
+        arr = load_tiff_as_array(str(path))
+        # Handle nodata
+        arr[arr < -20000] = 0.0
+        arr[arr > 20000] = 0.0
+        print(f"    Magellan shape: {arr.shape}, "
+              f"range: [{arr.min():.0f}, {arr.max():.0f}] m")
+        return arr
+    except Exception as e:
+        print(f"    GeoTIFF failed ({e}), trying PDS binary...")
+
+    # PDS fallback — try common dimension/format combinations
+    try:
+        path = ensure_cached(MAGELLAN_PDS_URL, MAGELLAN_PDS_FILE)
+        print(f"    loading Magellan PDS: {path}")
+        for w, h in [(4096, 2048), (2048, 1024), (8192, 4096)]:
+            try:
+                arr = load_raw_binary(str(path), w, h, dtype=">i2", offset=0)
+                arr[arr > 20000] = 0.0
+                arr[arr < -20000] = 0.0
+                print(f"    Magellan PDS: {w}x{h}, range: [{arr.min():.0f}, {arr.max():.0f}]")
+                return arr
+            except ValueError:
+                continue
+    except Exception as e3:
+        print(f"    PDS also failed ({e3})")
+
+    # All sources failed — fall through to procedural generation
+    print(f"    WARNING: all Magellan sources failed, using procedural")
+    return None
+
+
+def build_terrain(body_def: dict) -> dict:
+    """Build Venus terrain dict from Magellan data."""
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+    from planet_simulation import compute_biome
+
+    print("  Venus: loading Magellan data...")
+
+    # ── 1. Elevation ────────────────────────────────────────────────────
+    raw = _load_magellan()
+    if raw is None:
+        # Fall back to procedural simulation
+        import sys
+        sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+        from planet_simulation import simulate
+        return simulate(body_def)
+
+    from sol_data.shared import greenwich_to_dateline
+    shifted = greenwich_to_dateline(raw)
+    elevation_m = resample_to_grid(shifted, GRID_H, GRID_W, order=1)
+    elevation = normalize_01(elevation_m, VENUS_MIN_ELEV_M, VENUS_MAX_ELEV_M)
+
+    sea_level = 0.0
+    surface_water = np.zeros((GRID_H, GRID_W), dtype=bool)
+
+    # ── 2. Temperature ──────────────────────────────────────────────────
+    # Venus has nearly uniform surface temperature due to dense atmosphere
+    temperature_K = temperature_grid_analytical(
+        base_T_K=VENUS_SURFACE_TEMP_K,
+        elevation=elevation,
+        lapse_rate_K_per_unit=50.0,  # slight cooling at altitude
+        lat_gradient_K=5.0,          # almost no lat variation (thick atmo)
+    )
+
+    # ── 3. Moisture ─────────────────────────────────────────────────────
+    moisture = np.zeros((GRID_H, GRID_W), dtype=np.float32)
+
+    # ── 4. Biome ────────────────────────────────────────────────────────
+    biome = compute_biome(body_def, elevation, sea_level, surface_water,
+                          temperature_K, moisture)
+
+    # ── 5. Hillshade ────────────────────────────────────────────────────
+    hillshade = compute_hillshade(elevation)
+
+    return assemble_terrain(
+        elevation=elevation, temperature_K=temperature_K,
+        moisture=moisture, biome=biome,
+        surface_water=surface_water, hillshade=hillshade,
+        rivers=[], sea_level=sea_level,
+    )

tooling/planet-gen/sol_import.py

@@ -0,0 +1,370 @@
+#!/usr/bin/env python3
+"""
+sol_import.py — Import real-world data for the Sol system (GJ-0).
+
+Produces the same output format as generate.py (heightmap.png, globe.png,
+markers.json, terrain.npz) by constructing terrain dicts from real
+planetary science data instead of procedural simulation.
+
+Usage:
+    python3 sol_import.py                        # All Sol bodies
+    python3 sol_import.py --body GJ0d            # Earth only
+    python3 sol_import.py --body GJ0d --body GJ0e  # Earth + Mars
+    python3 sol_import.py --download-only        # Fetch data, skip rendering
+    python3 sol_import.py --heightmap-size 2048x1024 --globe-size 1024
+
+Data is cached in tooling/planet-gen/sol_data/.cache/ after first download.
+"""
+
+import argparse
+import json
+import os
+import sys
+import time
+
+# Venv bootstrap — re-exec into .venv/bin/python if not already there.
+from pathlib import Path
+TOOLING_DIR = Path(__file__).resolve().parent
+WORKTREE_ROOT = (TOOLING_DIR / ".." / "..").resolve()
+_venv_python = WORKTREE_ROOT / ".venv" / "bin" / "python"
+if _venv_python.exists() and Path(sys.executable).resolve() != _venv_python.resolve():
+    os.execv(str(_venv_python), [str(_venv_python)] + sys.argv)
+
+import numpy as np
+
+from planet_simulation import simulate, compute_biome, compute_hillshade
+from render_heightmap import render_heightmap
+from generate import _build_markers
+
+# Per-body importers (lazy-loaded)
+SOL_INDEX = WORKTREE_ROOT / "wiki" / "star-systems" / "GJ-0" / "index.md"
+SOL_OVERRIDES = TOOLING_DIR / "sol_overrides.json"
+SOL_BODIES_DIR = WORKTREE_ROOT / "wiki" / "star-systems" / "GJ-0" / "bodies"
+SOL_MARKERS_DIR = TOOLING_DIR / "sol_markers"
+
+# Bodies that use real-world data (keyed by body_id → importer module)
+REAL_DATA_BODIES = {
+    "GJ0b":   "mercury",
+    "GJ0c":   "venus",
+    "GJ0d":   "earth",
+    "GJ0d-1": "luna",
+    "GJ0e":   "mars",
+    "GJ0f-1": "io_moon",
+    "GJ0f-2": "ice_moons",
+    "GJ0f-3": "ice_moons",
+    "GJ0f-4": "ice_moons",
+    "GJ0g-1": "titan",
+    "GJ0g-2": "ice_moons",
+}
+
+# Bodies that fall through to procedural simulation
+PROCEDURAL_BODIES = {"GJ0e-1", "GJ0e-2"}
+
+# Non-renderable body types
+SKIP_TYPES = {"asteroid_belt", "oort_cloud"}
+
+
+def _load_importer(module_name: str):
+    """Lazy-import a sol_data.* module."""
+    import importlib
+    return importlib.import_module(f"sol_data.{module_name}")
+
+
+def _apply_named_features(markers: dict, body_id: str) -> dict:
+    """Overlay named features from sol_markers/ onto auto-detected markers."""
+    features_map = {
+        "GJ0d":   "earth_features.json",
+        "GJ0e":   "mars_features.json",
+        "GJ0d-1": "luna_features.json",
+    }
+    outer_bodies = {"GJ0f-1", "GJ0f-2", "GJ0f-3", "GJ0f-4",
+                    "GJ0g-1", "GJ0g-2"}
+
+    filename = features_map.get(body_id)
+    if not filename and body_id in outer_bodies:
+        filename = "outer_features.json"
+
+    if not filename:
+        return markers
+
+    features_path = SOL_MARKERS_DIR / filename
+    if not features_path.exists():
+        return markers
+
+    with open(features_path) as f:
+        features = json.load(f)
+
+    body_features = features.get(body_id, features)
+
+    # Name auto-detected oceans by matching center coordinates
+    if "oceans" in body_features:
+        for named_ocean in body_features["oceans"]:
+            best_match = None
+            best_dist = float("inf")
+            nc = named_ocean["center"]
+            for detected in markers["oceans"]:
+                dc = detected["center"]
+                dist = (dc[0] - nc[0])**2 + (dc[1] - nc[1])**2
+                if dist < best_dist:
+                    best_dist = dist
+                    best_match = detected
+            if best_match and best_dist < 2500:  # within ~50 cells
+                best_match["name"] = named_ocean["name"]
+
+    # Name auto-detected mountain ranges by matching peak coordinates
+    if "mountain_ranges" in body_features:
+        for named_range in body_features["mountain_ranges"]:
+            best_match = None
+            best_dist = float("inf")
+            nc = named_range.get("peak", named_range.get("center", [0, 0]))
+            for detected in markers["mountain_ranges"]:
+                dp = detected.get("peak", detected.get("center", [0, 0]))
+                dist = (dp[0] - nc[0])**2 + (dp[1] - nc[1])**2
+                if dist < best_dist:
+                    best_dist = dist
+                    best_match = detected
+            if best_match and best_dist < 1600:  # within ~40 cells
+                best_match["name"] = named_range["name"]
+
+    # Name rivers by matching start/end coordinates
+    if "rivers" in body_features:
+        for named_river in body_features["rivers"]:
+            best_match = None
+            best_dist = float("inf")
+            nc = named_river.get("mouth", named_river.get("center", [0, 0]))
+            for detected in markers["rivers"]:
+                if not detected["path"]:
+                    continue
+                # Check last point (mouth) of river path
+                dp = detected["path"][-1]
+                dist = (dp[0] - nc[0])**2 + (dp[1] - nc[1])**2
+                if dist < best_dist:
+                    best_dist = dist
+                    best_match = detected
+            if best_match and best_dist < 900:
+                best_match["name"] = named_river["name"]
+
+    # Add cities as POIs
+    if "cities" in body_features:
+        for city in body_features["cities"]:
+            markers["cities"].append({
+                "id": f"city_{city['name'].lower().replace(' ', '_')}",
+                "name": city["name"],
+                "center": city["center"],
+                "population": city.get("population"),
+            })
+
+    # Add POIs
+    if "pois" in body_features:
+        for poi in body_features["pois"]:
+            markers["pois"].append(poi)
+
+    return markers
+
+
+def _generate_body(body_def: dict, hmap_w: int, hmap_h: int,
+                   globe_size: int, render_mode: str, output_dir: Path,
+                   download_only: bool = False):
+    """Generate all outputs for a single Sol body."""
+    body_id = body_def["id"]
+    body_type = body_def.get("body_type", "planet")
+    planet_class = body_def.get("planet_class", "unknown")
+    name = body_def.get("name") or body_id
+
+    # Skip non-renderable types
+    if body_type in SKIP_TYPES:
+        print(f"\n  {body_id} ({name}) — skipped ({body_type})")
+        return
+
+    body_dir = output_dir / body_id
+    body_dir.mkdir(parents=True, exist_ok=True)
+
+    print(f"\n  {body_id} ({name}) — {planet_class}")
+
+    t0 = time.time()
+
+    # ── 1. Build terrain ────────────────────────────────────────────────
+    terrain = {}
+    is_gas = planet_class in ("gas_giant",) or body_type == "gas_giant"
+
+    if is_gas:
+        # Gas giants: no terrain, renderer handles bands procedurally
+        terrain = {}
+        print(f"    terrain: gas giant (procedural bands)")
+    elif body_id in REAL_DATA_BODIES:
+        # Real-world data import
+        module_name = REAL_DATA_BODIES[body_id]
+        print(f"    importing real data via sol_data.{module_name}...")
+        importer = _load_importer(module_name)
+        terrain = importer.build_terrain(body_def)
+        if download_only:
+            print(f"    download complete, skipping render")
+            return
+    elif body_id in PROCEDURAL_BODIES:
+        # Fall through to standard procedural simulation
+        print(f"    procedural simulation (irregular body)...")
+        terrain = simulate(body_def)
+    else:
+        print(f"    WARNING: no importer for {body_id}, using procedural")
+        terrain = simulate(body_def)
+
+    t_terrain = time.time()
+
+    if terrain:
+        print(f"    terrain: {t_terrain - t0:.1f}s  "
+              f"sea={terrain['sea_level']:.3f}  "
+              f"land={int((~terrain['surface_water']).sum())}  "
+              f"rivers={len(terrain['rivers'])}")
+    else:
+        print(f"    terrain: gas giant ({t_terrain - t0:.1f}s)")
+
+    # ── 2. Render heightmap ─────────────────────────────────────────────
+    t_hmap = t_terrain
+    if terrain:
+        hmap_img = render_heightmap(body_def, terrain,
+                                    out_w=hmap_w, out_h=hmap_h,
+                                    render_mode=render_mode, chrome=False)
+        hmap_img.save(str(body_dir / "heightmap.png"))
+        t_hmap = time.time()
+        print(f"    heightmap: {t_hmap - t_terrain:.1f}s  {hmap_w}x{hmap_h}")
+
+    # ── 3. Render globe ─────────────────────────────────────────────────
+    try:
+        from planet_renderer import render_globe
+        globe_img = render_globe(body_def, terrain, size=globe_size)
+        globe_img.save(str(body_dir / "globe.png"))
+        t_globe = time.time()
+        print(f"    globe: {t_globe - t_hmap:.1f}s  {globe_size}x{globe_size}")
+    except Exception as e:
+        print(f"    globe: FAILED — {e}")
+        t_globe = time.time()
+
+    # ── 4. Write data files ─────────────────────────────────────────────
+    if terrain:
+        # terrain.npz
+        save_dict = {}
+        for key in ("elevation", "temperature", "moisture", "hillshade",
+                    "biome", "surface_water", "river_grid"):
+            if key in terrain:
+                save_dict[key] = terrain[key]
+        save_dict["sea_level"] = np.array([terrain["sea_level"]])
+        np.savez_compressed(str(body_dir / "terrain.npz"), **save_dict)
+
+        # markers.json — auto-detected + named features overlay
+        markers = _build_markers(body_def, terrain)
+        markers = _apply_named_features(markers, body_id)
+        with open(body_dir / "markers.json", "w") as f:
+            json.dump(markers, f, indent=2)
+
+    # ── 5. Write index.md frontmatter ───────────────────────────────────
+    _write_index_md(body_def, body_dir)
+
+    elapsed = time.time() - t0
+    print(f"    total: {elapsed:.1f}s -> {body_dir}/")
+
+
+def _write_index_md(body_def: dict, body_dir: Path):
+    """Write body index.md with YAML frontmatter."""
+    import yaml
+
+    # Strip internal fields
+    bd = {k: v for k, v in body_def.items()
+          if not k.startswith("_") and k != "wiki"}
+
+    fm = yaml.dump(bd, default_flow_style=False, sort_keys=False,
+                   allow_unicode=True)
+
+    name = body_def.get("name") or body_def["id"]
+    planet_class = body_def.get("planet_class", "unknown")
+    system_link = "[GJ-0](../../index.md)"
+
+    md = f"""---
+{fm.rstrip()}
+---
+
+# {name}
+
+{planet_class.replace('_', ' ').title()} {'planet' if body_def.get('body_type') == 'planet' else body_def.get('body_type', 'body')}.
+
+**System:** {system_link}
+
+## Visual
+
+![Heightmap](heightmap.png)
+
+![Globe](globe.png)
+"""
+    with open(body_dir / "index.md", "w") as f:
+        f.write(md)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Sol system (GJ-0) real-world terrain importer")
+
+    parser.add_argument("--body", action="append", default=None,
+                        help="Specific body ID(s) to generate (repeatable)")
+    parser.add_argument("--download-only", action="store_true",
+                        help="Download source data without rendering")
+    parser.add_argument("--output-dir", default=None,
+                        help="Override output directory")
+    parser.add_argument("--heightmap-size", default="1024x512",
+                        help="Heightmap resolution (WxH)")
+    parser.add_argument("--globe-size", type=int, default=512,
+                        help="Globe resolution (square)")
+    parser.add_argument("--render-mode", choices=["cartographic", "photographic"],
+                        default="cartographic")
+
+    args = parser.parse_args()
+
+    # Parse heightmap size
+    try:
+        hw, hh = args.heightmap_size.lower().split("x")
+        hmap_w, hmap_h = int(hw), int(hh)
+    except ValueError:
+        print(f"error: invalid heightmap size '{args.heightmap_size}'",
+              file=sys.stderr)
+        sys.exit(1)
+
+    output_dir = Path(args.output_dir) if args.output_dir else SOL_BODIES_DIR
+
+    # Parse body definitions from GJ-0 index.md
+    from body_definition_parser import parse_system
+
+    overrides = {}
+    if SOL_OVERRIDES.exists():
+        with open(SOL_OVERRIDES) as f:
+            overrides = json.load(f)
+
+    body_defs = parse_system(str(SOL_INDEX), overrides=overrides)
+    print(f"Sol system: {len(body_defs)} bodies parsed")
+
+    # Filter to requested bodies
+    if args.body:
+        requested = set(args.body)
+        body_defs = [bd for bd in body_defs if bd["id"] in requested]
+        if not body_defs:
+            print(f"error: no matching bodies for {args.body}", file=sys.stderr)
+            sys.exit(1)
+
+    # Generate
+    t_total = time.time()
+    failed = []
+    for bd in body_defs:
+        try:
+            _generate_body(bd, hmap_w, hmap_h, args.globe_size,
+                           args.render_mode, output_dir,
+                           download_only=args.download_only)
+        except Exception as e:
+            print(f"\n  FAILED: {bd['id']} — {e}")
+            failed.append(bd["id"])
+
+    elapsed = time.time() - t_total
+    n_ok = len(body_defs) - len(failed)
+    print(f"\n  Done: {n_ok}/{len(body_defs)} bodies in {elapsed:.1f}s")
+    if failed:
+        print(f"  Failed: {', '.join(failed)}")
+
+
+if __name__ == "__main__":
+    main()

tooling/planet-gen/sol_markers/earth_features.json

@@ -0,0 +1,94 @@
+{
+    "GJ0d": {
+        "oceans": [
+            {"name": "Pacific Ocean",  "center": [128, 440]},
+            {"name": "Atlantic Ocean", "center": [128, 170]},
+            {"name": "Indian Ocean",   "center": [160, 330]},
+            {"name": "Arctic Ocean",   "center": [15, 256]},
+            {"name": "Southern Ocean", "center": [230, 256]}
+        ],
+        "mountain_ranges": [
+            {"name": "Himalayas",             "peak": [93, 350],  "center": [93, 348]},
+            {"name": "Andes",                 "peak": [118, 140], "center": [140, 140]},
+            {"name": "Rocky Mountains",       "peak": [88, 105],  "center": [85, 105]},
+            {"name": "Alps",                  "peak": [82, 264],  "center": [82, 264]},
+            {"name": "Urals",                 "peak": [68, 300],  "center": [72, 300]},
+            {"name": "Atlas Mountains",       "peak": [93, 254],  "center": [94, 254]},
+            {"name": "Great Dividing Range",  "peak": [160, 420], "center": [162, 420]}
+        ],
+        "rivers": [
+            {"name": "Danube",        "mouth": [82, 278]},
+            {"name": "Volga",         "mouth": [75, 298]},
+            {"name": "Rhine",         "mouth": [79, 264]},
+            {"name": "Mississippi",   "mouth": [96, 107]},
+            {"name": "St. Lawrence",  "mouth": [80, 132]},
+            {"name": "Amazon",        "mouth": [126, 165]},
+            {"name": "Paraná",        "mouth": [148, 155]},
+            {"name": "Nile",          "mouth": [98, 286]},
+            {"name": "Congo",         "mouth": [125, 268]},
+            {"name": "Tigris",        "mouth": [96, 303]},
+            {"name": "Yangtze",       "mouth": [97, 387]},
+            {"name": "Ganges",        "mouth": [102, 351]},
+            {"name": "Mekong",        "mouth": [114, 374]},
+            {"name": "Murray",        "mouth": [170, 417]}
+        ],
+        "cities": [
+            {"name": "London",           "center": [79, 260], "population": 9000000,   "region": "europe"},
+            {"name": "Istanbul",         "center": [83, 279], "population": 15000000,  "region": "europe"},
+            {"name": "Moscow",           "center": [72, 294], "population": 12700000,  "region": "europe"},
+            {"name": "Paris",            "center": [80, 261], "population": 11000000,  "region": "europe"},
+            {"name": "Berlin",           "center": [77, 269], "population": 3700000,   "region": "europe"},
+
+            {"name": "Mexico City",      "center": [107, 101], "population": 21800000, "region": "north_america"},
+            {"name": "New York",         "center": [87, 130],  "population": 20100000, "region": "north_america"},
+            {"name": "Los Angeles",      "center": [93, 95],   "population": 13200000, "region": "north_america"},
+            {"name": "Toronto",          "center": [84, 123],  "population": 6200000,  "region": "north_america"},
+            {"name": "Chicago",          "center": [85, 115],  "population": 9500000,  "region": "north_america"},
+
+            {"name": "São Paulo",        "center": [143, 164], "population": 22400000, "region": "south_america"},
+            {"name": "Lima",             "center": [133, 131], "population": 10700000, "region": "south_america"},
+            {"name": "Bogotá",           "center": [121, 135], "population": 11300000, "region": "south_america"},
+            {"name": "Rio de Janeiro",   "center": [142, 168], "population": 13500000, "region": "south_america"},
+            {"name": "Buenos Aires",     "center": [151, 153], "population": 15200000, "region": "south_america"},
+
+            {"name": "Lagos",            "center": [120, 262], "population": 15400000, "region": "africa"},
+            {"name": "Kinshasa",         "center": [124, 270], "population": 15600000, "region": "africa"},
+            {"name": "Cairo",            "center": [97, 286],  "population": 21300000, "region": "africa"},
+            {"name": "Johannesburg",     "center": [156, 279], "population": 6000000,  "region": "africa"},
+            {"name": "Nairobi",          "center": [128, 293], "population": 5100000,  "region": "africa"},
+
+            {"name": "Tehran",           "center": [92, 308],  "population": 9000000,  "region": "west_asia"},
+            {"name": "Baghdad",          "center": [94, 303],  "population": 8100000,  "region": "west_asia"},
+            {"name": "Riyadh",           "center": [103, 304], "population": 7700000,  "region": "west_asia"},
+            {"name": "Ankara",           "center": [87, 284],  "population": 5700000,  "region": "west_asia"},
+            {"name": "Karachi",          "center": [103, 327], "population": 16500000, "region": "west_asia"},
+
+            {"name": "Tokyo",            "center": [92, 400],  "population": 37400000, "region": "east_asia"},
+            {"name": "Delhi",            "center": [99, 339],  "population": 32900000, "region": "east_asia"},
+            {"name": "Shanghai",         "center": [97, 387],  "population": 28500000, "region": "east_asia"},
+            {"name": "Beijing",          "center": [87, 383],  "population": 21500000, "region": "east_asia"},
+            {"name": "Mumbai",           "center": [107, 333], "population": 21700000, "region": "east_asia"},
+
+            {"name": "Jakarta",          "center": [120, 374], "population": 34500000, "region": "fill"},
+            {"name": "Dhaka",            "center": [103, 351], "population": 23000000, "region": "fill"},
+            {"name": "Manila",           "center": [109, 388], "population": 14400000, "region": "fill"},
+            {"name": "Bangkok",          "center": [109, 370], "population": 11000000, "region": "fill"},
+            {"name": "Seoul",            "center": [90, 393],  "population": 9800000,  "region": "fill"},
+            {"name": "Osaka",            "center": [93, 398],  "population": 19300000, "region": "fill"},
+            {"name": "Chongqing",        "center": [97, 375],  "population": 17000000, "region": "fill"},
+            {"name": "Kolkata",          "center": [103, 349], "population": 15100000, "region": "fill"},
+            {"name": "Lahore",           "center": [97, 336],  "population": 14000000, "region": "fill"},
+            {"name": "Shenzhen",         "center": [104, 382], "population": 13400000, "region": "fill"},
+            {"name": "Bangalore",        "center": [111, 339], "population": 13200000, "region": "fill"},
+            {"name": "Ho Chi Minh City", "center": [113, 374], "population": 9300000,  "region": "fill"},
+            {"name": "Luanda",           "center": [132, 268], "population": 9000000,  "region": "fill"},
+            {"name": "Addis Ababa",      "center": [119, 292], "population": 5500000,  "region": "fill"},
+            {"name": "Santiago",         "center": [147, 137], "population": 7000000,  "region": "fill"},
+            {"name": "Taipei",           "center": [103, 388], "population": 7000000,  "region": "fill"},
+            {"name": "Hong Kong",        "center": [104, 382], "population": 7500000,  "region": "fill"},
+            {"name": "Singapore",        "center": [119, 372], "population": 5900000,  "region": "fill"},
+            {"name": "Sydney",           "center": [161, 421], "population": 5300000,  "region": "fill"},
+            {"name": "Casablanca",       "center": [93, 249],  "population": 3800000,  "region": "fill"}
+        ]
+    }
+}

tooling/planet-gen/sol_markers/luna_features.json

@@ -0,0 +1,16 @@
+{
+    "GJ0d-1": {
+        "pois": [
+            {"id": "poi_mare_tranquillitatis",  "name": "Mare Tranquillitatis",  "center": [119, 282], "kind": "mare"},
+            {"id": "poi_mare_imbrium",           "name": "Mare Imbrium",           "center": [93, 247],  "kind": "mare"},
+            {"id": "poi_oceanus_procellarum",    "name": "Oceanus Procellarum",    "center": [107, 230], "kind": "mare"},
+            {"id": "poi_mare_serenitatis",       "name": "Mare Serenitatis",       "center": [104, 277], "kind": "mare"},
+            {"id": "poi_mare_crisium",           "name": "Mare Crisium",           "center": [108, 302], "kind": "mare"},
+            {"id": "poi_mare_nubium",            "name": "Mare Nubium",            "center": [134, 248], "kind": "mare"},
+            {"id": "poi_mare_fecunditatis",      "name": "Mare Fecunditatis",      "center": [124, 299], "kind": "mare"},
+            {"id": "poi_south_pole_aitken",      "name": "South Pole-Aitken Basin","center": [213, 330], "kind": "basin"},
+            {"id": "poi_tycho",                  "name": "Tycho",                  "center": [163, 249], "kind": "crater"},
+            {"id": "poi_copernicus",             "name": "Copernicus",             "center": [118, 243], "kind": "crater"}
+        ]
+    }
+}

tooling/planet-gen/sol_markers/mars_features.json

@@ -0,0 +1,23 @@
+{
+    "GJ0e": {
+        "mountain_ranges": [
+            {"name": "Olympus Mons",     "peak": [107, 358], "center": [107, 358]},
+            {"name": "Tharsis Bulge",    "peak": [115, 365], "center": [118, 362]},
+            {"name": "Elysium Mons",     "peak": [103, 413], "center": [103, 413]},
+            {"name": "Ascraeus Mons",    "peak": [108, 367], "center": [108, 367]},
+            {"name": "Arsia Mons",       "peak": [118, 363], "center": [118, 363]}
+        ],
+        "oceans": [
+            {"name": "Hellas Basin",    "center": [148, 329]},
+            {"name": "Utopia Planitia", "center": [80, 385]},
+            {"name": "Isidis Planitia", "center": [112, 343]}
+        ],
+        "pois": [
+            {"id": "poi_valles_marineris", "name": "Valles Marineris", "center": [118, 380], "kind": "canyon"},
+            {"id": "poi_north_polar_cap",  "name": "North Polar Cap",  "center": [10, 256],  "kind": "ice_cap"},
+            {"id": "poi_south_polar_cap",  "name": "South Polar Cap",  "center": [245, 256], "kind": "ice_cap"},
+            {"id": "poi_chryse_planitia",  "name": "Chryse Planitia",  "center": [100, 392], "kind": "plain"},
+            {"id": "poi_acidalia_planitia","name": "Acidalia Planitia","center": [80, 395],  "kind": "plain"}
+        ]
+    }
+}

tooling/planet-gen/sol_markers/outer_features.json

@@ -0,0 +1,48 @@
+{
+    "GJ0f-1": {
+        "pois": [
+            {"id": "poi_loki_patera",    "name": "Loki Patera",    "center": [115, 295], "kind": "volcano"},
+            {"id": "poi_pele",           "name": "Pele",           "center": [140, 358], "kind": "volcano"},
+            {"id": "poi_tvashtar",       "name": "Tvashtar Paterae","center": [46, 354],  "kind": "volcano"},
+            {"id": "poi_prometheus",     "name": "Prometheus",     "center": [128, 412], "kind": "volcano"},
+            {"id": "poi_masubi",         "name": "Masubi",         "center": [152, 370], "kind": "volcano"}
+        ]
+    },
+    "GJ0f-2": {
+        "pois": [
+            {"id": "poi_conamara_chaos", "name": "Conamara Chaos", "center": [118, 328], "kind": "chaos"},
+            {"id": "poi_pwyll_crater",   "name": "Pwyll Crater",   "center": [155, 328], "kind": "crater"},
+            {"id": "poi_thera_macula",   "name": "Thera Macula",   "center": [145, 340], "kind": "macula"},
+            {"id": "poi_tyre",           "name": "Tyre",           "center": [93, 357],  "kind": "multi_ring"}
+        ]
+    },
+    "GJ0f-3": {
+        "pois": [
+            {"id": "poi_galileo_regio",  "name": "Galileo Regio",  "center": [90, 375],  "kind": "dark_terrain"},
+            {"id": "poi_uruk_sulcus",    "name": "Uruk Sulcus",    "center": [108, 310], "kind": "grooved"},
+            {"id": "poi_gilgamesh",      "name": "Gilgamesh",      "center": [178, 370], "kind": "crater"}
+        ]
+    },
+    "GJ0f-4": {
+        "pois": [
+            {"id": "poi_valhalla",       "name": "Valhalla",       "center": [108, 310], "kind": "multi_ring"},
+            {"id": "poi_asgard",         "name": "Asgard",         "center": [93, 370],  "kind": "multi_ring"}
+        ]
+    },
+    "GJ0g-1": {
+        "pois": [
+            {"id": "poi_kraken_mare",    "name": "Kraken Mare",    "center": [25, 340],  "kind": "methane_sea"},
+            {"id": "poi_ligeia_mare",    "name": "Ligeia Mare",    "center": [20, 370],  "kind": "methane_sea"},
+            {"id": "poi_punga_mare",     "name": "Punga Mare",     "center": [30, 350],  "kind": "methane_sea"},
+            {"id": "poi_xanadu",         "name": "Xanadu",         "center": [117, 375], "kind": "bright_terrain"},
+            {"id": "poi_shangri_la",     "name": "Shangri-La",     "center": [130, 310], "kind": "dune_field"}
+        ]
+    },
+    "GJ0g-2": {
+        "pois": [
+            {"id": "poi_tiger_stripes",  "name": "Tiger Stripes",  "center": [220, 256], "kind": "fracture"},
+            {"id": "poi_baghdad_sulcus",  "name": "Baghdad Sulcus", "center": [218, 270], "kind": "fracture"},
+            {"id": "poi_samarkand_sulcus","name": "Samarkand Sulcus","center": [215, 240], "kind": "fracture"}
+        ]
+    }
+}

tooling/planet-gen/sol_overrides.json

@@ -0,0 +1,108 @@
+{
+    "GJ0b": {
+        "_comment": "Mercury — barren rock, tidally locked, extreme temps",
+        "orbit": { "axial_tilt_deg": 0.034 },
+        "terrain": { "land_fraction": 1.0, "tectonics": "none" },
+        "environment": { "geothermal_flux": "low", "substrate": "silicate" }
+    },
+    "GJ0c": {
+        "_comment": "Venus — thick sulfuric acid clouds hide the surface completely",
+        "orbit": { "axial_tilt_deg": 177.4 },
+        "terrain": { "land_fraction": 1.0, "tectonics": "active" },
+        "physical": { "atmosphere_color": [0.92, 0.85, 0.55] },
+        "environment": { "geothermal_flux": "high", "substrate": "silicate" },
+        "clouds": { "enabled": true, "coverage_base": 0.95 }
+    },
+    "GJ0d": {
+        "_comment": "Earth — use real-world data pipeline",
+        "orbit": { "axial_tilt_deg": 23.44 },
+        "terrain": { "land_fraction": 0.29, "polar_ice_lat": 0.85, "tectonics": "active" },
+        "environment": { "hydrosphere": "ocean", "geothermal_flux": "low" },
+        "clouds": { "enabled": true, "coverage_base": 0.5 }
+    },
+    "GJ0d-1": {
+        "_comment": "Luna — barren, tidally locked",
+        "orbit": { "axial_tilt_deg": 6.68 },
+        "terrain": { "land_fraction": 1.0, "tectonics": "none" },
+        "environment": { "geothermal_flux": "low", "substrate": "silicate" }
+    },
+    "GJ0e": {
+        "_comment": "Mars — thin atmo, partially terraformed in lore (800 years)",
+        "orbit": { "axial_tilt_deg": 25.19 },
+        "terrain": { "land_fraction": 0.98, "polar_ice_lat": 0.65, "tectonics": "none" },
+        "environment": { "hydrosphere": "ice", "geothermal_flux": "low", "substrate": "silicate" }
+    },
+    "GJ0f": {
+        "_comment": "Jupiter — gas giant, Great Red Spot",
+        "gas_giant": {
+            "band_palette": "jovian",
+            "storm_count": 2,
+            "storm_max_size": 0.12
+        },
+        "rings": false
+    },
+    "GJ0f-1": {
+        "_comment": "Io — volcanic moon of Jupiter, tidal heating",
+        "terrain": { "land_fraction": 1.0, "tectonics": "extreme" },
+        "environment": { "geothermal_flux": "high", "substrate": "silicate", "chemosynthetic": false }
+    },
+    "GJ0f-2": {
+        "_comment": "Europa — ice moon, subsurface ocean",
+        "terrain": { "land_fraction": 1.0, "tectonics": "low" },
+        "environment": { "geothermal_flux": "low", "substrate": "ice" }
+    },
+    "GJ0f-3": {
+        "_comment": "Ganymede — largest moon, ice/rock dichotomy",
+        "terrain": { "land_fraction": 1.0, "tectonics": "none" },
+        "environment": { "geothermal_flux": "low", "substrate": "ice" }
+    },
+    "GJ0f-4": {
+        "_comment": "Callisto — heavily cratered ice moon",
+        "terrain": { "land_fraction": 1.0, "tectonics": "none" },
+        "environment": { "geothermal_flux": "low", "substrate": "ice" }
+    },
+    "GJ0g": {
+        "_comment": "Saturn — gas giant with prominent ring system",
+        "gas_giant": {
+            "band_palette": "saturnian",
+            "storm_count": 1,
+            "storm_max_size": 0.06
+        },
+        "rings": {
+            "enabled": true,
+            "inner_radius_factor": 1.12,
+            "outer_radius_factor": 2.65,
+            "opacity_base": 0.68,
+            "ring_color": [0.88, 0.78, 0.55]
+        }
+    },
+    "GJ0g-1": {
+        "_comment": "Titan — dense atmosphere, methane cycle",
+        "terrain": { "land_fraction": 0.60, "tectonics": "low" },
+        "environment": { "hydrosphere": "rivers", "geothermal_flux": "low", "substrate": "ice" }
+    },
+    "GJ0g-2": {
+        "_comment": "Enceladus — small ice moon, geysers",
+        "terrain": { "land_fraction": 1.0, "tectonics": "low" },
+        "environment": { "geothermal_flux": "moderate", "substrate": "ice" }
+    },
+    "GJ0h": {
+        "_comment": "Uranus — ice giant, extreme axial tilt",
+        "orbit": { "axial_tilt_deg": 97.8 },
+        "gas_giant": {
+            "band_palette": "icy",
+            "storm_count": 1,
+            "storm_max_size": 0.04
+        },
+        "rings": false
+    },
+    "GJ0i": {
+        "_comment": "Neptune — ice giant, active storms",
+        "gas_giant": {
+            "band_palette": "neptunian",
+            "storm_count": 3,
+            "storm_max_size": 0.08
+        },
+        "rings": false
+    }
+}