settled-reach/tooling/planet-gen/sol_data/luna.py

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