settled-reach/spikes/heightmap-pipeline/generate_kallast.py

#!/usr/bin/env python3
"""
Heightmap Spike: Kallast (GJ144d)

Produces one annotated heightmap for Kallast — Ran system's inner habitable
world. Wave 1 agricultural planet, breathable atmosphere, 0.95g, temperate.

Pipeline:
  1. Tectonic simulation (pyplatec) → elevation grid
  2. Hydraulic erosion (numpy/scipy) → soften ridges, carve valleys
  3. Climate pass → moisture/temperature from latitude + elevation
  4. Terrain classification → biome zones
  5. Hydrology → flow accumulation → river network
  6. Settlement placement → cities along rivers + fertile plains
  7. Road network → minimum spanning connections between major cities
  8. Annotated render → PNG export

Usage:
    python3 generate_kallast.py [--output path] [--seed N] [--fast]

Outputs:
    kallast_heightmap.png   — annotated world map
    kallast_terrain.npy     — raw terrain grid (numpy, for batch reuse)
    kallast_rivers.npy      — river network mask
    kallast_settlements.json — city coordinates and names
"""

import argparse
import json
import math
import random
import time
import sys

import numpy as np
from PIL import Image, ImageDraw, ImageFont

# ─────────────────────────────────────────────────────────────────────────────
# Planet parameters (from wiki/systems.db: GJ144d Kallast)
# ─────────────────────────────────────────────────────────────────────────────

PLANET = {
    "name": "Kallast",
    "system": "Ran (GJ 144)",
    "body_id": "GJ144d",
    "gravity": 0.95,        # g — affects tectonic force scaling
    "atmosphere": "breathable",
    "biome_summary": "temperate",
    "hydrosphere": "ocean",
    "population": 2_000_000_000,
    "settlement_wave": 1,
    "settlement_age_years": 580,
    "industrial": "Agricultural_Syndic",
    # Tectonic profile: Earth-like, high activity (mature world, well-eroded)
    "plate_count": 10,
    "land_fraction": 0.40,  # 40% land coverage — continental grain belt world
    # Note: SEA_LEVEL is computed dynamically from land_fraction after
    # tectonic simulation, since pyplatec produces a skewed distribution
    # that does not map linearly to target coverage percentages.
}

# ─────────────────────────────────────────────────────────────────────────────
# Grid settings
# ─────────────────────────────────────────────────────────────────────────────

GRID_W = 1024
GRID_H = 512
OUTPUT_W = 4096
OUTPUT_H = 2048

# ─────────────────────────────────────────────────────────────────────────────
# Color palette (Settled Reach visual grammar: muted, earthy, legible)
# ─────────────────────────────────────────────────────────────────────────────

PALETTE = {
    "ocean_deep":    (18,  32,  58),
    "ocean_mid":     (28,  52,  90),
    "ocean_shallow": (42,  80, 110),
    "coast":         (80, 105,  75),
    "lowland":       (95, 115,  65),
    "plains":        (130, 145,  80),
    "grassland":     (165, 155,  75),   # temperate plains / savanna (amber tone)
    "hills":         (120, 110,  80),
    "forest":        (60,  90,  55),
    "highland":      (110, 100,  90),
    "mountain":      (140, 130, 120),
    "peak":          (195, 190, 185),
    "snow":          (230, 228, 225),
    # Annotation colors
    "river":         (80, 140, 200),
    "road_major":    (180, 155,  90),
    "road_minor":    (160, 140,  85),
    "city_major":    (220,  60,  50),
    "city_minor":    (200, 110,  60),
    "city_label":    (240, 235, 220),
    "freight_elev":  (200, 180, 100),   # freight elevator pads
    "grid_line":     (255, 255, 255),
    "title_bg":      (15,  18,  25),
    "title_text":    (200, 208, 224),   # insert chrome: #c8d0e0
}


def step(msg: str):
    print(f"  [{time.strftime('%H:%M:%S')}] {msg}", flush=True)


