settled-reach/tooling/domains/atlas/planet/sol_data/titan.py

"""
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 scipy.ndimage import gaussian_filter

from tooling.domains.atlas.planet.sol_data.download import ensure_cached
from tooling.domains.atlas.planet.sol_data.shared import (
    GRID_W, GRID_H,
    load_image_as_elevation, resample_to_grid, normalize_01,
    compute_hillshade, assemble_terrain,
)
from tooling.core import console

# 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)
        console.event(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:
        console.event(f"Titan mosaic unavailable ({e}), synthetic", level="warn")
        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."""
    from tooling.domains.atlas.planet.planet_simulation import compute_biome

    console.event("  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 tooling.domains.atlas.planet.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

    console.event(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)
    console.event(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