feat(p1): city flood stages and zones from CMU Water Center's 2024 flood map

CMU Water Center's Chiang Mai flood hazard map (CMFlood) replaces the older
municipal map. It is built on the October 2024 flood and has the Ping leave
its banks above P.1 4.20 m, with five stages: 4.30, 4.50, 4.70, 5.00 and
5.30 m (the 2024 record). The old map's seven stages (3.70-4.60 m) said the
city floods from 3.70 m; CMU calls that map out of date.

- P1_FLOOD_STAGES: the five CMU stages, discharge on RID's 2024 rating.
  /forecast's P.1 `stages` now has five entries (same format). Thresholds
  are unchanged: 4.20 m danger is the bank top, 3.70 m warning half a metre
  below (as on CMU's gauge page). flood_config.json and the Rust parity
  fixtures regenerated (only the stages changed).
- flood-zones-p1.geojson: CMU's KML (L1..L5) as five disjoint zones, zone k
  floods first at stage k; scripts/build_flood_zones.py rebuilds it.
- Dashboard: CMU's colours and credit (map attribution); the stage grid follows
  the number of stages; zone lookup for markers honours holes. The city
  verdict uses statusOf's classes: flooding (over the banks) from 4.20 m,
  high water from 3.70 m. Zones show by themselves when the river is over its
  banks or stage 1 is >= 50 % likely. Alert, replay and verdict texts no
  longer call 3.70 m city flooding (EN + TH).
- Docs: DATA_SOURCES 5.1 (CMFlood: the map, flood marks, poles, city gauges),
  FLOOD_FORECASTING, DASHBOARD_REDESIGN, NOTIFICATIONS, README credits.

scripts/build_flood_zones.py

@@ -0,0 +1,151 @@
+# /// script
+# requires-python = ">=3.11"
+# dependencies = ["shapely>=2.0"]
+# ///
+"""Build src/static/flood-zones-p1.geojson from CMU Water Center's flood map.
+
+CMU's Chiang Mai flood hazard map (watercenter.cmu.ac.th, CMFlood, built from
+the October 2024 flood) publishes the flooded area at five P.1 levels as KML:
+/data/KML5Layer/L1.kml (P.1 = 4.30 m) to L5.kml (5.30 m, the 2024 record).
+Each file is the whole extent at that level. The dashboard wants one ring per
+stage instead: zone k is what floods first when P.1 reaches stage k's level
+(the extent at level k minus everything flooded at a lower level).
+
+    uv run scripts/build_flood_zones.py              # fetch the KML from CMU
+    uv run scripts/build_flood_zones.py --kml-dir D  # or read L1..L5.kml from D
+
+The levels must match P1_FLOOD_STAGES in src/ml/features.py.
+"""
+
+import argparse
+import json
+import re
+import urllib.request
+from pathlib import Path
+
+import shapely
+from shapely import make_valid
+from shapely.geometry import MultiPolygon, Polygon, mapping
+from shapely.ops import unary_union
+
+REPO = Path(__file__).resolve().parents[1]
+OUT = REPO / "src" / "static" / "flood-zones-p1.geojson"
+BASE = "https://watercenter.cmu.ac.th/data/KML5Layer"
+LEVELS = [4.30, 4.50, 4.70, 5.00, 5.30]  # L1..L5, P.1 gauge level in m
+# ~3 m at Chiang Mai: invisible at the zoom levels the map allows, and it
+# keeps the file small. Slivers below MIN_AREA (about 300 m²) come from the
+# layers' borders not lining up exactly and are dropped.
+SIMPLIFY_DEG = 3e-5
+MIN_AREA_DEG2 = 2.5e-8
+GRID = 1e-5  # coordinates snapped to 5 decimals (~1 m) with valid geometry
+SOURCE = (
+    "CMU Water Center (ศูนย์วิชาการสนับสนุนการบริหารจัดการน้ำ มหาวิทยาลัยเชียงใหม่), "
+    "Chiang Mai flood hazard map from the October 2024 flood, "
+    "https://watercenter.cmu.ac.th/chiangmai/cmflood/warning/floodmap"
+)
+
+
+def fetch(i: int, kml_dir: Path | None) -> str:
+    if kml_dir:
+        return (kml_dir / f"L{i}.kml").read_text(encoding="utf-8")
+    req = urllib.request.Request(
+        f"{BASE}/L{i}.kml", headers={"User-Agent": "Mozilla/5.0"}
+    )
+    with urllib.request.urlopen(req, timeout=60) as r:
+        return r.read().decode("utf-8")
+
+
+def ring(text: str) -> list:
+    pts = []
+    for tok in text.split():
+        x, y, *_ = tok.split(",")
+        pts.append((float(x), float(y)))
+    return pts
+
+
+def kml_extent(text: str):
+    polys = []
+    for body in re.findall(r"<Polygon>(.*?)</Polygon>", text, re.S):
+        outer = re.search(
+            r"<outerBoundaryIs>.*?<coordinates>(.*?)</coordinates>", body, re.S
+        )
+        inner = re.findall(
+            r"<innerBoundaryIs>.*?<coordinates>(.*?)</coordinates>", body, re.S
+        )
+        polys.append(
+            make_valid(Polygon(ring(outer.group(1)), [ring(s) for s in inner]))
+        )
+    return unary_union(polys)
+
+
+def polygons(geom) -> list:
+    if isinstance(geom, Polygon):
+        return [geom]
+    return [
+        g for g in getattr(geom, "geoms", []) if isinstance(g, (Polygon, MultiPolygon))
+    ]
+
+
+def clean(geom):
+    parts = []
+    for g in polygons(geom):
+        for p in polygons(g):
+            if p.area >= MIN_AREA_DEG2:
+                parts.append(p)
+    return unary_union(parts)
+
+
+def rounded(coords):
+    # the geometry is already on the grid; this only drops float noise
+    if isinstance(coords, (list, tuple)) and coords and isinstance(coords[0], float):
+        return [round(coords[0], 5), round(coords[1], 5)]
+    return [rounded(c) for c in coords]
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description=__doc__.splitlines()[0])
+    ap.add_argument("--kml-dir", type=Path, help="read L1..L5.kml from here")
+    args = ap.parse_args()
+
+    flooded = None
+    features = []
+    for i, level in enumerate(LEVELS, 1):
+        extent = kml_extent(fetch(i, args.kml_dir))
+        extent = shapely.set_precision(
+            make_valid(extent.simplify(SIMPLIFY_DEG, preserve_topology=True)), GRID
+        )
+        if flooded is None:
+            new, flooded = extent, extent
+        else:
+            new = shapely.difference(extent, flooded, grid_size=GRID)
+            flooded = shapely.union(flooded, extent, grid_size=GRID)
+        zone = shapely.set_precision(clean(new), GRID)
+        assert zone.is_valid, f"zone {i} is not a valid geometry"
+        geo = mapping(zone)
+        features.append(
+            {
+                "type": "Feature",
+                "properties": {"zone": i, "level": level},
+                "geometry": {
+                    "type": geo["type"],
+                    "coordinates": rounded(geo["coordinates"]),
+                },
+            }
+        )
+        print(f"zone {i} (P.1 {level:.2f} m): {len(polygons(zone))} polygons")
+    doc = {
+        "type": "FeatureCollection",
+        "source": SOURCE,
+        "note": "Zone k floods first when P.1 reaches its level; zones are disjoint.",
+        # highest zone first, so the zones that flood first are drawn on top
+        "features": features[::-1],
+    }
+    OUT.write_text(
+        json.dumps(doc, ensure_ascii=False, separators=(",", ":")) + "\n",
+        encoding="utf-8",
+    )
+    print(f"wrote {OUT.relative_to(REPO)}: {OUT.stat().st_size / 1024:.0f} KB")
+
+
+if __name__ == "__main__":
+    main()