170 lines
7.4 KiB
Python
170 lines
7.4 KiB
Python
#!/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()
|
|
|