# ─────────────────────────────────────────────────────────────────────────────
# 1. Tectonic simulation
# ─────────────────────────────────────────────────────────────────────────────

def run_tectonics(seed: int, fast: bool = False) -> np.ndarray:
    """Run pyplatec tectonic simulation. Returns normalised float32 grid."""
    step("Running tectonic simulation (pyplatec)…")
    import platec

    sim_steps = 50 if fast else 200

    p = platec.create(
        seed,
        GRID_W, GRID_H,
        sea_level=PLANET["land_fraction"],
        erosion_period=60,
        folding_ratio=0.02,
        aggr_overlap_abs=1_000_000,
        aggr_overlap_rel=0.33,
        cycle_count=2,
        num_plates=PLANET["plate_count"],
    )

    for _ in range(sim_steps):
        platec.step(p)

    hmap_raw = platec.get_heightmap(p)
    platec.destroy(p)

    arr = np.array(hmap_raw, dtype=np.float32).reshape(GRID_H, GRID_W)
    # Normalise to [0, 1]
    lo, hi = arr.min(), arr.max()
    arr = (arr - lo) / (hi - lo + 1e-9)
    step(f"Tectonic done. Elevation range: [{lo:.1f}, {hi:.1f}]")
    return arr


# ─────────────────────────────────────────────────────────────────────────────
# 2. Hydraulic erosion (simplified — scipy gaussian smoothing on steep slopes)
# ─────────────────────────────────────────────────────────────────────────────

def erode(terrain: np.ndarray, passes: int = 3) -> np.ndarray:
    """Simplified erosion: smooth with slope-weighted kernel."""
    step(f"Applying erosion ({passes} passes)…")
    from scipy.ndimage import gaussian_filter, uniform_filter

    result = terrain.copy()
    for i in range(passes):
        # Identify steep slopes
        gy, gx = np.gradient(result)
        slope = np.sqrt(gx**2 + gy**2)
        # Smooth strongly on steep areas (erosion), less on plains
        smooth = gaussian_filter(result, sigma=1.5)
        weight = np.clip(slope * 8, 0, 1)
        result = result * (1 - weight * 0.4) + smooth * (weight * 0.4)
    return result


# ─────────────────────────────────────────────────────────────────────────────
# 3. Terrain classification
# ─────────────────────────────────────────────────────────────────────────────

SEA_LEVEL = None  # set dynamically after tectonic simulation

def compute_sea_level(terrain: np.ndarray, land_fraction: float) -> float:
    """
    Compute sea level as the percentile that yields the target land fraction.
    pyplatec produces a skewed elevation distribution (most area at low elevation,
    peaks only at plate boundaries), so we cannot use a fixed fraction of the
    0-1 normalised range.
    """
    ocean_fraction = 1.0 - land_fraction
    sl = float(np.percentile(terrain, ocean_fraction * 100))
    step(f"Sea level computed: {sl:.4f} (target land={land_fraction*100:.0f}%, "
         f"actual≈{(terrain >= sl).sum() / terrain.size * 100:.1f}%)")
    return sl


