#!/usr/bin/env python3
"""Analyze first-pass landforms for the Ground Zero tile."""

from __future__ import annotations

import json
from datetime import datetime, timezone
from pathlib import Path

import numpy as np
from osgeo import gdal

gdal.UseExceptions()


PROJECT_ROOT = Path(__file__).resolve().parents[1]
TILE_ID = "gz_us_ca_pacifica_utm10n_e544_n4160"
DEM_PATH = PROJECT_ROOT / "Data" / "Terrain" / "Extracted" / TILE_ID / f"{TILE_ID}_1m_dem_subset.tif"
WATER_ANALYSIS_PATH = PROJECT_ROOT / "Data" / "Terrain" / "Analysis" / TILE_ID / f"{TILE_ID}_water_shoreline_analysis.json"
OUTPUT_DIR = PROJECT_ROOT / "Data" / "Terrain" / "Analysis" / TILE_ID
OUTPUT_PATH = OUTPUT_DIR / f"{TILE_ID}_landform_analysis.json"

PIXEL_SIZE_M = 1.0


def percent(count: int, total: int) -> float:
    return count / total * 100.0 if total else 0.0


def finite_stats(values: np.ndarray) -> dict[str, float | int | None]:
    finite = values[np.isfinite(values)]
    if finite.size == 0:
        return {
            "sample_count": 0,
            "min": None,
            "max": None,
            "mean": None,
            "p05": None,
            "p50": None,
            "p95": None,
        }

    return {
        "sample_count": int(finite.size),
        "min": float(np.min(finite)),
        "max": float(np.max(finite)),
        "mean": float(np.mean(finite)),
        "p05": float(np.percentile(finite, 5)),
        "p50": float(np.percentile(finite, 50)),
        "p95": float(np.percentile(finite, 95)),
    }


def summarize_mask(mask: np.ndarray, total: int) -> dict[str, float | int]:
    count = int(np.count_nonzero(mask))
    return {
        "samples": count,
        "percent": percent(count, total),
    }


def summarize_edge(values: np.ndarray, slope_degrees: np.ndarray) -> dict[str, object]:
    return {
        "elevation_m": finite_stats(values),
        "slope_degrees": finite_stats(slope_degrees),
    }


