#!/usr/bin/env python3 """Analyze first-pass water and shoreline applicability 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" OUTPUT_DIR = PROJECT_ROOT / "Data" / "Terrain" / "Analysis" / TILE_ID OUTPUT_PATH = OUTPUT_DIR / f"{TILE_ID}_water_shoreline_analysis.json" SEA_LEVEL_M = 0.0 SHORELINE_BUFFER_M = 2.0 LOW_COASTAL_BUFFER_M = 5.0 def summarize_edge(values: np.ndarray) -> dict[str, float | int]: finite = values[np.isfinite(values)] return { "sample_count": int(finite.size), "min_elevation_m": float(np.min(finite)) if finite.size else None, "max_elevation_m": float(np.max(finite)) if finite.size else None, "mean_elevation_m": float(np.mean(finite)) if finite.size else None, "sea_level_or_below_samples": int(np.count_nonzero(finite <= SEA_LEVEL_M)), "near_sea_level_samples": int(np.count_nonzero(finite <= SHORELINE_BUFFER_M)), "low_coastal_samples": int(np.count_nonzero(finite <= LOW_COASTAL_BUFFER_M)), } 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) data = band.ReadAsArray().astype(float) nodata = band.GetNoDataValue() if nodata is not None: data[data == nodata] = np.nan finite = data[np.isfinite(data)] if finite.size == 0: raise RuntimeError(f"DEM has no finite elevation samples: {DEM_PATH}") total_samples = int(finite.size) sea_or_below = int(np.count_nonzero(finite <= SEA_LEVEL_M)) near_sea = int(np.count_nonzero(finite <= SHORELINE_BUFFER_M)) low_coastal = int(np.count_nonzero(finite <= LOW_COASTAL_BUFFER_M)) edge_width = 10 edges = { "north": summarize_edge(data[:edge_width, :].reshape(-1)), "south": summarize_edge(data[-edge_width:, :].reshape(-1)), "west": summarize_edge(data[:, :edge_width].reshape(-1)), "east": summarize_edge(data[:, -edge_width:].reshape(-1)), } has_open_water = sea_or_below > 0 has_shoreline = near_sea > 0 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)), "thresholds": { "sea_level_m": SEA_LEVEL_M, "shoreline_buffer_m": SHORELINE_BUFFER_M, "low_coastal_buffer_m": LOW_COASTAL_BUFFER_M, }, "elevation_summary": { "sample_count": total_samples, "min_elevation_m": float(np.min(finite)), "max_elevation_m": float(np.max(finite)), "mean_elevation_m": float(np.mean(finite)), "sea_level_or_below_samples": sea_or_below, "near_sea_level_samples": near_sea, "low_coastal_samples": low_coastal, "sea_level_or_below_percent": sea_or_below / total_samples * 100.0, "near_sea_level_percent": near_sea / total_samples * 100.0, "low_coastal_percent": low_coastal / total_samples * 100.0, }, "edge_summary": edges, "classification": { "contains_open_water": has_open_water, "contains_shoreline": has_shoreline, "water_depth_required_for_mvp_tile": False, "shoreline_mesh_required_for_mvp_tile": False, "freshwater_source_required_for_gameplay": True, "recommended_first_pass": "No ocean/shoreline actor required inside the current Ground Zero tile. Add a small freshwater source for gameplay separately, and defer ocean/bathymetry work to west/southwest neighbor coastal tiles.", }, "notes": [ "This analysis uses the extracted 1-meter USGS DEM for the selected 1 km tile.", "The selected Ground Zero tile is coastal-influenced but the DEM minimum is above sea level.", "Ocean depth/bathymetry remains required for adjacent coastal/ocean tiles before multi-tile expansion.", ], } 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["classification"], indent=2)) if __name__ == "__main__": main()