def classify(terrain: np.ndarray) -> np.ndarray:
    """
    Returns integer class array relative to the dynamic sea level.
    Thresholds are expressed as fractions of the land elevation range
    (sea_level to max) rather than fixed offsets, so they scale correctly
    regardless of the pyplatec output distribution.

      0 = ocean_deep   (lowest)
      1 = ocean_mid
      2 = ocean_shallow
      3 = coast        (just above SL)
      4 = lowland
      5 = plains
      6 = grassland    (temperate plains / savanna — natural terrain class)
      7 = hills
      8 = forest
      9 = highland
     10 = mountain
     11 = peak
     12 = snow         (highest)
    """
    sl = SEA_LEVEL
    land_max = terrain.max()
    land_range = max(land_max - sl, 1e-6)

    c = np.zeros_like(terrain, dtype=np.int8)
    # Ocean bands (below sea level)
    ocean_range = max(sl - terrain.min(), 1e-6)
    c[terrain >= terrain.min()]          = 0   # ocean_deep (baseline)
    c[terrain >= sl - ocean_range * 0.5] = 1   # ocean_mid
    c[terrain >= sl - ocean_range * 0.2] = 2   # ocean_shallow
    # Land bands (above sea level, as fraction of land_range)
    c[terrain >= sl]                     = 3   # coast
    c[terrain >= sl + land_range * 0.05] = 4   # lowland
    c[terrain >= sl + land_range * 0.15] = 5   # plains
    c[terrain >= sl + land_range * 0.28] = 6   # grassland
    c[terrain >= sl + land_range * 0.42] = 7   # hills
    c[terrain >= sl + land_range * 0.53] = 8   # forest
    c[terrain >= sl + land_range * 0.63] = 9   # highland
    c[terrain >= sl + land_range * 0.74] = 10  # mountain
    c[terrain >= sl + land_range * 0.85] = 11  # peak
    c[terrain >= sl + land_range * 0.93] = 12  # snow
    return c


CLASS_COLORS = [
    PALETTE["ocean_deep"],
    PALETTE["ocean_mid"],
    PALETTE["ocean_shallow"],
    PALETTE["coast"],
    PALETTE["lowland"],
    PALETTE["plains"],
    PALETTE["grassland"],
    PALETTE["hills"],
    PALETTE["forest"],
    PALETTE["highland"],
    PALETTE["mountain"],
    PALETTE["peak"],
    PALETTE["snow"],
]


# ─────────────────────────────────────────────────────────────────────────────
# 4. Hydrology — flow accumulation → river network
# ─────────────────────────────────────────────────────────────────────────────

def compute_rivers(terrain: np.ndarray, threshold: int = 30) -> np.ndarray:
    """
    Simple D8 flow accumulation. Returns boolean mask of river cells.
    Not physically accurate but produces plausible branching networks.
    """
    step("Computing river network…")
    H, W = terrain.shape
    # D8 direction offsets
    dirs = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]

    # For each land cell, find steepest descent
    flow_acc = np.zeros((H, W), dtype=np.int32)
    land = terrain >= SEA_LEVEL

    # Simplified: accumulate flow by draining from high to low
    # Sort cells by elevation descending
    ys, xs = np.where(land)
    order = np.argsort(terrain[ys, xs])[::-1]
    ys_sorted = ys[order]
    xs_sorted = xs[order]

    for y, x in zip(ys_sorted, xs_sorted):
        flow_acc[y, x] += 1
        # Find steepest downslope neighbour
        best_drop = 0
        best_ny, best_nx = -1, -1
        for dy, dx in dirs:
            ny, nx = y + dy, x + dx
            if 0 <= ny < H and 0 <= nx < W:
                drop = terrain[y, x] - terrain[ny, nx]
                if drop > best_drop:
                    best_drop = drop
                    best_ny, best_nx = ny, nx
        if best_ny >= 0:
            flow_acc[best_ny, best_nx] += flow_acc[y, x]

    rivers = (flow_acc > threshold) & land
    step(f"River network: {rivers.sum()} cells above threshold {threshold}")
    return rivers


# ─────────────────────────────────────────────────────────────────────────────
# 5. Settlement placement
# ─────────────────────────────────────────────────────────────────────────────

KALLAST_CITIES = [
    # Name, (relative grid position), tier (1=capital, 2=regional, 3=node)
    # Placed on fertile plains near river confluences.
    # Temperate plains world → most cities on the grassland continental shelf.
    ("Kallast Prime",   (0.45, 0.48), 1),   # capital, central continent
    ("Ardenvall",       (0.30, 0.40), 2),   # western grain province
    ("Thessmark",       (0.60, 0.52), 2),   # eastern province
    ("Coldwater",       (0.25, 0.62), 2),   # southern coast, fishing + export
    ("Brightfield",     (0.50, 0.36), 2),   # northern plains
    ("Vorn's Crossing", (0.38, 0.55), 3),   # river crossing, freight node
    ("Saltmere",        (0.68, 0.42), 3),   # coast + processing node
    ("Kaspel",          (0.20, 0.50), 3),   # western interior node
    ("Drenmark",        (0.72, 0.58), 3),   # southeastern node
    ("New Farrow",      (0.55, 0.64), 3),   # southern freight hub
    ("Ossenfield",      (0.42, 0.30), 3),   # northern highland approach
    ("Tyne Station",    (0.33, 0.45), 3),   # freight elevator ground station
]