def main() -> None:
    dataset = gdal.Open(str(DEM_PATH), gdal.GA_ReadOnly)
    if dataset is None:
        raise RuntimeError(f"Could not open DEM: {DEM_PATH}")

    band = dataset.GetRasterBand(1)
    elevation = band.ReadAsArray().astype(float)
    nodata = band.GetNoDataValue()
    if nodata is not None:
        elevation[elevation == nodata] = np.nan

    finite = elevation[np.isfinite(elevation)]
    if finite.size == 0:
        raise RuntimeError(f"DEM has no finite elevation samples: {DEM_PATH}")

    dz_dy, dz_dx = np.gradient(elevation, PIXEL_SIZE_M, PIXEL_SIZE_M)
    slope_degrees = np.degrees(np.arctan(np.sqrt((dz_dx * dz_dx) + (dz_dy * dz_dy))))
    total_samples = int(finite.size)

    elevation_bands = {
        "low_coastal_floor_0_to_5m": summarize_mask((elevation >= 0.0) & (elevation <= 5.0), total_samples),
        "valley_floor_5_to_15m": summarize_mask((elevation > 5.0) & (elevation <= 15.0), total_samples),
        "lower_slope_15_to_30m": summarize_mask((elevation > 15.0) & (elevation <= 30.0), total_samples),
        "hillside_30_to_60m": summarize_mask((elevation > 30.0) & (elevation <= 60.0), total_samples),
        "upper_hillside_over_60m": summarize_mask(elevation > 60.0, total_samples),
    }

    slope_classes = {
        "flat_0_to_3deg": summarize_mask((slope_degrees >= 0.0) & (slope_degrees <= 3.0), total_samples),
        "gentle_3_to_8deg": summarize_mask((slope_degrees > 3.0) & (slope_degrees <= 8.0), total_samples),
        "moderate_8_to_15deg": summarize_mask((slope_degrees > 8.0) & (slope_degrees <= 15.0), total_samples),
        "steep_15_to_30deg": summarize_mask((slope_degrees > 15.0) & (slope_degrees <= 30.0), total_samples),
        "very_steep_over_30deg": summarize_mask(slope_degrees > 30.0, total_samples),
    }

    playable_terrain = np.isfinite(elevation) & (slope_degrees <= 15.0)
    difficult_terrain = np.isfinite(elevation) & (slope_degrees > 15.0)
    likely_drainage_corridor = np.isfinite(elevation) & (elevation <= 8.0) & (slope_degrees <= 8.0)

    edge_width = 10
    edges = {
        "north": summarize_edge(elevation[:edge_width, :].reshape(-1), slope_degrees[:edge_width, :].reshape(-1)),
        "south": summarize_edge(elevation[-edge_width:, :].reshape(-1), slope_degrees[-edge_width:, :].reshape(-1)),
        "west": summarize_edge(elevation[:, :edge_width].reshape(-1), slope_degrees[:, :edge_width].reshape(-1)),
        "east": summarize_edge(elevation[:, -edge_width:].reshape(-1), slope_degrees[:, -edge_width:].reshape(-1)),
    }

    water_analysis = json.loads(WATER_ANALYSIS_PATH.read_text(encoding="utf-8")) if WATER_ANALYSIS_PATH.exists() else {}
    contains_coastline = bool(water_analysis.get("classification", {}).get("contains_shoreline", False))

    result = {
        "schema_version": 1,
        "tile_id": TILE_ID,
        "generated_at_utc": datetime.now(timezone.utc).replace(microsecond=0).isoformat().replace("+00:00", "Z"),
        "source_dem": str(DEM_PATH.relative_to(PROJECT_ROOT)),
        "source_pixel_size_m": PIXEL_SIZE_M,
        "elevation_summary_m": finite_stats(elevation),
        "slope_summary_degrees": finite_stats(slope_degrees),
        "elevation_bands": elevation_bands,
        "slope_classes": slope_classes,
        "edge_summary": edges,
        "landform_classification": {
            "contains_mountains": False,
            "contains_hills": bool(np.nanmax(elevation) >= 30.0),
            "contains_steep_slopes": bool(np.nanmax(slope_degrees) > 30.0),
            "contains_river": False,
            "contains_stream": False,
            "contains_stream_candidate_from_dem": bool(np.count_nonzero(likely_drainage_corridor) > 0),
            "contains_lake": False,
            "contains_coastline": contains_coastline,
            "coastline_absence_confirmed_by_water_pass": not contains_coastline,
        },
        "gameplay_masks": {
            "walkable_or_buildable_first_pass": summarize_mask(playable_terrain, total_samples),
            "difficult_or_slow_travel_first_pass": summarize_mask(difficult_terrain, total_samples),
            "possible_drainage_or_freshwater_search_zone": summarize_mask(likely_drainage_corridor, total_samples),
        },
        "handling_decisions": {
            "hills": "Present. Use slope/elevation classes for movement modifiers, resource placement, and visual foliage density.",
            "mountains": "Not present in this 1 km MVP tile. Defer mountain rules to higher-elevation future tiles.",
            "river": "Not confirmed from current DEM-only pass. Require NHD or equivalent hydrography before placing real river geometry.",
            "stream": "DEM suggests low drainage/search zones, but stream presence is not confirmed. Use these zones as candidates for later NHD validation and MVP freshwater placement.",
            "lake": "Not present in current tile.",
            "coastline": "Absent inside the current tile. Coastline and bathymetry remain neighbor-tile work.",
        },
        "notes": [
            "This is a DEM-only first pass and intentionally does not invent rivers, lakes, or coastline.",
            "Use confirmed hydrography before committing permanent watercourse geometry.",
            "The tile has enough hills and slope variation to support movement-cost and resource-placement rules.",
        ],
    }

    OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
    OUTPUT_PATH.write_text(json.dumps(result, indent=2) + "\n", encoding="utf-8")
    print(f"Wrote {OUTPUT_PATH.relative_to(PROJECT_ROOT)}")
    print(json.dumps(result["landform_classification"], indent=2))


if __name__ == "__main__":
    main()