# Freight elevator locations (visible from orbit per wiki)
FREIGHT_ELEVATORS = [
    ("Kallast Anchor",  (0.45, 0.46)),
    ("Ardenvall Lift",  (0.29, 0.38)),
    ("Thessmark Riser", (0.61, 0.50)),
]


def place_settlements(terrain: np.ndarray, cities: list) -> list:
    """
    Snap city positions to nearest suitable terrain cell.
    Suitable = plains or grassland class, preferably near river.
    Returns list of (name, grid_y, grid_x, tier) tuples.
    """
    step("Placing settlements…")
    classified = classify(terrain)
    # Debug: show land cell distribution
    for cls_id in range(13):
        n = (classified == cls_id).sum()
        if n > 0:
            step(f"  class {cls_id}: {n} cells")
    H, W = terrain.shape
    placed = []

    for name, (rx, ry), tier in cities:
        cx = int(rx * W)
        cy = int(ry * H)
        # Search in expanding radius for valid terrain (up to 60 cells)
        best_y, best_x = cy, cx
        best_score = -1
        for radius in [20, 40, 60, 80]:
            for dy in range(-radius, radius + 1):
                for dx in range(-radius, radius + 1):
                    ty, tx = cy + dy, cx + dx
                    if 0 <= ty < H and 0 <= tx < W:
                        c = classified[ty, tx]
                        # Score: higher for plains/grain_belt, acceptable for coast/lowland
                        score = 0
                        if c in (5, 6):   # plains/grassland (ideal settlement terrain)
                            score = 20 - (abs(dy) + abs(dx)) * 0.2
                        elif c == 4:      # lowland
                            score = 12 - (abs(dy) + abs(dx)) * 0.2
                        elif c == 3:      # coast (ports, export hubs)
                            score = 8 - (abs(dy) + abs(dx)) * 0.2
                        elif c == 7:      # hills (defensible / highland cities)
                            score = 5 - (abs(dy) + abs(dx)) * 0.2
                        if score > best_score:
                            best_score = score
                            best_y, best_x = ty, tx
            if best_score > 0:
                break   # found valid terrain at this radius, stop expanding
        placed.append((name, best_y, best_x, tier))
        step(f"  {name} → ({best_x}, {best_y}) terrain={classified[best_y, best_x]}")

    return placed


# ─────────────────────────────────────────────────────────────────────────────
# 6. Road network (simple greedy connections)
# ─────────────────────────────────────────────────────────────────────────────

def build_roads(settlements: list) -> list:
    """
    Connect tier-1 and tier-2 cities with major roads.
    Connect tier-3 nodes to nearest tier-1 or tier-2.
    Returns list of (y1, x1, y2, x2, road_type) tuples.
    """
    roads = []
    tier12 = [(n, y, x) for (n, y, x, t) in settlements if t <= 2]
    tier3  = [(n, y, x) for (n, y, x, t) in settlements if t == 3]

    # Connect all tier-1/2 cities in order (simple chain + cross-links)
    for i in range(len(tier12) - 1):
        _, y1, x1 = tier12[i]
        _, y2, x2 = tier12[i + 1]
        roads.append((y1, x1, y2, x2, "major"))
    # Capital spurs to each tier-2
    cap = tier12[0]
    for city in tier12[1:]:
        roads.append((cap[1], cap[2], city[1], city[2], "major"))

    # Tier-3 nodes to nearest tier-1/2
    for (n3, y3, x3) in tier3:
        best_d = 1e9
        best = tier12[0]
        for city in tier12:
            d = (city[1]-y3)**2 + (city[2]-x3)**2
            if d < best_d:
                best_d = d
                best = city
        roads.append((y3, x3, best[1], best[2], "minor"))

    return roads


# ─────────────────────────────────────────────────────────────────────────────
# 7. Render
# ─────────────────────────────────────────────────────────────────────────────

def render(
    terrain: np.ndarray,
    rivers: np.ndarray,
    output_path: str,
):
    """
    Renders a geographic heightmap: terrain colors, rivers, coastlines,
    biome zones, lat/lon grid, title panel.
    Human layer (settlements, roads, irrigation) is stored in the JSON sidecar only.

    River rendering order matters: dilate at source resolution to preserve flow
    topology, then paint AFTER upscaling via NEAREST neighbor to avoid LANCZOS
    blurring the river network into invisibility.
    """
    step("Rendering geographic heightmap…")
    from scipy.ndimage import binary_dilation
    from PIL import ImageFont
    H, W = terrain.shape
    classified = classify(terrain)

    # UI scale relative to 1024×512 reference resolution
    UI_SCALE = OUTPUT_W / 1024

    # Base terrain color layer
    rgb = np.zeros((H, W, 3), dtype=np.uint8)
    for cls_id, color in enumerate(CLASS_COLORS):
        mask = classified == cls_id
        rgb[mask] = color

    # Slight elevation shading on land (ambient occlusion approximation)
    land_mask = terrain >= SEA_LEVEL
    elev_norm = np.clip((terrain - SEA_LEVEL) / (1.0 - SEA_LEVEL + 1e-9), 0, 1)
    shade = (0.85 + 0.15 * elev_norm)[..., np.newaxis]
    rgb = np.where(land_mask[..., np.newaxis], (rgb * shade).astype(np.uint8), rgb)

    # Dilate river mask at source resolution (preserves branching topology)
    rivers_drawn = binary_dilation(rivers, iterations=2)
    step(f"River cells after dilation: {rivers_drawn.sum()}")

    # Upscale terrain with LANCZOS (smooth gradients) — WITHOUT rivers painted in yet
    img = Image.fromarray(rgb).resize((OUTPUT_W, OUTPUT_H), Image.LANCZOS)

    # Paint rivers AFTER upscale using NEAREST neighbor — no blur, sharp edges
    rivers_up = Image.fromarray((rivers_drawn.astype(np.uint8) * 255)).resize(
        (OUTPUT_W, OUTPUT_H), Image.NEAREST
    )
    img_arr = np.array(img)
    img_arr[np.array(rivers_up) > 0] = PALETTE["river"]
    img = Image.fromarray(img_arr)
    draw = ImageDraw.Draw(img)

    # Lat/lon grid lines (every 30°)
    grid_w = max(1, int(UI_SCALE))
    for lat_pct in [1/6, 2/6, 3/6, 4/6, 5/6]:
        y = int(lat_pct * OUTPUT_H)
        draw.line([(0, y), (OUTPUT_W, y)], fill=(80, 90, 100), width=grid_w)
    for lon_pct in [1/6, 2/6, 3/6, 4/6, 5/6]:
        x = int(lon_pct * OUTPUT_W)
        draw.line([(x, 0), (x, OUTPUT_H)], fill=(80, 90, 100), width=grid_w)

    # Fonts — load_default(size=N) requires Pillow 10+
    try:
        font_title = ImageFont.load_default(size=int(14 * UI_SCALE))
        font_sub   = ImageFont.load_default(size=int(12 * UI_SCALE))
        font_small = ImageFont.load_default(size=int(11 * UI_SCALE))
    except TypeError:
        font_title = font_sub = font_small = ImageFont.load_default()

    # Title / metadata panel (insert chrome aesthetic)
    panel_h = int(56 * UI_SCALE)
    panel = Image.new('RGBA', (OUTPUT_W, panel_h), (15, 18, 25, 200))
    img = img.convert('RGBA')
    img.paste(panel, (0, 0), panel)
    img = img.convert('RGB')
    draw = ImageDraw.Draw(img)

    px = int(12 * UI_SCALE)
    py = int(8 * UI_SCALE)
    lh = int(18 * UI_SCALE)
    draw.text((px, py),        f"KALLAST  ·  GJ 144d  ·  Ran System", fill=PALETTE["title_text"], font=font_title)
    draw.text((px, py + lh),   f"temperate / ocean / breathable / 0.95g / {PLANET['population']//1_000_000_000:.1f}B pop / Wave {PLANET['settlement_wave']} / {PLANET['settlement_age_years']}y settled", fill=(130, 140, 160), font=font_sub)
    draw.text((px, py + lh*2), "HEIGHTMAP SPIKE v0.1 — Settled Reach Phase 1", fill=(80, 90, 110), font=font_small)

    # Legend (bottom strip) — natural geographic features only
    sw = int(16 * UI_SCALE)
    legend_y = OUTPUT_H - int(40 * UI_SCALE)
    legend_items = [
        ("ocean",          PALETTE["ocean_deep"]),
        ("coast",          PALETTE["coast"]),
        ("plains",         PALETTE["plains"]),
        ("grassland",      PALETTE["grassland"]),
        ("mountain",       PALETTE["mountain"]),
        ("river",          PALETTE["river"]),
    ]
    lx = px
    for label, color in legend_items:
        draw.rectangle([(lx, legend_y + int(4*UI_SCALE)), (lx+sw, legend_y + sw + int(4*UI_SCALE))], fill=color)
        draw.text((lx + sw + int(4*UI_SCALE), legend_y + int(3*UI_SCALE)), label, fill=(180, 185, 200), font=font_small)
        lx += int(110 * UI_SCALE)

    img.save(output_path, format="PNG", optimize=False)
    step(f"Saved: {output_path}")


# ─────────────────────────────────────────────────────────────────────────────
# Main
# ─────────────────────────────────────────────────────────────────────────────

def main():
    parser = argparse.ArgumentParser(description="Generate Kallast heightmap spike")
    parser.add_argument("--output", default="spikes/heightmap-pipeline/kallast_heightmap.png")
    parser.add_argument("--terrain-out", default="spikes/heightmap-pipeline/kallast_terrain.npy")
    parser.add_argument("--seed", type=int, default=144042)  # GJ144d seed
    parser.add_argument("--fast", action="store_true", help="Fewer tectonic steps (quicker, less detail)")
    args = parser.parse_args()

    print(f"\nKallast Heightmap Spike")
    print(f"  Planet: {PLANET['name']} ({PLANET['body_id']})")
    print(f"  Grid:   {GRID_W}×{GRID_H}")
    print(f"  Seed:   {args.seed}")
    print(f"  Mode:   {'fast' if args.fast else 'standard'}\n")

    t0 = time.time()

    terrain = run_tectonics(args.seed, fast=args.fast)
    terrain = erode(terrain, passes=4)

    # Set dynamic sea level based on target land coverage
    global SEA_LEVEL
    SEA_LEVEL = compute_sea_level(terrain, PLANET["land_fraction"])

    rivers = compute_rivers(terrain, threshold=60)  # calibrated for 1024×512 grid; flat world needs lower relative threshold

    settlements = place_settlements(terrain, KALLAST_CITIES)

    render(terrain, rivers, args.output)

    np.save(args.terrain_out, terrain)
    step(f"Terrain grid saved: {args.terrain_out}")

    # Write settlements JSON
    sj_path = args.output.replace(".png", "_settlements.json")
    with open(sj_path, "w") as f:
        json.dump([
            {"name": n, "grid_y": int(gy), "grid_x": int(gx), "tier": t}
            for (n, gy, gx, t) in settlements
        ], f, indent=2)
    step(f"Settlement data saved: {sj_path}")

    elapsed = time.time() - t0
    print(f"\nDone in {elapsed:.1f}s")
    print(f"Output: {args.output}")


if __name__ == "__main__":
    main()