{
  "subcatchments": [
    {
      "subcatchment_id": "D1",
      "dem_slope_pct": 2.8,
      "dem_flow_length_m": 62.0
    },
    {
      "subcatchment_id": "D2",
      "dem_elev_max_m": 108.0,
      "dem_elev_min_m": 101.0,
      "dem_flow_length_m": 75.0
    },
    {
      "subcatchment_id": "D3",
      "dem_elev_mean_m": 104.5,
      "dem_elev_outlet_m": 102.5,
      "dem_flow_length_m": 85.0
    }
  ]
}

#!/usr/bin/env python3
from __future__ import annotations

import argparse
import csv
import json
import os
from pathlib import Path
from typing import Any

import geopandas as gpd


REPO_ROOT = Path(__file__).resolve().parents[3]

# Sibling-skill seam (issue #246): this script defaults its landuse/soil
# lookup CSVs from the swmm-params skill's reference data. Resolve the
# skills/ root through --skills-root > AISWMM_SKILLS_ROOT env var > the
# repo-relative default, so a relocated/standalone deployment can point
# elsewhere without changing default behavior.
DEFAULT_SKILLS_ROOT = REPO_ROOT / "skills"
SKILLS_ROOT_ENV = "AISWMM_SKILLS_ROOT"
NUMERIC_LANDUSE_FIELDS = [
    "imperv_pct",
    "n_imperv",
    "n_perv",
    "dstore_imperv_in",
    "dstore_perv_in",
    "zero_imperv_pct",
    "pct_routed",
]
NUMERIC_SOIL_FIELDS = ["suction_mm", "ksat_mm_per_hr", "imdmax"]


def resolve_skills_root(cli_value: Path | None) -> Path:
    """Resolve the skills/ root directory holding sibling skills (e.g. swmm-params).

    Precedence: --skills-root flag > AISWMM_SKILLS_ROOT env var > default
    (repo-relative skills/ next to this skill's checkout).
    """
    if cli_value is not None:
        return cli_value
    env_value = os.environ.get(SKILLS_ROOT_ENV)
    if env_value:
        return Path(env_value)
    return DEFAULT_SKILLS_ROOT


def normalize_key(value: Any) -> str:
    return " ".join(str(value).strip().lower().split())


def load_csv_rows(path: Path) -> list[dict[str, str]]:
    with path.open(newline="", encoding="utf-8") as f:
        rows = list(csv.DictReader(f))
    if not rows:
        raise ValueError(f"CSV has no data rows: {path}")
    return rows


def parse_float(value: Any, *, field: str, csv_path: Path, row_number: int) -> float:
    if value is None or str(value).strip() == "":
        raise ValueError(f"Missing numeric value for '{field}' at {csv_path}:{row_number}")
    try:
        return float(str(value).strip())
    except ValueError as exc:
        raise ValueError(f"Invalid float for '{field}' at {csv_path}:{row_number}: {value}") from exc


def write_json(path: Path, obj: Any) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(json.dumps(obj, indent=2), encoding="utf-8")


def write_csv(path: Path, rows: list[dict[str, Any]], fieldnames: list[str]) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", newline="", encoding="utf-8") as f:
        writer = csv.DictWriter(f, fieldnames=fieldnames)
        writer.writeheader()
        for row in rows:
            writer.writerow({field: row.get(field, "") for field in fieldnames})


def load_landuse_lookup(path: Path) -> tuple[dict[str, dict[str, Any]], dict[str, Any] | None]:
    lookup: dict[str, dict[str, Any]] = {}
    default: dict[str, Any] | None = None
    for i, row in enumerate(load_csv_rows(path), start=2):
        raw = (row.get("landuse_class") or "").strip()
        if not raw:
            raise ValueError(f"Missing 'landuse_class' in lookup at {path}:{i}")
        rec = {
            "landuse_class": raw,
            "route_to": (row.get("route_to") or "").strip(),
            "notes": (row.get("notes") or "").strip(),
        }
        for field in NUMERIC_LANDUSE_FIELDS:
            rec[field] = parse_float(row.get(field), field=field, csv_path=path, row_number=i)
        lookup[normalize_key(raw)] = rec
        if raw.upper() == "DEFAULT":
            default = rec
    return lookup, default


def load_soil_lookup(path: Path) -> tuple[dict[str, dict[str, Any]], dict[str, Any] | None]:
    lookup: dict[str, dict[str, Any]] = {}
    default: dict[str, Any] | None = None
    for i, row in enumerate(load_csv_rows(path), start=2):
        raw = (row.get("texture") or "").strip()
        if not raw:
            raise ValueError(f"Missing 'texture' in lookup at {path}:{i}")
        rec = {"texture": raw, "notes": (row.get("notes") or "").strip()}
        for field in NUMERIC_SOIL_FIELDS:
            rec[field] = parse_float(row.get(field), field=field, csv_path=path, row_number=i)
        lookup[normalize_key(raw)] = rec
        if raw in {"-", "DEFAULT", "default"}:
            default = rec
    if default is None:
        default = lookup.get(normalize_key("-")) or lookup.get(normalize_key("default"))
    return lookup, default


def read_vector(path: Path, *, layer_name: str) -> gpd.GeoDataFrame:
    gdf = gpd.read_file(path)
    if gdf.empty:
        raise ValueError(f"{layer_name} layer has no features: {path}")
    if gdf.crs is None:
        raise ValueError(f"{layer_name} layer has no CRS: {path}")
    gdf = gdf[gdf.geometry.notna()].copy()
    gdf["geometry"] = gdf.geometry.buffer(0)
    gdf = gdf[~gdf.geometry.is_empty].copy()
    if gdf.empty:
        raise ValueError(f"{layer_name} layer has no valid polygon geometry after cleanup: {path}")
    return gdf


def ensure_projected(gdf: gpd.GeoDataFrame, *, layer_name: str) -> None:
    if not gdf.crs or not gdf.crs.is_projected:
        raise ValueError(f"{layer_name} CRS must be projected for area weighting, got {gdf.crs}")


def class_area_fractions(
    *,
    subcatchments: gpd.GeoDataFrame,
    thematic: gpd.GeoDataFrame,
    id_field: str,
    class_field: str,
    default_class: str,
    strict: bool,
    label: str,
) -> tuple[list[dict[str, Any]], list[dict[str, Any]]]:
    if class_field not in thematic.columns:
        raise ValueError(f"{label} field '{class_field}' not found. Available fields: {list(thematic.columns)}")

    sub = subcatchments[[id_field, "geometry"]].copy()
    sub["sub_area_m2"] = sub.geometry.area
    thematic = thematic[[class_field, "geometry"]].copy()
    if thematic.crs != sub.crs:
        thematic = thematic.to_crs(sub.crs)

    inter = gpd.overlay(sub, thematic, how="intersection", keep_geom_type=False)
    rows: list[dict[str, Any]] = []
    issues: list[dict[str, Any]] = []
    grouped: dict[str, dict[str, float]] = {}
    sub_areas = {str(row[id_field]): float(row["sub_area_m2"]) for _, row in sub.iterrows()}

    if not inter.empty:
        inter["intersect_area_m2"] = inter.geometry.area
        for _, row in inter.iterrows():
            sid = str(row[id_field])
            raw_class = str(row.get(class_field) or default_class).strip() or default_class
            area = float(row["intersect_area_m2"])
            if area <= 0:
                continue
            grouped.setdefault(sid, {})
            grouped[sid][raw_class] = grouped[sid].get(raw_class, 0.0) + area

    for sid, sub_area in sub_areas.items():
        if sub_area <= 0:
            raise ValueError(f"Subcatchment '{sid}' has non-positive area")
        class_areas = grouped.get(sid, {})
        covered = sum(class_areas.values())
        if covered <= 0:
            if strict:
                raise ValueError(f"Subcatchment '{sid}' has no {label} overlay coverage")
            class_areas = {default_class: sub_area}
            covered = sub_area
            issues.append({"severity": "warning", "id": sid, "message": f"No {label} coverage; used {default_class}."})
        elif covered < sub_area * 0.999:
            remainder = sub_area - covered
            class_areas[default_class] = class_areas.get(default_class, 0.0) + remainder
            issues.append(
                {
                    "severity": "warning",
                    "id": sid,
                    "message": f"{label} coverage below subcatchment area; assigned uncovered area to {default_class}.",
                    "coverage_fraction": covered / sub_area,
                }
            )
        elif covered > sub_area * 1.001:
            issues.append(
                {
                    "severity": "warning",
                    "id": sid,
                    "message": f"{label} overlay areas exceed subcatchment area; normalized fractions.",
                    "coverage_fraction": covered / sub_area,
                }
            )

        denominator = sum(class_areas.values())
        for class_name, area in sorted(class_areas.items()):
            fraction = area / denominator if denominator > 0 else 0.0
            if fraction <= 0:
                continue
            rows.append(
                {
                    "subcatchment_id": sid,
                    "class": class_name,
                    "area_m2": area,
                    "fraction": fraction,
                }
            )

    return rows, issues


def weighted_landuse(
    rows: list[dict[str, Any]],
    lookup: dict[str, dict[str, Any]],
    default: dict[str, Any] | None,
    *,
    strict: bool,
) -> tuple[dict[str, Any], list[dict[str, Any]], set[str]]:
    by_id: dict[str, list[dict[str, Any]]] = {}
    for row in rows:
        by_id.setdefault(row["subcatchment_id"], []).append(row)

    records: list[dict[str, Any]] = []
    subcatchment_section: list[dict[str, Any]] = []
    subarea_section: list[dict[str, Any]] = []
    audit: list[dict[str, Any]] = []
    unmatched: set[str] = set()

    for sid in sorted(by_id):
        weighted = {field: 0.0 for field in NUMERIC_LANDUSE_FIELDS}
        route_votes: dict[str, float] = {}
        components = []
        for row in by_id[sid]:
            class_name = row["class"]
            fraction = float(row["fraction"])
            rec = lookup.get(normalize_key(class_name))
            used_default = False
            if rec is None:
                if strict or default is None:
                    raise ValueError(f"Unmapped landuse class '{class_name}' for subcatchment '{sid}'")
                rec = default
                used_default = True
                unmatched.add(class_name)
            for field in NUMERIC_LANDUSE_FIELDS:
                weighted[field] += rec[field] * fraction
            route_to = rec["route_to"] or "OUTLET"
            route_votes[route_to] = route_votes.get(route_to, 0.0) + fraction
            components.append(
                {
                    "class": class_name,
                    "lookup_class": rec["landuse_class"],
                    "fraction": fraction,
                    "area_m2": row["area_m2"],
                    "used_default": used_default,
                }
            )
            audit.append(
                {
                    "subcatchment_id": sid,
                    "landuse_class": class_name,
                    "lookup_landuse_class": rec["landuse_class"],
                    "area_m2": round(float(row["area_m2"]), 6),
                    "fraction": round(fraction, 8),
                    "used_default": used_default,
                }
            )
        route_to = max(route_votes.items(), key=lambda item: item[1])[0]
        subcatchment_entry = {"id": sid, "pct_imperv": round(weighted["imperv_pct"], 6)}
        subarea_entry = {
            "id": sid,
            "n_imperv": round(weighted["n_imperv"], 6),
            "n_perv": round(weighted["n_perv"], 6),
            "dstore_imperv_in": round(weighted["dstore_imperv_in"], 6),
            "dstore_perv_in": round(weighted["dstore_perv_in"], 6),
            "zero_imperv_pct": round(weighted["zero_imperv_pct"], 6),
            "route_to": route_to,
            "pct_routed": round(weighted["pct_routed"], 6),
        }
        subcatchment_section.append(subcatchment_entry)
        subarea_section.append(subarea_entry)
        records.append(
            {
                "subcatchment_id": sid,
                "method": "area_weighted_landuse_parameters",
                "components": components,
                "subcatchment": subcatchment_entry,
                "subarea": subarea_entry,
            }
        )

    return (
        {
            "ok": True,
            "mapping": "area_weighted_landuse_to_runoff_subarea",
            "sections": {"subcatchments": subcatchment_section, "subareas": subarea_section},
            "records": records,
            "unmatched_landuse_classes": sorted(unmatched),
        },
        audit,
        unmatched,
    )


def weighted_soil(
    rows: list[dict[str, Any]],
    lookup: dict[str, dict[str, Any]],
    default: dict[str, Any] | None,
    *,
    strict: bool,
) -> tuple[dict[str, Any], list[dict[str, Any]], set[str]]:
    by_id: dict[str, list[dict[str, Any]]] = {}
    for row in rows:
        by_id.setdefault(row["subcatchment_id"], []).append(row)

    records: list[dict[str, Any]] = []
    infiltration_section: list[dict[str, Any]] = []
    audit: list[dict[str, Any]] = []
    unmatched: set[str] = set()

    for sid in sorted(by_id):
        weighted = {field: 0.0 for field in NUMERIC_SOIL_FIELDS}
        components = []
        for row in by_id[sid]:
            texture = row["class"]
            fraction = float(row["fraction"])
            rec = lookup.get(normalize_key(texture))
            used_default = False
            if rec is None:
                if strict or default is None:
                    raise ValueError(f"Unmapped soil texture '{texture}' for subcatchment '{sid}'")
                rec = default
                used_default = True
                unmatched.add(texture)
            for field in NUMERIC_SOIL_FIELDS:
                weighted[field] += rec[field] * fraction
            components.append(
                {
                    "texture": texture,
                    "lookup_texture": rec["texture"],
                    "fraction": fraction,
                    "area_m2": row["area_m2"],
                    "used_default": used_default,
                }
            )
            audit.append(
                {
                    "subcatchment_id": sid,
                    "soil_texture": texture,
                    "lookup_texture": rec["texture"],
                    "area_m2": round(float(row["area_m2"]), 6),
                    "fraction": round(fraction, 8),
                    "used_default": used_default,
                }
            )
        entry = {
            "id": sid,
            "suction_mm": round(weighted["suction_mm"], 6),
            "ksat_mm_per_hr": round(weighted["ksat_mm_per_hr"], 6),
            "imdmax": round(weighted["imdmax"], 6),
        }
        infiltration_section.append(entry)
        records.append(
            {
                "subcatchment_id": sid,
                "method": "area_weighted_soil_green_ampt_parameters",
                "components": components,
                "infiltration": entry,
            }
        )

    return (
        {
            "ok": True,
            "mapping": "area_weighted_soil_to_green_ampt",
            "sections": {"infiltration": infiltration_section},
            "records": records,
            "unmatched_soil_textures": sorted(unmatched),
        },
        audit,
        unmatched,
    )


def merge_params(landuse: dict[str, Any], soil: dict[str, Any]) -> dict[str, Any]:
    subcatchments = {row["id"]: row for row in landuse["sections"]["subcatchments"]}
    subareas = {row["id"]: row for row in landuse["sections"]["subareas"]}
    infiltration = {row["id"]: row for row in soil["sections"]["infiltration"]}
    all_ids = sorted(set(subcatchments) | set(subareas) | set(infiltration))
    incomplete = []
    by_subcatchment = []
    for sid in all_ids:
        rec: dict[str, Any] = {"id": sid}
        missing = []
        for key, source, section in [
            ("subcatchment", subcatchments, "subcatchments"),
            ("subarea", subareas, "subareas"),
            ("infiltration", infiltration, "infiltration"),
        ]:
            if sid in source:
                rec[key] = source[sid]
            else:
                missing.append(section)
        if missing:
            rec["missing_sections"] = missing
            incomplete.append({"id": sid, "missing_sections": missing})
        by_subcatchment.append(rec)

    return {
        "ok": True,
        "mapping": "merged_area_weighted_swmm_params",
        "counts": {
            "subcatchment_count": len(all_ids),
            "subcatchments_with_subcatchment_section": len(subcatchments),
            "subcatchments_with_subarea_section": len(subareas),
            "subcatchments_with_infiltration_section": len(infiltration),
            "incomplete_subcatchment_count": len(incomplete),
        },
        "incomplete_ids": incomplete,
        "sections": {
            "subcatchments": [subcatchments[sid] for sid in sorted(subcatchments)],
            "subareas": [subareas[sid] for sid in sorted(subareas)],
            "infiltration": [infiltration[sid] for sid in sorted(infiltration)],
        },
        "by_subcatchment": by_subcatchment,
    }


def main() -> None:
    ap = argparse.ArgumentParser(description="Build area-weighted SWMM params from subcatchment/landuse/soil polygon overlays.")
    ap.add_argument("--subcatchments", type=Path, required=True)
    ap.add_argument("--landuse", type=Path, required=True)
    ap.add_argument("--soil", type=Path, required=True)
    ap.add_argument("--out-dir", type=Path, required=True)
    ap.add_argument("--id-field", default="basin_id")
    ap.add_argument("--landuse-field", default="CLASS")
    ap.add_argument("--soil-field", default="TEXTURE")
    ap.add_argument(
        "--skills-root",
        type=Path,
        default=None,
        help=(
            "Root directory containing sibling skills (e.g. swmm-params). "
            f"Overrides: flag > {SKILLS_ROOT_ENV} env var > default '<repo>/skills'."
        ),
    )
    ap.add_argument(
        "--landuse-lookup",
        type=Path,
        default=None,
        help="Lookup CSV for land use mapping (default: <skills-root>/swmm-params/references/landuse_class_to_subcatch_params.csv).",
    )
    ap.add_argument(
        "--soil-lookup",
        type=Path,
        default=None,
        help="Lookup CSV for soil texture mapping (default: <skills-root>/swmm-params/references/soil_texture_to_greenampt.csv).",
    )
    ap.add_argument("--strict", action="store_true", help="Fail on missing overlay coverage or missing lookup classes.")
    args = ap.parse_args()

    skills_root = resolve_skills_root(args.skills_root)
    params_dir = skills_root / "swmm-params"
    landuse_lookup_path = args.landuse_lookup or (params_dir / "references/landuse_class_to_subcatch_params.csv")
    soil_lookup_path = args.soil_lookup or (params_dir / "references/soil_texture_to_greenampt.csv")

    sub = read_vector(args.subcatchments, layer_name="subcatchments")
    landuse = read_vector(args.landuse, layer_name="landuse")
    soil = read_vector(args.soil, layer_name="soil")
    if args.id_field == "__feature_id__":
        sub = sub.reset_index(drop=True)
        sub[args.id_field] = [f"S{i + 1}" for i in range(len(sub))]
    elif args.id_field not in sub.columns:
        raise ValueError(f"Subcatchment id field '{args.id_field}' not found. Available fields: {list(sub.columns)}")
    ensure_projected(sub, layer_name="subcatchments")
    sub[args.id_field] = sub[args.id_field].astype(str)
    duplicate_ids = sorted(sub.loc[sub[args.id_field].duplicated(), args.id_field].unique())
    if duplicate_ids:
        sample = ", ".join(duplicate_ids[:10])
        raise ValueError(
            f"Subcatchment id field '{args.id_field}' is not unique; duplicate values include: {sample}. "
            "Use a unique id field or pass --id-field __feature_id__ to generate one id per feature."
        )

    land_lookup, land_default = load_landuse_lookup(landuse_lookup_path)
    soil_lookup, soil_default = load_soil_lookup(soil_lookup_path)
    if land_default is None and not args.strict:
        raise ValueError(f"Landuse lookup has no DEFAULT row: {landuse_lookup_path}")
    if soil_default is None and not args.strict:
        raise ValueError(f"Soil lookup has no '-' or DEFAULT row: {soil_lookup_path}")

    land_rows, land_issues = class_area_fractions(
        subcatchments=sub,
        thematic=landuse,
        id_field=args.id_field,
        class_field=args.landuse_field,
        default_class="DEFAULT",
        strict=args.strict,
        label="landuse",
    )
    soil_rows, soil_issues = class_area_fractions(
        subcatchments=sub,
        thematic=soil,
        id_field=args.id_field,
        class_field=args.soil_field,
        default_class="-",
        strict=args.strict,
        label="soil",
    )

    land_payload, land_audit, unmatched_land = weighted_landuse(land_rows, land_lookup, land_default, strict=args.strict)
    soil_payload, soil_audit, unmatched_soil = weighted_soil(soil_rows, soil_lookup, soil_default, strict=args.strict)
    merged = merge_params(land_payload, soil_payload)
    merged["sources"] = {
        "subcatchments": str(args.subcatchments),
        "landuse": str(args.landuse),
        "soil": str(args.soil),
        "landuse_lookup": str(landuse_lookup_path),
        "soil_lookup": str(soil_lookup_path),
    }
    merged["area_weighting"] = {
        "method": "polygon_intersection_area_fraction",
        "subcatchment_id_field": args.id_field,
        "landuse_field": args.landuse_field,
        "soil_field": args.soil_field,
        "missing_landuse_area_policy": "DEFAULT",
        "missing_soil_area_policy": "-",
        "soil_ksat_policy": "linear_area_weighted_first_draft",
    }
    merged["issues"] = land_issues + soil_issues
    merged["unmatched_landuse_classes"] = sorted(unmatched_land)
    merged["unmatched_soil_textures"] = sorted(unmatched_soil)

    # Build structured warnings: per-unmatched-class, sum of area routed
    # through the DEFAULT row. Lets framework_mcp_manifest promote these
    # automatically into missing_or_fallback_inputs.
    def _summarise_default_use(audit_rows, source_field, fallback_class):
        per_class: dict[str, float] = {}
        for row in audit_rows:
            if not row.get("used_default"):
                continue
            raw_class = row.get(source_field) or "<missing>"
            per_class[raw_class] = per_class.get(raw_class, 0.0) + float(row.get("area_m2") or 0.0)
        return [
            {
                "code": f"{source_field}_unmatched",
                "value": cls,
                "fallback_class": fallback_class,
                "fallback_area_m2": area,
            }
            for cls, area in sorted(per_class.items())
        ]

    warnings: list[dict[str, Any]] = []
    warnings.extend(_summarise_default_use(land_audit, "landuse_class", "DEFAULT"))
    warnings.extend(_summarise_default_use(soil_audit, "soil_texture", "DEFAULT"))
    merged["warnings"] = warnings

    out_dir = args.out_dir
    write_json(out_dir / "landuse_weighted_params.json", land_payload)
    write_json(out_dir / "soil_weighted_params.json", soil_payload)
    write_json(out_dir / "weighted_params.json", merged)
    write_csv(
        out_dir / "landuse_area_weights.csv",
        land_audit,
        ["subcatchment_id", "landuse_class", "lookup_landuse_class", "area_m2", "fraction", "used_default"],
    )
    write_csv(
        out_dir / "soil_area_weights.csv",
        soil_audit,
        ["subcatchment_id", "soil_texture", "lookup_texture", "area_m2", "fraction", "used_default"],
    )

    print(
        json.dumps(
            {
                "ok": True,
                "out_dir": str(out_dir),
                "weighted_params_json": str(out_dir / "weighted_params.json"),
                "landuse_area_weights_csv": str(out_dir / "landuse_area_weights.csv"),
                "soil_area_weights_csv": str(out_dir / "soil_area_weights.csv"),
                "subcatchment_count": merged["counts"]["subcatchment_count"],
                "issue_count": len(merged["issues"]),
                "unmatched_landuse_classes": merged["unmatched_landuse_classes"],
                "unmatched_soil_textures": merged["unmatched_soil_textures"],
                "warnings": warnings,
            },
            indent=2,
        )
    )


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Convert a municipal basin shapefile into SWMM-ready subcatchments.

The cold-start agent and the operator baseline both had to hand-pick one
or more polygons from a raw DrainageBasinBoundary shapefile, attach a
``subcatchment_id`` field, compute area, and synthesise width and slope.
This tool packages that into a single MCP-callable step.

Selection strategies (the ``--mode`` flag):

- ``by_id_field``: select features where ``--id-field`` equals
  ``--id-value`` (e.g. OBJECTID = 100). Single match expected.
- ``by_index``: select the feature at ``--index`` (0-based) in the
  layer's iteration order.
- ``largest``: select the single feature with the largest projected area.
- ``all``: every feature in the layer becomes its own subcatchment.

Width and slope are synthesised when the source layer does not carry
them. The defaults match the smoke runs:

- ``width_m = sqrt(area_m2)`` (geometric proxy)
- ``slope_pct = 1.0`` (flat-by-default placeholder; replace once a
  DEM-based slope tool exists; F12 in BACKLOG).
"""
from __future__ import annotations

import argparse
import csv
import json
import sys
from pathlib import Path

import geopandas as gpd


MODES = ("by_id_field", "by_index", "largest", "all")


def _ensure_projected(gdf: gpd.GeoDataFrame, layer_label: str) -> gpd.GeoDataFrame:
    if gdf.crs is None:
        raise ValueError(f"{layer_label} has no CRS; cannot compute area")
    if not gdf.crs.is_projected:
        raise ValueError(
            f"{layer_label} CRS {gdf.crs} is geographic; reproject to a projected CRS "
            "before computing subcatchment area / width."
        )
    return gdf


def _select(gdf: gpd.GeoDataFrame, mode: str, args: argparse.Namespace) -> gpd.GeoDataFrame:
    if mode == "by_id_field":
        if not args.id_field or args.id_value is None:
            raise ValueError("mode=by_id_field requires --id-field and --id-value")
        if args.id_field not in gdf.columns:
            raise ValueError(
                f"--id-field '{args.id_field}' not in basin layer columns: {list(gdf.columns)}"
            )
        # Match either as string or as the column's native dtype if it parses.
        col = gdf[args.id_field]
        match = gdf[col.astype(str) == str(args.id_value)]
        if len(match) == 0:
            raise ValueError(
                f"no basin matched {args.id_field}={args.id_value}; "
                f"sample values in field: {list(col.head(5).astype(str))}"
            )
        return match.reset_index(drop=True)
    if mode == "by_index":
        if args.index is None:
            raise ValueError("mode=by_index requires --index")
        if args.index < 0 or args.index >= len(gdf):
            raise ValueError(f"--index {args.index} out of range [0,{len(gdf) - 1}]")
        return gdf.iloc[[args.index]].reset_index(drop=True)
    if mode == "largest":
        gdf_sorted = gdf.copy()
        gdf_sorted["__area_m2__"] = gdf_sorted.geometry.area
        gdf_sorted = gdf_sorted.sort_values("__area_m2__", ascending=False)
        return gdf_sorted.iloc[[0]].drop(columns=["__area_m2__"]).reset_index(drop=True)
    if mode == "all":
        return gdf.reset_index(drop=True)
    raise ValueError(f"unknown mode: {mode}")


def _subcatchment_id(i: int, custom_prefix: str | None) -> str:
    prefix = custom_prefix if custom_prefix else "S"
    return f"{prefix}{i + 1}"


def parse_args() -> argparse.Namespace:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--basin-shp", required=True, help="Source basin shapefile or geojson path.")
    ap.add_argument("--mode", choices=MODES, default="by_id_field")
    ap.add_argument("--id-field", default="OBJECTID")
    ap.add_argument("--id-value", default=None)
    ap.add_argument("--index", type=int, default=None)
    ap.add_argument("--id-prefix", default="S", help="Prefix for generated subcatchment IDs.")
    ap.add_argument("--outlet-node-id", default="OUT1")
    ap.add_argument("--rain-gage-id", default="RG1")
    ap.add_argument("--default-slope-pct", type=float, default=1.0)
    ap.add_argument(
        "--width-method",
        choices=("sqrt_area",),
        default="sqrt_area",
        help="How to synthesise width when the source layer lacks one.",
    )
    ap.add_argument("--out-geojson", required=True)
    ap.add_argument("--out-csv", required=True)
    return ap.parse_args()


def main() -> None:
    args = parse_args()
    basin_path = Path(args.basin_shp)
    if not basin_path.exists():
        raise FileNotFoundError(basin_path)

    gdf = gpd.read_file(basin_path)
    if len(gdf) == 0:
        raise ValueError(f"basin layer is empty: {basin_path}")
    gdf = _ensure_projected(gdf, "basin")

    selected = _select(gdf, args.mode, args)
    if len(selected) == 0:
        raise ValueError("selection produced no features")

    rows: list[dict] = []
    out_features: list[dict] = []
    for i, (_, row) in enumerate(selected.iterrows()):
        sid = _subcatchment_id(i, args.id_prefix)
        geom = row.geometry
        area_m2 = float(geom.area)
        area_ha = area_m2 / 10000.0
        if args.width_method == "sqrt_area":
            width_m = float(area_m2) ** 0.5
        else:  # pragma: no cover  (only one supported value today)
            raise ValueError(f"unsupported width_method: {args.width_method}")
        rows.append({
            "subcatchment_id": sid,
            "outlet": args.outlet_node_id,
            "area_ha": area_ha,
            "width_m": width_m,
            "slope_pct": args.default_slope_pct,
            "rain_gage": args.rain_gage_id,
        })
        out_features.append({
            "type": "Feature",
            "properties": {"subcatchment_id": sid},
            "geometry": json.loads(gpd.GeoSeries([geom], crs=selected.crs).to_json())["features"][0]["geometry"],
        })

    out_geojson = Path(args.out_geojson)
    out_csv = Path(args.out_csv)
    out_geojson.parent.mkdir(parents=True, exist_ok=True)
    out_csv.parent.mkdir(parents=True, exist_ok=True)

    geojson_obj = {
        "type": "FeatureCollection",
        "name": "subcatchments",
        "crs": {"type": "name", "properties": {"name": f"urn:ogc:def:crs:EPSG::{selected.crs.to_epsg()}"}}
        if selected.crs.to_epsg() else None,
        "features": out_features,
    }
    # Drop crs key if it ended up as None (some CRS don't have EPSG codes).
    geojson_obj = {k: v for k, v in geojson_obj.items() if v is not None}
    out_geojson.write_text(json.dumps(geojson_obj, indent=2), encoding="utf-8")

    with out_csv.open("w", newline="", encoding="utf-8") as f:
        writer = csv.DictWriter(
            f, fieldnames=["subcatchment_id", "outlet", "area_ha", "width_m", "slope_pct", "rain_gage"]
        )
        writer.writeheader()
        writer.writerows(rows)

    report = {
        "ok": True,
        "skill": "swmm-gis",
        "tool": "basin_shp_to_subcatchments",
        "mode": args.mode,
        "counts": {
            "subcatchments_emitted": len(rows),
            "source_features_total": len(gdf),
        },
        "outlet_node_id": args.outlet_node_id,
        "rain_gage_id": args.rain_gage_id,
        "width_method": args.width_method,
        "default_slope_pct": args.default_slope_pct,
        "outputs": {
            "subcatchments_geojson": str(out_geojson),
            "subcatchments_csv": str(out_csv),
        },
        "inputs": {"basin_shp": str(basin_path)},
    }
    print(json.dumps(report, indent=2))


if __name__ == "__main__":
    try:
        main()
    except Exception as exc:
        print(f"basin_shp_to_subcatchments failed: {exc}", file=sys.stderr)
        raise

#!/usr/bin/env python3
"""Find DEM-based pour point candidates for watershed outlet selection.

Methods:
- boundary_min_elev: choose the minimum elevation cell on the DEM boundary.
- boundary_max_accum: compute D8 flow accumulation (with depression fill + flat resolution),
  then choose the boundary cell with maximum accumulation.

Outputs:
- GeoJSON point (same CRS as DEM)
- Preview PNG (DEM + outlet marker)

This is intended as a reproducible preprocessing step for SWMM experiments.
"""

from __future__ import annotations

import argparse
import json
from pathlib import Path

import numpy as np
import rasterio
from rasterio.transform import xy


def boundary_mask(shape: tuple[int, int]) -> np.ndarray:
    h, w = shape
    b = np.zeros((h, w), dtype=bool)
    b[0, :] = True
    b[-1, :] = True
    b[:, 0] = True
    b[:, -1] = True
    return b


def find_boundary_min_elev(dem: np.ma.MaskedArray) -> tuple[int, int, float]:
    b = boundary_mask(dem.shape)
    valid = b & (~dem.mask)
    vals = dem.data[valid]
    if vals.size == 0:
        raise RuntimeError("No valid border cells")
    minval = float(vals.min())
    idx = np.argwhere(valid & (dem.data == minval))[0]
    r, c = int(idx[0]), int(idx[1])
    return r, c, minval


def find_boundary_max_accum(dem_path: Path) -> tuple[int, int, float]:
    # Lazy import (heavier deps)
    from pysheds.grid import Grid

    grid = Grid.from_raster(str(dem_path))
    dem = grid.read_raster(str(dem_path))

    filled = grid.fill_depressions(dem)
    filled = grid.resolve_flats(filled)
    fdir = grid.flowdir(filled)
    acc = grid.accumulation(fdir)

    b = boundary_mask(acc.shape)
    valid = b & np.isfinite(acc)
    acc_border = np.where(valid, acc, -np.inf)
    r, c = np.unravel_index(np.argmax(acc_border), acc_border.shape)
    return int(r), int(c), float(acc_border[r, c])


def write_geojson_point(out_path: Path, x: float, y: float, props: dict):
    fc = {
        "type": "FeatureCollection",
        "features": [
            {
                "type": "Feature",
                "properties": props,
                "geometry": {"type": "Point", "coordinates": [x, y]},
            }
        ],
    }
    out_path.parent.mkdir(parents=True, exist_ok=True)
    out_path.write_text(json.dumps(fc, indent=2), encoding="utf-8")


def plot_preview_png(dem_path: Path, x: float, y: float, out_png: Path, title: str | None = None):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt

    with rasterio.open(dem_path) as ds:
        dem = ds.read(1, masked=True)
        extent = [ds.bounds.left, ds.bounds.right, ds.bounds.bottom, ds.bounds.top]
        crs = ds.crs

    vals = dem.compressed()
    vmin, vmax = np.quantile(vals, [0.02, 0.98])

    plt.rcParams.update({'font.family':'Arial','font.size':12})
    fig, ax = plt.subplots(figsize=(7, 5), dpi=200)
    im = ax.imshow(dem, extent=extent, origin='upper', cmap='terrain', vmin=vmin, vmax=vmax)
    ax.scatter([x], [y], c='red', s=50, marker='x', linewidths=2, label='Pour point')
    ax.set_xlabel('Easting (m)')
    ax.set_ylabel('Northing (m)')
    if title:
        ax.set_title(title)
    ax.legend(loc='lower left', framealpha=0.9)
    cb = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
    cb.set_label('Elevation')
    fig.tight_layout()
    out_png.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(out_png)


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument('--dem', required=True, type=Path)
    ap.add_argument('--method', required=True, choices=['boundary_min_elev', 'boundary_max_accum'])
    ap.add_argument('--out-geojson', required=True, type=Path)
    ap.add_argument('--out-png', required=True, type=Path)
    ap.add_argument('--name', default='pour_point')
    args = ap.parse_args()

    with rasterio.open(args.dem) as ds:
        dem = ds.read(1, masked=True)
        transform = ds.transform
        crs = ds.crs

    if args.method == 'boundary_min_elev':
        r, c, score = find_boundary_min_elev(dem)
        props = {"name": args.name, "method": args.method, "score": score}
    else:
        r, c, score = find_boundary_max_accum(args.dem)
        props = {"name": args.name, "method": args.method, "score": score}

    x, y = xy(transform, r, c, offset='center')
    props.update({"row": r, "col": c, "crs": crs.to_string() if crs else None})

    write_geojson_point(args.out_geojson, float(x), float(y), props)
    plot_preview_png(args.dem, float(x), float(y), args.out_png, title=None)

    print(json.dumps({"x": float(x), "y": float(y), **props}, indent=2))


if __name__ == '__main__':
    main()

#!/usr/bin/env python3
from __future__ import annotations

import argparse
from pathlib import Path

import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import rasterio
from matplotlib.lines import Line2D
from pyproj import Transformer


def raster_extent(dataset: rasterio.DatasetReader) -> tuple[float, float, float, float]:
    bounds = dataset.bounds
    return bounds.left, bounds.right, bounds.bottom, bounds.top


def add_lonlat_ticks(ax: plt.Axes, crs: str, extent: tuple[float, float, float, float]) -> None:
    left, right, bottom, top = extent
    xs = np.linspace(left, right, 5)
    ys = np.linspace(bottom, top, 5)
    transformer = Transformer.from_crs(crs, "EPSG:4326", always_xy=True)
    lon_labels = [transformer.transform(float(x), bottom)[0] for x in xs]
    lat_labels = [transformer.transform(left, float(y))[1] for y in ys]
    ax.set_xticks(xs)
    ax.set_yticks(ys)
    ax.set_xticklabels([f"{lon:.3f}°".replace("-", "−") for lon in lon_labels], fontsize=11)
    ax.set_yticklabels([f"{lat:.3f}°".replace("-", "−") for lat in lat_labels], fontsize=11)
    ax.tick_params(top=True, right=True, labeltop=False, labelright=False, direction="in", length=5, width=1.0)
    ax.set_xlabel("Longitude", fontsize=12)
    ax.set_ylabel("Latitude", fontsize=12)


def plot_overview(
    *,
    slope: Path,
    subcatchments: Path,
    flow: Path,
    outfall: Path,
    out_png: Path,
    title: str,
) -> None:
    with rasterio.open(slope) as ds:
        arr = ds.read(1, masked=True).astype(float)
        extent = raster_extent(ds)
        crs = ds.crs.to_string()

    data = arr.compressed()
    vmax = float(np.percentile(data, 98)) if data.size else 1.0
    vmin = float(np.percentile(data, 2)) if data.size else 0.0

    sub = gpd.read_file(subcatchments)
    flow_gdf = gpd.read_file(flow)
    outfall_gdf = gpd.read_file(outfall)
    if sub.crs and sub.crs.to_string() != crs:
        sub = sub.to_crs(crs)
    if flow_gdf.crs and flow_gdf.crs.to_string() != crs:
        flow_gdf = flow_gdf.to_crs(crs)
    if outfall_gdf.crs and outfall_gdf.crs.to_string() != crs:
        outfall_gdf = outfall_gdf.to_crs(crs)

    plt.rcParams.update(
        {
            "font.family": ["Arial", "Helvetica", "DejaVu Sans"],
            "axes.linewidth": 0.8,
            "figure.dpi": 180,
            "savefig.dpi": 240,
        }
    )
    fig, ax = plt.subplots(figsize=(7.2, 9.0), constrained_layout=True)
    image = ax.imshow(
        arr,
        extent=extent,
        origin="upper",
        cmap="RdYlGn_r",
        vmin=vmin,
        vmax=vmax,
        interpolation="nearest",
        alpha=0.62,
    )
    sub.boundary.plot(ax=ax, color="#343a40", linewidth=1.15, alpha=0.96)
    if not flow_gdf.empty:
        flow_gdf.plot(ax=ax, color="#0072B2", linewidth=1.85, alpha=0.96)
    if not outfall_gdf.empty:
        outfall_gdf.plot(ax=ax, marker="*", color="#D55E00", edgecolor="white", linewidth=0.6, markersize=130)

    if title:
        ax.set_title(title, fontsize=14, pad=8)
    ax.set_xlim(extent[0], extent[1])
    ax.set_ylim(extent[2], extent[3])
    ax.set_aspect("equal")
    add_lonlat_ticks(ax, crs, extent)
    ax.grid(color="#b8bec6", linewidth=0.35, alpha=0.6)

    cbar = fig.colorbar(image, ax=ax, fraction=0.034, pad=0.018)
    cbar.set_label("Slope (%)", fontsize=12)
    cbar.ax.tick_params(labelsize=11)

    handles = [
        Line2D([0], [0], color="#343a40", linewidth=2.1, label="Subcatchment boundary"),
        Line2D([0], [0], color="#0072B2", linewidth=3.0, label="Flow path"),
        Line2D([0], [0], marker="*", color="none", markerfacecolor="#D55E00", markeredgecolor="white", markersize=13, label="Outfall"),
    ]
    ax.legend(handles=handles, loc="lower right", frameon=True, framealpha=0.94, fontsize=11)
    out_png.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(out_png, bbox_inches="tight")
    plt.close(fig)


def main() -> None:
    parser = argparse.ArgumentParser(description="Plot clean QGIS/GRASS standard watershed final layers.")
    parser.add_argument("--slope", type=Path, required=True)
    parser.add_argument("--subcatchments", type=Path, required=True)
    parser.add_argument("--flow", type=Path, required=True)
    parser.add_argument("--outfall", type=Path, required=True)
    parser.add_argument("--out-png", type=Path, required=True)
    parser.add_argument("--title", default="")
    args = parser.parse_args()
    plot_overview(
        slope=args.slope,
        subcatchments=args.subcatchments,
        flow=args.flow,
        outfall=args.outfall,
        out_png=args.out_png,
        title=args.title,
    )


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
from __future__ import annotations

import argparse
import csv
import json
import math
from pathlib import Path
from typing import Any


FloatCandidate = tuple[str, dict[str, Any], str]


def load_json(path: Path) -> Any:
    return json.loads(path.read_text(encoding="utf-8"))


def write_json(path: Path, obj: Any) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(json.dumps(obj, indent=2), encoding="utf-8")


def write_csv(path: Path, rows: list[dict[str, Any]], headers: list[str]) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", encoding="utf-8", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=headers)
        writer.writeheader()
        writer.writerows(rows)


def parse_xy(value: Any) -> tuple[float, float]:
    if not isinstance(value, list) or len(value) < 2:
        raise ValueError(f"Invalid coordinate value: {value}")
    return float(value[0]), float(value[1])


def ensure_closed_ring(ring: list[list[float]]) -> list[list[float]]:
    if len(ring) < 3:
        raise ValueError("Polygon ring must have at least 3 points")
    first = ring[0]
    last = ring[-1]
    if float(first[0]) == float(last[0]) and float(first[1]) == float(last[1]):
        return ring
    return ring + [first]


def ring_metrics(ring: list[list[float]]) -> tuple[float, float, float, float]:
    closed = ensure_closed_ring(ring)

    double_area = 0.0
    cx_num = 0.0
    cy_num = 0.0
    perimeter = 0.0

    for i in range(len(closed) - 1):
        x0, y0 = parse_xy(closed[i])
        x1, y1 = parse_xy(closed[i + 1])
        cross = x0 * y1 - x1 * y0
        double_area += cross
        cx_num += (x0 + x1) * cross
        cy_num += (y0 + y1) * cross
        perimeter += math.hypot(x1 - x0, y1 - y0)

    area = 0.5 * double_area
    if abs(area) < 1e-12:
        xs = [parse_xy(pt)[0] for pt in closed[:-1]]
        ys = [parse_xy(pt)[1] for pt in closed[:-1]]
        return 0.0, sum(xs) / len(xs), sum(ys) / len(ys), perimeter

    cx = cx_num / (6.0 * area)
    cy = cy_num / (6.0 * area)
    return area, cx, cy, perimeter


def polygon_metrics(polygon_coords: list[Any]) -> tuple[float, float, float, float]:
    if not polygon_coords:
        raise ValueError("Polygon has no rings")

    outer = polygon_coords[0]
    a_outer, cx_outer, cy_outer, outer_perimeter = ring_metrics(outer)
    area_outer = abs(a_outer)
    if area_outer < 1e-12:
        raise ValueError("Polygon outer ring has zero area")

    hole_area_sum = 0.0
    hole_cx_sum = 0.0
    hole_cy_sum = 0.0
    for hole in polygon_coords[1:]:
        a_hole, cx_hole, cy_hole, _ = ring_metrics(hole)
        ah = abs(a_hole)
        hole_area_sum += ah
        hole_cx_sum += cx_hole * ah
        hole_cy_sum += cy_hole * ah

    net_area = area_outer - hole_area_sum
    if net_area <= 1e-12:
        raise ValueError("Polygon net area is <= 0 after holes")

    cx = (cx_outer * area_outer - hole_cx_sum) / net_area
    cy = (cy_outer * area_outer - hole_cy_sum) / net_area
    return net_area, cx, cy, outer_perimeter


def geometry_metrics(geometry: dict[str, Any]) -> tuple[float, float, float, float]:
    gtype = geometry.get("type")
    coords = geometry.get("coordinates")
    if gtype == "Polygon":
        return polygon_metrics(coords)

    if gtype == "MultiPolygon":
        total_area = 0.0
        total_perimeter = 0.0
        cx_weighted = 0.0
        cy_weighted = 0.0
        for poly in coords:
            area, cx, cy, perimeter = polygon_metrics(poly)
            total_area += area
            total_perimeter += perimeter
            cx_weighted += cx * area
            cy_weighted += cy * area
        if total_area <= 1e-12:
            raise ValueError("MultiPolygon total area is <= 0")
        return total_area, cx_weighted / total_area, cy_weighted / total_area, total_perimeter

    raise ValueError(f"Unsupported geometry type: {gtype}")


def parse_optional_float_value(raw: Any, *, context: str) -> float | None:
    if raw is None:
        return None
    if isinstance(raw, bool):
        raise ValueError(f"Invalid float in '{context}': {raw}")
    if isinstance(raw, (int, float)):
        return float(raw)
    text = str(raw).strip()
    if not text:
        return None
    try:
        return float(text)
    except ValueError as exc:
        raise ValueError(f"Invalid float in '{context}': {raw}") from exc


def parse_optional_float(mapping: dict[str, Any], field: str, *, mapping_name: str = "properties") -> float | None:
    if field not in mapping:
        return None
    return parse_optional_float_value(mapping.get(field), context=f"{mapping_name}.{field}")


def pick_first_float(candidates: list[FloatCandidate]) -> tuple[float, str] | None:
    for mapping_name, mapping, field in candidates:
        value = parse_optional_float(mapping, field, mapping_name=mapping_name)
        if value is None:
            continue
        return value, f"{mapping_name}:{field}"
    return None


def normalize_non_blank(value: Any) -> str | None:
    if value is None:
        return None
    token = str(value).strip()
    return token or None


def increment_count(counter: dict[str, int], key: str) -> None:
    counter[key] = counter.get(key, 0) + 1


def first_present_id(record: dict[str, Any], fields: list[str], *, context: str) -> str:
    for field in fields:
        value = normalize_non_blank(record.get(field))
        if value is not None:
            return value
    raise ValueError(f"{context} missing id field; tried {fields}")


def build_dem_stats_index(
    dem_stats_json: Path | None,
    *,
    id_field: str,
    dem_stats_id_field: str,
) -> tuple[dict[str, dict[str, Any]], str | None]:
    if dem_stats_json is None:
        return {}, None

    raw = load_json(dem_stats_json)
    candidates = [dem_stats_id_field, id_field, "subcatchment_id", "id"]

    records: list[tuple[str, dict[str, Any]]] = []

    def append_record(record_id: str, record: dict[str, Any], *, context: str) -> None:
        if not isinstance(record, dict):
            raise ValueError(f"{context} must be an object")
        rid = normalize_non_blank(record_id)
        if rid is None:
            rid = first_present_id(record, candidates, context=context)
        records.append((rid, record))

    if isinstance(raw, dict) and "subcatchments" in raw:
        body = raw.get("subcatchments")
        if isinstance(body, list):
            for idx, entry in enumerate(body, start=1):
                if not isinstance(entry, dict):
                    raise ValueError(
                        f"DEM stats entry {idx} in {dem_stats_json} must be an object"
                    )
                rid = first_present_id(entry, candidates, context=f"DEM stats entry {idx}")
                records.append((rid, entry))
        elif isinstance(body, dict):
            for key, entry in body.items():
                append_record(str(key), entry, context=f"DEM stats entry '{key}'")
        else:
            raise ValueError(
                f"DEM stats file {dem_stats_json} field 'subcatchments' must be a list or object"
            )
    elif isinstance(raw, list):
        for idx, entry in enumerate(raw, start=1):
            if not isinstance(entry, dict):
                raise ValueError(f"DEM stats entry {idx} in {dem_stats_json} must be an object")
            rid = first_present_id(entry, candidates, context=f"DEM stats entry {idx}")
            records.append((rid, entry))
    elif isinstance(raw, dict) and all(isinstance(v, dict) for v in raw.values()):
        for key, entry in raw.items():
            append_record(str(key), entry, context=f"DEM stats entry '{key}'")
    else:
        raise ValueError(
            "DEM stats JSON must be one of: list[object], "
            "{'subcatchments': list|object}, or object keyed by subcatchment id"
        )

    out: dict[str, dict[str, Any]] = {}
    for rid, entry in records:
        if rid in out:
            raise ValueError(f"Duplicate DEM stats id '{rid}' in {dem_stats_json}")
        out[rid] = dict(entry)

    return out, str(dem_stats_json)


def build_node_index(network_json: Path) -> dict[str, tuple[float, float]]:
    network = load_json(network_json)
    nodes: dict[str, tuple[float, float]] = {}

    for node in list(network.get("junctions") or []) + list(network.get("outfalls") or []):
        node_id = str(node.get("id") or "").strip()
        if not node_id:
            raise ValueError(f"Node with missing id in {network_json}")
        coords = node.get("coordinates")
        if not isinstance(coords, dict):
            raise ValueError(f"Node '{node_id}' missing coordinates in {network_json}")
        x = float(coords.get("x"))
        y = float(coords.get("y"))
        nodes[node_id] = (x, y)

    if not nodes:
        raise ValueError(f"No nodes found in network JSON: {network_json}")
    return nodes


def nearest_node(
    centroid_x: float,
    centroid_y: float,
    nodes: dict[str, tuple[float, float]],
) -> tuple[str, float]:
    best_id: str | None = None
    best_dist = float("inf")
    for node_id, (nx, ny) in nodes.items():
        d = math.hypot(centroid_x - nx, centroid_y - ny)
        if d < best_dist:
            best_dist = d
            best_id = node_id
    if best_id is None:
        raise ValueError("No nodes available for nearest-node search")
    return best_id, best_dist


def clamp_with_min(
    value: float,
    *,
    min_value: float,
    source: str,
    metric: str,
    diagnostics: list[str],
) -> tuple[float, str]:
    if value >= min_value:
        return value, source
    diagnostics.append(
        f"{metric} from {source}={value:.6f} below min {min_value:.6f}; clamped."
    )
    return min_value, f"{source}|clamped_min"


def resolve_dem_flow_length(
    props: dict[str, Any],
    dem_stats: dict[str, Any],
    diagnostics: list[str],
) -> tuple[float, str] | None:
    dem_flow = pick_first_float(
        [
            ("dem_stats", dem_stats, "dem_flow_length_m"),
            ("dem_stats", dem_stats, "flow_length_m"),
            ("properties", props, "dem_flow_length_m"),
            ("properties", props, "flow_length_m"),
        ]
    )
    if dem_flow is None:
        return None

    value, source = dem_flow
    if value <= 0:
        diagnostics.append(f"Ignoring non-positive flow length from {source}={value:.6f}.")
        return None
    return value, source


def estimate_width_m(
    area_m2: float,
    perimeter_m: float,
    *,
    min_width_m: float,
    props: dict[str, Any],
    dem_stats: dict[str, Any],
    dem_flow_length: tuple[float, str] | None,
    diagnostics: list[str],
) -> tuple[float, str]:
    width_direct = pick_first_float(
        [
            ("properties", props, "width_m"),
            ("properties", props, "hydraulic_width_m"),
            ("dem_stats", dem_stats, "dem_width_m"),
            ("dem_stats", dem_stats, "width_m"),
        ]
    )
    if width_direct is not None:
        width, source = width_direct
        if width > 0:
            return clamp_with_min(
                width,
                min_value=min_width_m,
                source=source,
                metric="width_m",
                diagnostics=diagnostics,
            )
        diagnostics.append(f"Ignoring non-positive width from {source}={width:.6f}.")

    if dem_flow_length is not None:
        flow_length_m, flow_source = dem_flow_length
        width = area_m2 / max(flow_length_m, 1e-9)
        width, width_source = clamp_with_min(
            width,
            min_value=min_width_m,
            source=f"derived:area_m2/{flow_source}",
            metric="width_m",
            diagnostics=diagnostics,
        )
        return width, width_source

    # Deterministic surrogate: equivalent hydraulic width from area and perimeter.
    width_geom = 2.0 * area_m2 / max(perimeter_m, 1e-9)
    width_geom, width_source = clamp_with_min(
        width_geom,
        min_value=min_width_m,
        source="derived:2*area_m2/perimeter_m",
        metric="width_m",
        diagnostics=diagnostics,
    )
    return width_geom, width_source


def estimate_flow_length_m(
    area_m2: float,
    width_m: float,
    *,
    width_source: str,
    dem_flow_length: tuple[float, str] | None,
) -> tuple[float, str]:
    if dem_flow_length is not None:
        value, source = dem_flow_length
        return value, source
    return area_m2 / max(width_m, 1e-9), f"derived:area_m2/width_m ({width_source})"


def estimate_slope_pct(
    props: dict[str, Any],
    dem_stats: dict[str, Any],
    flow_length_m: float,
    *,
    flow_length_source: str,
    default_slope_pct: float,
    min_slope_pct: float,
    diagnostics: list[str],
) -> tuple[float, str]:
    slope_direct = pick_first_float([("properties", props, "slope_pct")])
    if slope_direct is not None:
        slope, source = slope_direct
        return clamp_with_min(
            slope,
            min_value=min_slope_pct,
            source=source,
            metric="slope_pct",
            diagnostics=diagnostics,
        )

    dem_slope_direct = pick_first_float(
        [
            ("dem_stats", dem_stats, "dem_slope_pct"),
            ("dem_stats", dem_stats, "raster_slope_pct"),
            ("dem_stats", dem_stats, "mean_slope_pct"),
            ("dem_stats", dem_stats, "slope_pct"),
            ("properties", props, "dem_slope_pct"),
            ("properties", props, "raster_slope_pct"),
        ]
    )
    if dem_slope_direct is not None:
        slope, source = dem_slope_direct
        return clamp_with_min(
            slope,
            min_value=min_slope_pct,
            source=source,
            metric="slope_pct",
            diagnostics=diagnostics,
        )

    dem_mean = pick_first_float(
        [
            ("dem_stats", dem_stats, "dem_elev_mean_m"),
            ("dem_stats", dem_stats, "elev_mean_m"),
            ("properties", props, "dem_elev_mean_m"),
        ]
    )
    dem_outlet = pick_first_float(
        [
            ("dem_stats", dem_stats, "dem_elev_outlet_m"),
            ("dem_stats", dem_stats, "elev_outlet_m"),
            ("properties", props, "dem_elev_outlet_m"),
        ]
    )
    if dem_mean is not None and dem_outlet is not None:
        mean_value, mean_source = dem_mean
        outlet_value, outlet_source = dem_outlet
        slope = (mean_value - outlet_value) / max(flow_length_m, 1e-9) * 100.0
        return clamp_with_min(
            slope,
            min_value=min_slope_pct,
            source=f"derived:({mean_source}-{outlet_source})/{flow_length_source}*100",
            metric="slope_pct",
            diagnostics=diagnostics,
        )

    dem_max = pick_first_float(
        [
            ("dem_stats", dem_stats, "dem_elev_max_m"),
            ("dem_stats", dem_stats, "elev_max_m"),
            ("properties", props, "dem_elev_max_m"),
        ]
    )
    dem_min = pick_first_float(
        [
            ("dem_stats", dem_stats, "dem_elev_min_m"),
            ("dem_stats", dem_stats, "elev_min_m"),
            ("properties", props, "dem_elev_min_m"),
        ]
    )
    if dem_max is not None and dem_min is not None:
        max_value, max_source = dem_max
        min_value, min_source = dem_min
        slope = (max_value - min_value) / max(flow_length_m, 1e-9) * 100.0
        return clamp_with_min(
            slope,
            min_value=min_slope_pct,
            source=f"derived:({max_source}-{min_source})/{flow_length_source}*100",
            metric="slope_pct",
            diagnostics=diagnostics,
        )

    if dem_mean is not None and dem_min is not None:
        mean_value, mean_source = dem_mean
        min_value, min_source = dem_min
        slope = (mean_value - min_value) / max(flow_length_m, 1e-9) * 100.0
        return clamp_with_min(
            slope,
            min_value=min_slope_pct,
            source=f"derived:({mean_source}-{min_source})/{flow_length_source}*100",
            metric="slope_pct",
            diagnostics=diagnostics,
        )

    elev_mean = pick_first_float([("properties", props, "elev_mean_m")])
    elev_outlet = pick_first_float([("properties", props, "elev_outlet_m")])
    if elev_mean is not None and elev_outlet is not None:
        mean_value, mean_source = elev_mean
        outlet_value, outlet_source = elev_outlet
        slope = (mean_value - outlet_value) / max(flow_length_m, 1e-9) * 100.0
        return clamp_with_min(
            slope,
            min_value=min_slope_pct,
            source=f"derived:({mean_source}-{outlet_source})/{flow_length_source}*100",
            metric="slope_pct",
            diagnostics=diagnostics,
        )

    return clamp_with_min(
        default_slope_pct,
        min_value=min_slope_pct,
        source="default_slope_pct",
        metric="slope_pct",
        diagnostics=diagnostics,
    )


def link_outlet(
    *,
    subcatchment_id: str,
    props: dict[str, Any],
    outlet_hint_field: str,
    centroid_x: float,
    centroid_y: float,
    node_index: dict[str, tuple[float, float]],
) -> tuple[str, float, str, str, list[str]]:
    diagnostics: list[str] = []
    hint = str(props.get(outlet_hint_field) or "").strip()

    if hint:
        if hint in node_index:
            nx, ny = node_index[hint]
            outlet_distance = math.hypot(centroid_x - nx, centroid_y - ny)
            return hint, outlet_distance, f"hint:{outlet_hint_field}", hint, diagnostics

        diagnostics.append(
            f"Feature '{subcatchment_id}' has unknown outlet hint '{hint}' in "
            f"properties.{outlet_hint_field}; used nearest node fallback."
        )
        outlet, outlet_distance = nearest_node(centroid_x, centroid_y, node_index)
        return (
            outlet,
            outlet_distance,
            f"nearest_node_fallback:invalid_hint:{outlet_hint_field}",
            hint,
            diagnostics,
        )

    diagnostics.append(
        f"Feature '{subcatchment_id}' has blank properties.{outlet_hint_field}; used nearest node fallback."
    )
    outlet, outlet_distance = nearest_node(centroid_x, centroid_y, node_index)
    return outlet, outlet_distance, "nearest_node", "", diagnostics


def is_dem_assisted_source(source: str) -> bool:
    token = source.lower()
    return "dem_stats" in token or "dem_" in token or "raster" in token


def main() -> None:
    ap = argparse.ArgumentParser(
        description=(
            "Preprocess subcatchment polygons into builder-ready CSV with deterministic "
            "width/slope/outlet linking and optional DEM-assisted metrics."
        )
    )
    ap.add_argument("--subcatchments-geojson", type=Path, required=True)
    ap.add_argument("--network-json", type=Path, required=True)
    ap.add_argument("--out-csv", type=Path, required=True)
    ap.add_argument("--out-json", type=Path, required=True)
    ap.add_argument("--id-field", default="subcatchment_id")
    ap.add_argument("--outlet-hint-field", default="outlet_hint")
    ap.add_argument("--dem-stats-json", type=Path, default=None)
    ap.add_argument("--dem-stats-id-field", default="subcatchment_id")
    ap.add_argument("--default-slope-pct", type=float, default=1.0)
    ap.add_argument("--min-slope-pct", type=float, default=0.1)
    ap.add_argument("--min-width-m", type=float, default=10.0)
    ap.add_argument("--default-curb-length-m", type=float, default=0.0)
    ap.add_argument("--default-rain-gage", default="")
    ap.add_argument("--max-link-distance-m", type=float, default=None)
    args = ap.parse_args()

    if args.default_slope_pct <= 0:
        raise ValueError("--default-slope-pct must be > 0")
    if args.min_slope_pct <= 0:
        raise ValueError("--min-slope-pct must be > 0")
    if args.min_width_m <= 0:
        raise ValueError("--min-width-m must be > 0")

    gj = load_json(args.subcatchments_geojson)
    if gj.get("type") != "FeatureCollection":
        raise ValueError("subcatchments GeoJSON must be a FeatureCollection")
    features = gj.get("features") or []
    if not isinstance(features, list) or not features:
        raise ValueError("subcatchments GeoJSON has no features")

    node_index = build_node_index(args.network_json)
    dem_stats_index, dem_stats_source = build_dem_stats_index(
        args.dem_stats_json,
        id_field=args.id_field,
        dem_stats_id_field=args.dem_stats_id_field,
    )

    seen_ids: set[str] = set()
    used_dem_ids: set[str] = set()
    diagnostics: list[dict[str, str]] = []
    csv_rows: list[dict[str, Any]] = []
    detail_rows: list[dict[str, Any]] = []

    width_source_counts: dict[str, int] = {}
    slope_source_counts: dict[str, int] = {}
    outlet_method_counts: dict[str, int] = {}

    for idx, feature in enumerate(features):
        props = feature.get("properties") or {}
        geom = feature.get("geometry")
        if not isinstance(geom, dict):
            raise ValueError(f"Feature index {idx} has no geometry")

        feature_id = props.get(args.id_field)
        if feature_id is None:
            feature_id = props.get("id")
        if feature_id is None:
            feature_id = feature.get("id")
        subcatchment_id = str(feature_id or "").strip()
        if not subcatchment_id:
            raise ValueError(f"Feature index {idx} missing id field '{args.id_field}'")
        if subcatchment_id in seen_ids:
            raise ValueError(f"Duplicate subcatchment id '{subcatchment_id}'")
        seen_ids.add(subcatchment_id)

        dem_stats = dem_stats_index.get(subcatchment_id, {})
        if dem_stats:
            used_dem_ids.add(subcatchment_id)

        row_diagnostics: list[str] = []
        if args.dem_stats_json is not None and not dem_stats:
            row_diagnostics.append(
                f"No DEM stats record matched id '{subcatchment_id}'; using deterministic fallback where needed."
            )

        area_m2, cx, cy, perimeter_m = geometry_metrics(geom)
        area_ha = area_m2 / 10000.0

        dem_flow_length = resolve_dem_flow_length(props, dem_stats, row_diagnostics)
        width_m, width_source = estimate_width_m(
            area_m2,
            perimeter_m,
            min_width_m=args.min_width_m,
            props=props,
            dem_stats=dem_stats,
            dem_flow_length=dem_flow_length,
            diagnostics=row_diagnostics,
        )
        flow_length_m, flow_length_source = estimate_flow_length_m(
            area_m2,
            width_m,
            width_source=width_source,
            dem_flow_length=dem_flow_length,
        )
        slope_pct, slope_source = estimate_slope_pct(
            props,
            dem_stats,
            flow_length_m,
            flow_length_source=flow_length_source,
            default_slope_pct=args.default_slope_pct,
            min_slope_pct=args.min_slope_pct,
            diagnostics=row_diagnostics,
        )

        outlet, outlet_distance, outlet_method, outlet_hint, outlet_diag = link_outlet(
            subcatchment_id=subcatchment_id,
            props=props,
            outlet_hint_field=args.outlet_hint_field,
            centroid_x=cx,
            centroid_y=cy,
            node_index=node_index,
        )
        row_diagnostics.extend(outlet_diag)

        if args.max_link_distance_m is not None and outlet_distance > args.max_link_distance_m:
            raise ValueError(
                f"Feature '{subcatchment_id}' linked outlet '{outlet}' via {outlet_method} at distance "
                f"{outlet_distance:.3f} m exceeding --max-link-distance-m={args.max_link_distance_m}"
            )

        curb_length_m = parse_optional_float(props, "curb_length_m")
        if curb_length_m is None:
            curb_length_m = args.default_curb_length_m

        snow_pack = str(props.get("snow_pack") or "").strip()
        rain_gage = str(props.get("rain_gage") or args.default_rain_gage).strip()

        dem_assisted = bool(dem_stats) and (
            is_dem_assisted_source(width_source)
            or is_dem_assisted_source(flow_length_source)
            or is_dem_assisted_source(slope_source)
        )

        increment_count(width_source_counts, width_source)
        increment_count(slope_source_counts, slope_source)
        increment_count(outlet_method_counts, outlet_method)

        for message in row_diagnostics:
            diagnostics.append({"id": subcatchment_id, "message": message})

        csv_row = {
            "subcatchment_id": subcatchment_id,
            "outlet": outlet,
            "area_ha": round(area_ha, 6),
            "width_m": round(width_m, 6),
            "slope_pct": round(slope_pct, 6),
            "curb_length_m": round(curb_length_m, 6),
            "snow_pack": snow_pack,
            "rain_gage": rain_gage,
            "area_source": "geometry:planar_polygon",
            "width_source": width_source,
            "flow_length_source": flow_length_source,
            "slope_source": slope_source,
            "outlet_method": outlet_method,
            "outlet_distance_m": round(outlet_distance, 6),
            "outlet_hint": outlet_hint,
            "outlet_diagnostic": " | ".join(outlet_diag),
            "dem_assisted": "yes" if dem_assisted else "no",
        }
        csv_rows.append(csv_row)

        detail_rows.append(
            {
                "id": subcatchment_id,
                "area_m2": area_m2,
                "perimeter_m": perimeter_m,
                "centroid": {"x": cx, "y": cy},
                "flow_length_m": flow_length_m,
                "flow_length_source": flow_length_source,
                "width_m": width_m,
                "width_source": width_source,
                "slope_pct": slope_pct,
                "slope_source": slope_source,
                "outlet": outlet,
                "outlet_distance_m": outlet_distance,
                "outlet_method": outlet_method,
                "outlet_hint": outlet_hint,
                "dem_stats_used": bool(dem_stats),
                "dem_assisted": dem_assisted,
                "diagnostics": row_diagnostics,
                "sources": {
                    "area": "geometry:planar_polygon",
                    "width": width_source,
                    "flow_length": flow_length_source,
                    "slope": slope_source,
                    "outlet": outlet_method,
                },
            }
        )

    csv_rows.sort(key=lambda r: str(r["subcatchment_id"]))
    detail_rows.sort(key=lambda r: str(r["id"]))

    csv_headers = [
        "subcatchment_id",
        "outlet",
        "area_ha",
        "width_m",
        "slope_pct",
        "curb_length_m",
        "snow_pack",
        "rain_gage",
        "area_source",
        "width_source",
        "flow_length_source",
        "slope_source",
        "outlet_method",
        "outlet_distance_m",
        "outlet_hint",
        "outlet_diagnostic",
        "dem_assisted",
    ]
    write_csv(args.out_csv, csv_rows, csv_headers)

    unmatched_dem_ids = sorted(set(dem_stats_index) - used_dem_ids)
    for dem_id in unmatched_dem_ids:
        diagnostics.append(
            {
                "id": dem_id,
                "message": "DEM stats record did not match any subcatchment feature id.",
            }
        )

    report = {
        "ok": True,
        "skill": "swmm-gis",
        "inputs": {
            "subcatchments_geojson": str(args.subcatchments_geojson),
            "network_json": str(args.network_json),
            "dem_stats_json": dem_stats_source,
        },
        "assumptions": {
            "planar_coordinates": True,
            "coordinate_units": "meters",
            "width_priority": [
                "properties.width_m / properties.hydraulic_width_m",
                "DEM flow length -> area_m2 / flow_length_m",
                "2 * area_m2 / perimeter_m",
            ],
            "flow_length_priority": [
                "DEM flow length fields",
                "area_m2 / width_m",
            ],
            "slope_priority": [
                "properties.slope_pct",
                "DEM direct slope fields",
                "(DEM elevation stats) / flow_length_m",
                "(properties.elev_mean_m - properties.elev_outlet_m) / flow_length_m * 100",
                "default_slope_pct",
            ],
            "outlet_link_priority": [
                f"properties.{args.outlet_hint_field} (if valid)",
                "nearest network node fallback",
            ],
        },
        "parameters": {
            "id_field": args.id_field,
            "outlet_hint_field": args.outlet_hint_field,
            "dem_stats_id_field": args.dem_stats_id_field,
            "default_slope_pct": args.default_slope_pct,
            "min_slope_pct": args.min_slope_pct,
            "min_width_m": args.min_width_m,
            "default_curb_length_m": args.default_curb_length_m,
            "default_rain_gage": args.default_rain_gage,
            "max_link_distance_m": args.max_link_distance_m,
        },
        "counts": {
            "feature_count": len(features),
            "subcatchment_count": len(detail_rows),
            "network_node_count": len(node_index),
            "dem_stats_count": len(dem_stats_index),
            "dem_stats_matched_count": len(used_dem_ids),
            "diagnostic_count": len(diagnostics),
        },
        "method_counts": {
            "width_source": width_source_counts,
            "slope_source": slope_source_counts,
            "outlet_method": outlet_method_counts,
        },
        "diagnostics": diagnostics,
        "subcatchments": detail_rows,
        "outputs": {
            "builder_csv": str(args.out_csv),
        },
    }
    write_json(args.out_json, report)

    print(
        json.dumps(
            {
                "ok": True,
                "out_csv": str(args.out_csv),
                "out_json": str(args.out_json),
                "subcatchment_count": len(detail_rows),
                "network_node_count": len(node_index),
                "dem_stats_count": len(dem_stats_index),
                "dem_stats_matched_count": len(used_dem_ids),
                "diagnostic_count": len(diagnostics),
            },
            indent=2,
        )
    )


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
from __future__ import annotations

import argparse
import json
import shutil
from pathlib import Path
from typing import Any

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.transform import xy
from shapely.geometry import LineString, Point

from plot_qgis_standard_layers import plot_overview


SHAPEFILE_SIDECARS = (".shp", ".shx", ".dbf", ".prj", ".cpg")


def copy_shapefile(src_shp: Path, dst_shp: Path) -> list[str]:
    if src_shp.suffix.lower() != ".shp":
        raise ValueError(f"Expected .shp source: {src_shp}")
    if not src_shp.exists():
        raise FileNotFoundError(src_shp)
    dst_shp.parent.mkdir(parents=True, exist_ok=True)
    copied: list[str] = []
    for suffix in SHAPEFILE_SIDECARS:
        src = src_shp.with_suffix(suffix)
        if src.exists():
            dst = dst_shp.with_suffix(suffix)
            shutil.copy2(src, dst)
            copied.append(str(dst))
    return copied


def copy_raster(src: Path, dst: Path) -> list[str]:
    if not src.exists():
        raise FileNotFoundError(src)
    dst.parent.mkdir(parents=True, exist_ok=True)
    shutil.copy2(src, dst)
    copied = [str(dst)]
    aux = src.with_name(src.name + ".aux.xml")
    if aux.exists():
        dst_aux = dst.with_name(dst.name + ".aux.xml")
        shutil.copy2(aux, dst_aux)
        copied.append(str(dst_aux))
    return copied


def clear_layer_stem(folder: Path, stem: str) -> None:
    for path in folder.glob(f"{stem}.*"):
        path.unlink()


def derive_flow_and_outfall(*, stream: Path, accumulation: Path, flow_shp: Path, outfall_shp: Path) -> dict[str, Any]:
    clear_layer_stem(flow_shp.parent, flow_shp.stem)
    clear_layer_stem(outfall_shp.parent, outfall_shp.stem)

    with rasterio.open(stream) as stream_ds, rasterio.open(accumulation) as acc_ds:
        stream_arr = stream_ds.read(1, masked=True)
        acc_arr = acc_ds.read(1, masked=True).astype(float)
        transform = stream_ds.transform
        crs = stream_ds.crs

        stream_mask = (~stream_arr.mask) & (stream_arr.filled(0) > 0)
        rows, cols = np.where(stream_mask)
        stream_cells = set(zip(rows.tolist(), cols.tolist()))

        geometries: list[LineString] = []
        records: list[dict[str, Any]] = []
        outfall_cell: tuple[int, int] | None = None
        outfall_acc = float("-inf")

        for r, c in zip(rows.tolist(), cols.tolist()):
            current_acc = float(acc_arr[r, c]) if not acc_arr.mask[r, c] else float("-inf")
            if current_acc > outfall_acc:
                outfall_acc = current_acc
                outfall_cell = (r, c)

            candidates: list[tuple[float, int, int]] = []
            for dr in (-1, 0, 1):
                for dc in (-1, 0, 1):
                    if dr == 0 and dc == 0:
                        continue
                    nr, nc = r + dr, c + dc
                    if (nr, nc) in stream_cells and not acc_arr.mask[nr, nc]:
                        neighbor_acc = float(acc_arr[nr, nc])
                        if neighbor_acc > current_acc:
                            candidates.append((neighbor_acc, nr, nc))

            if not candidates:
                continue

            _, nr, nc = max(candidates)
            x1, y1 = xy(transform, r, c, offset="center")
            x2, y2 = xy(transform, nr, nc, offset="center")
            geometries.append(LineString([(x1, y1), (x2, y2)]))
            records.append(
                {
                    "from_row": int(r),
                    "from_col": int(c),
                    "to_row": int(nr),
                    "to_col": int(nc),
                    "acc_from": current_acc,
                    "acc_to": float(acc_arr[nr, nc]),
                }
            )

        flow = gpd.GeoDataFrame(records, geometry=geometries, crs=crs)
        flow.to_file(flow_shp)

        if outfall_cell is None:
            outfall = gpd.GeoDataFrame([{"role": "outfall", "accum": None}], geometry=[Point()], crs=crs)
        else:
            r, c = outfall_cell
            x, y = xy(transform, r, c, offset="center")
            outfall = gpd.GeoDataFrame(
                [{"role": "outfall", "row": int(r), "col": int(c), "accum": float(outfall_acc)}],
                geometry=[Point(x, y)],
                crs=crs,
            )
        outfall.to_file(outfall_shp)

    return {
        "flow_segments": int(len(records)),
        "outfalls": 1 if outfall_cell is not None else 0,
        "outfall_accumulation": None if outfall_cell is None else float(outfall_acc),
    }


def package_final_layers(args: argparse.Namespace) -> dict[str, Any]:
    final_dir: Path = args.final_dir
    final_dir.mkdir(parents=True, exist_ok=True)

    subcatchments_dst = final_dir / "subcatchments.shp"
    flow_dst = final_dir / "flow.shp"
    slope_dst = final_dir / "slope_percent.tif"
    outfall_dst = final_dir / "outfall.shp"
    overview_dst = final_dir / "overview.png"
    manifest_dst = final_dir / "manifest.json"

    for stem in ("subcatchments", "flow", "outfall"):
        clear_layer_stem(final_dir, stem)
    for path in (slope_dst, slope_dst.with_name(slope_dst.name + ".aux.xml"), overview_dst, manifest_dst):
        if path.exists():
            path.unlink()

    copied_subcatchments = copy_shapefile(args.subcatchments, subcatchments_dst)
    copied_slope = copy_raster(args.slope, slope_dst)
    flow_summary = derive_flow_and_outfall(
        stream=args.stream,
        accumulation=args.accumulation,
        flow_shp=flow_dst,
        outfall_shp=outfall_dst,
    )

    if not args.no_overview:
        plot_overview(
            slope=slope_dst,
            subcatchments=subcatchments_dst,
            flow=flow_dst,
            outfall=outfall_dst,
            out_png=overview_dst,
            title=args.title or "",
        )

    subcatchments = gpd.read_file(subcatchments_dst)
    manifest = {
        "ok": True,
        "case_id": args.case_id,
        "final_layers": {
            "subcatchments": str(subcatchments_dst),
            "flow": str(flow_dst),
            "slope": str(slope_dst),
            "outfall": str(outfall_dst),
            "overview": None if args.no_overview else str(overview_dst),
        },
        "sources": {
            "subcatchments": str(args.subcatchments),
            "stream": str(args.stream),
            "accumulation": str(args.accumulation),
            "slope": str(args.slope),
        },
        "counts": {
            "subcatchments": int(len(subcatchments)),
            **flow_summary,
        },
        "copied_files": {
            "subcatchments": copied_subcatchments,
            "slope": copied_slope,
        },
        "crs": str(subcatchments.crs) if subcatchments.crs else None,
        "note": "Audit artifacts remain in the parent run directories; this folder contains only user-facing GIS/SWMM layers.",
    }
    manifest_dst.write_text(json.dumps(manifest, indent=2), encoding="utf-8")
    return manifest


def main() -> None:
    parser = argparse.ArgumentParser(description="Package QGIS/GRASS watershed outputs into a clean final_layers folder.")
    parser.add_argument("--case-id", required=True)
    parser.add_argument("--subcatchments", type=Path, required=True)
    parser.add_argument("--stream", type=Path, required=True)
    parser.add_argument("--accumulation", type=Path, required=True)
    parser.add_argument("--slope", type=Path, required=True)
    parser.add_argument("--final-dir", type=Path, required=True)
    parser.add_argument("--title", default="")
    parser.add_argument("--no-overview", action="store_true")
    args = parser.parse_args()
    print(json.dumps(package_final_layers(args), indent=2))


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
from __future__ import annotations

import argparse
import csv
import json
import os
import re
import shutil
import subprocess
import sys
from pathlib import Path
from typing import Any


REPO_ROOT = Path(__file__).resolve().parents[3]
GIS_DIR = Path(__file__).resolve().parents[1]

# Sibling-skill seam (issue #246): export-swmm-intermediates and
# import-drainage-assets subprocess-shell into swmm-params/swmm-network
# scripts. Resolve the skills/ root through --skills-root >
# AISWMM_SKILLS_ROOT env var > the repo-relative default, so a
# relocated/standalone deployment can point elsewhere without changing
# default behavior.
DEFAULT_SKILLS_ROOT = REPO_ROOT / "skills"
SKILLS_ROOT_ENV = "AISWMM_SKILLS_ROOT"

DEFAULT_QGIS_PROCESS = "/Applications/QGIS-final-4_0_2.app/Contents/MacOS/qgis_process"
DEFAULT_PROJ_LIB = "/Applications/QGIS-final-4_0_2.app/Contents/Resources/qgis/proj"
DEFAULT_GISBASE = "/Applications/GRASS-8.4.app/Contents/Resources"


def resolve_skills_root(cli_value: Path | None) -> Path:
    """Resolve the skills/ root directory holding sibling skills (e.g. swmm-params, swmm-network).

    Precedence: --skills-root flag > AISWMM_SKILLS_ROOT env var > default
    (repo-relative skills/ next to this skill's checkout).
    """
    if cli_value is not None:
        return cli_value
    env_value = os.environ.get(SKILLS_ROOT_ENV)
    if env_value:
        return Path(env_value)
    return DEFAULT_SKILLS_ROOT


def load_json(path: Path) -> Any:
    return json.loads(path.read_text(encoding="utf-8"))


def write_json(path: Path, obj: Any) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(json.dumps(obj, indent=2), encoding="utf-8")


def write_csv(path: Path, rows: list[dict[str, Any]], headers: list[str]) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", encoding="utf-8", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=headers)
        writer.writeheader()
        writer.writerows(rows)


def run_python(args: list[str]) -> dict[str, Any]:
    proc = subprocess.run(
        [sys.executable, *args],
        cwd=REPO_ROOT,
        check=True,
        capture_output=True,
        text=True,
    )
    text = proc.stdout.strip()
    if not text:
        return {"stdout": ""}
    try:
        return json.loads(text)
    except json.JSONDecodeError:
        return {"stdout": text}


def qgis_env(*, proj_lib: Path | None, gisbase: Path | None) -> dict[str, str]:
    env = os.environ.copy()
    if proj_lib:
        env["PROJ_LIB"] = str(proj_lib)
    if gisbase:
        env["GISBASE"] = str(gisbase)
    return env


def run_qgis(
    qgis_process: Path,
    algorithm: str,
    params: list[tuple[str, Any]],
    *,
    env: dict[str, str],
    cwd: Path,
    audit: list[dict[str, Any]],
) -> str:
    cmd = [str(qgis_process), "run", algorithm, "--"]
    cmd.extend(f"{key}={value}" for key, value in params if value is not None)
    proc = subprocess.run(cmd, cwd=cwd, env=env, text=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    audit.append(
        {
            "algorithm": algorithm,
            "cmd": cmd,
            "returncode": proc.returncode,
            "stdout_tail": proc.stdout[-4000:],
            "stderr_tail": proc.stderr[-4000:],
        }
    )
    if proc.returncode != 0:
        raise RuntimeError(f"QGIS algorithm failed rc={proc.returncode}: {algorithm}\n{proc.stderr}")
    return proc.stdout


def layer_sidecars(path: Path) -> list[str]:
    if path.suffix.lower() != ".shp":
        return []
    sidecars = []
    for suffix in [".shx", ".dbf", ".prj"]:
        candidate = path.with_suffix(suffix)
        if candidate.exists():
            sidecars.append(str(candidate))
    return sidecars


def read_crs_hint(path: Path) -> dict[str, Any]:
    suffix = path.suffix.lower()
    if suffix == ".shp":
        prj = path.with_suffix(".prj")
        if prj.exists():
            return {"source": str(prj), "kind": "wkt", "text": prj.read_text(encoding="utf-8", errors="ignore").strip()}
        return {"source": None, "kind": "missing_prj", "text": ""}
    if suffix in {".geojson", ".json"}:
        try:
            obj = load_json(path)
        except Exception:
            return {"source": str(path), "kind": "unknown_json", "text": ""}
        crs = obj.get("crs")
        if crs:
            return {"source": str(path), "kind": "geojson_crs", "text": json.dumps(crs, sort_keys=True)}
        return {"source": str(path), "kind": "geojson_default_or_unspecified", "text": ""}
    return {"source": str(path), "kind": "not_inspected", "text": ""}


def _epsg_from_wkt(wkt: str) -> str | None:
    """Return 'EPSG:NNNNN' for the OUTER CRS of a WKT string, or None.

    WKT nests authority nodes: the spheroid/datum/unit each carry their
    own ``AUTHORITY["EPSG",...]`` (WKT1) or ``ID["EPSG",...]`` (WKT2)
    and the outer CRS's node is written LAST. Taking the FIRST match
    returned the spheroid's code (e.g. 7019/7030) instead of the CRS,
    so a WGS84/UTM layer read as "EPSG:7019" (found 2026-08-08). The
    last match is the outer CRS in both WKT1 and WKT2.
    """
    matches = re.findall(r'AUTHORITY\["EPSG",\s*"(\d+)"\]', wkt)
    if not matches:
        matches = re.findall(r'ID\["EPSG",\s*(\d+)\]', wkt)
    if matches:
        return f"EPSG:{matches[-1]}"
    return None


def get_layer_epsg(path: Path) -> str | None:
    """Return 'EPSG:NNNNN' for a layer file, or None if undetermined."""
    suffix = path.suffix.lower()
    if suffix in {".tif", ".tiff"}:
        try:
            result = subprocess.run(
                ["gdalinfo", "-json", str(path)],
                capture_output=True, text=True, timeout=30,
            )
            if result.returncode == 0:
                info = json.loads(result.stdout)
                wkt = info.get("coordinateSystem", {}).get("wkt", "")
                return _epsg_from_wkt(wkt)
        except Exception:
            pass
        return None
    if suffix == ".shp":
        prj = path.with_suffix(".prj")
        if prj.exists():
            wkt = prj.read_text(encoding="utf-8", errors="ignore")
            return _epsg_from_wkt(wkt)
        return None
    if suffix in {".geojson", ".json"}:
        try:
            obj = load_json(path)
            name = obj.get("crs", {}).get("properties", {}).get("name", "")
            m = re.search(r'EPSG[::]+(\d+)', name)
            if m:
                return f"EPSG:{m.group(1)}"
        except Exception:
            pass
        return None
    return None


def resolve_target_epsg(target_crs_ref: str) -> str | None:
    """Return 'EPSG:NNNNN' from an EPSG string like 'EPSG:32610' or a file path."""
    epsg_match = re.match(r'epsg:(\d+)', target_crs_ref, re.IGNORECASE)
    if epsg_match:
        return f"EPSG:{epsg_match.group(1)}"
    p = Path(target_crs_ref)
    return get_layer_epsg(p) if p.exists() else None


def build_layers_manifest(layers: dict[str, str | None], out_path: Path) -> dict[str, Any]:
    records = []
    issues = []
    for role, raw_path in layers.items():
        if raw_path is None:
            continue
        path = Path(raw_path)
        record = {
            "role": role,
            "path": str(path),
            "exists": path.exists(),
            "suffix": path.suffix.lower(),
            "sidecars": layer_sidecars(path),
            "crs": read_crs_hint(path) if path.exists() else None,
        }
        if not path.exists():
            issues.append({"severity": "error", "role": role, "message": f"Layer path does not exist: {path}"})
        if path.suffix.lower() == ".shp":
            missing = [suffix for suffix in [".shx", ".dbf"] if not path.with_suffix(suffix).exists()]
            for suffix in missing:
                issues.append({"severity": "error", "role": role, "message": f"Missing shapefile sidecar {suffix}: {path}"})
            if not path.with_suffix(".prj").exists():
                issues.append({"severity": "warning", "role": role, "message": f"Missing shapefile CRS sidecar .prj: {path}"})
        records.append(record)

    manifest = {
        "ok": not any(issue["severity"] == "error" for issue in issues),
        "skill": "swmm-gis",
        "adapter": "qgis_data_prep",
        "layers": records,
        "issue_count": len(issues),
        "issues": issues,
    }
    write_json(out_path, manifest)
    return manifest


def _canonical_crs_key(text: str) -> str:
    """Reduce a CRS description to a comparable identity.

    Two exports of the SAME CRS routinely differ as raw text: ESRI vs
    OGC WKT1 spelling ("WGS_1984_UTM_Zone_10N" vs "WGS 84 / UTM zone
    10N"), WKT1 vs WKT2, parameter formatting. Comparing raw strings
    flagged such pairs as "different CRS" and demanded a reprojection
    that was not needed (found 2026-08-08). Canonicalize to the EPSG
    code whenever one is present (WKT authority nodes, GeoJSON
    ``urn:ogc:def:crs:EPSG::N`` / ``EPSG:N`` names); layers whose text
    carries no EPSG identity fall back to whitespace-normalized text,
    where a mismatch stays a report-worthy difference.
    """
    epsg = _epsg_from_wkt(text)
    if epsg:
        return epsg
    m = re.search(r'EPSG[:]{1,2}(\d+)', text)
    if m:
        return f"EPSG:{m.group(1)}"
    return f"wkt:{text}"


def validate_crs(layers_manifest: dict[str, Any], out_path: Path) -> dict[str, Any]:
    comparable = []
    issues = []
    for layer in layers_manifest["layers"]:
        crs = layer.get("crs") or {}
        text = " ".join(str(crs.get("text") or "").split())
        record = {"role": layer["role"], "path": layer["path"], "kind": crs.get("kind"), "text": text}
        record["canonical"] = _canonical_crs_key(text) if text else None
        comparable.append(record)
        if not text and crs.get("kind") not in {"geojson_default_or_unspecified", "not_inspected"}:
            issues.append({"severity": "warning", "role": layer["role"], "message": "CRS could not be confirmed from the source layer."})

    explicit = [item for item in comparable if item["text"]]
    unique = sorted({item["canonical"] for item in explicit})
    if len(unique) > 1:
        issues.append(
            {
                "severity": "error",
                "message": "Input layers expose different CRS definitions. Reproject in QGIS before export.",
                "explicit_crs_count": len(unique),
                "distinct_crs": unique,
            }
        )

    report = {
        "ok": not any(issue["severity"] == "error" for issue in issues),
        "skill": "swmm-gis",
        "adapter": "qgis_data_prep",
        "crs_records": comparable,
        "issue_count": len(issues),
        "issues": issues,
        "assumption": "GeoJSON without explicit CRS is treated as already exported in the project CRS.",
    }
    write_json(out_path, report)
    return report


def is_raster(path: Path) -> bool:
    return path.suffix.lower() in {".tif", ".tiff", ".img", ".vrt", ".asc"}


def is_vector(path: Path) -> bool:
    return path.suffix.lower() in {".shp", ".geojson", ".json", ".gpkg"}


def normalize_vector_layer(
    *,
    qgis_process: Path,
    source: Path,
    boundary: Path,
    target_crs: str,
    out_path: Path,
    work_dir: Path,
    env: dict[str, str],
    audit: list[dict[str, Any]],
    skip_reproject: bool = False,
) -> dict[str, str]:
    work_dir.mkdir(parents=True, exist_ok=True)
    if skip_reproject:
        vector_for_clip = source
        reprojected_str = str(source)
    else:
        reprojected = work_dir / f"{out_path.stem}_reprojected.shp"
        run_qgis(
            qgis_process,
            "native:reprojectlayer",
            [
                ("INPUT", source),
                ("TARGET_CRS", target_crs),
                ("OUTPUT", reprojected),
            ],
            env=env,
            cwd=REPO_ROOT,
            audit=audit,
        )
        vector_for_clip = reprojected
        reprojected_str = str(reprojected)

    run_qgis(
        qgis_process,
        "native:clip",
        [
            ("INPUT", vector_for_clip),
            ("OVERLAY", boundary),
            ("OUTPUT", out_path),
        ],
        env=env,
        cwd=REPO_ROOT,
        audit=audit,
    )
    return {"reprojected": reprojected_str, "clipped": str(out_path)}


def normalize_raster_layer(
    *,
    qgis_process: Path,
    source: Path,
    boundary: Path,
    target_crs: str,
    out_path: Path,
    work_dir: Path,
    resampling: int,
    target_resolution: float | None,
    env: dict[str, str],
    audit: list[dict[str, Any]],
    skip_reproject: bool = False,
) -> dict[str, str]:
    work_dir.mkdir(parents=True, exist_ok=True)
    if skip_reproject:
        raster_for_clip = source
        reprojected_str = str(source)
    else:
        reprojected = work_dir / f"{out_path.stem}_reprojected.tif"
        warp_params: list[tuple[str, Any]] = [
            ("INPUT", source),
            ("TARGET_CRS", target_crs),
            ("RESAMPLING", resampling),
            ("MULTITHREADING", "true"),
            ("OUTPUT", reprojected),
        ]
        if target_resolution:
            warp_params.insert(3, ("TARGET_RESOLUTION", target_resolution))
        run_qgis(qgis_process, "gdal:warpreproject", warp_params, env=env, cwd=REPO_ROOT, audit=audit)
        raster_for_clip = reprojected
        reprojected_str = str(reprojected)

    clip_params: list[tuple[str, Any]] = [
        ("INPUT", raster_for_clip),
        ("MASK", boundary),
        ("TARGET_CRS", target_crs),
        ("CROP_TO_CUTLINE", "true"),
        ("KEEP_RESOLUTION", "true"),
        ("OUTPUT", out_path),
    ]
    run_qgis(qgis_process, "gdal:cliprasterbymasklayer", clip_params, env=env, cwd=REPO_ROOT, audit=audit)
    return {"reprojected": reprojected_str, "clipped": str(out_path)}


def normalize_layers(args: argparse.Namespace) -> dict[str, Any]:
    out_dir: Path = args.out_dir
    work_dir = out_dir / "_work"
    out_dir.mkdir(parents=True, exist_ok=True)
    commands: list[dict[str, Any]] = []
    env = qgis_env(proj_lib=args.proj_lib, gisbase=args.gisbase)

    sources = {
        "dem": args.dem,
        "boundary": args.boundary,
        "landuse": args.landuse,
        "soil": args.soil,
    }
    for role, path in sources.items():
        if not path.exists():
            raise FileNotFoundError(f"Missing {role}: {path}")

    target_crs = args.target_crs or str(args.boundary)
    target_epsg = resolve_target_epsg(target_crs)

    boundary_out = out_dir / "boundary.shp"
    boundary_epsg = get_layer_epsg(args.boundary)
    boundary_needs_reproject = not (target_epsg and boundary_epsg and boundary_epsg == target_epsg)
    if boundary_needs_reproject:
        run_qgis(
            args.qgis_process,
            "native:reprojectlayer",
            [
                ("INPUT", args.boundary),
                ("TARGET_CRS", target_crs),
                ("OUTPUT", boundary_out),
            ],
            env=env,
            cwd=REPO_ROOT,
            audit=commands,
        )
        commands.append({"action": "reprojectlayer", "role": "boundary", "source": str(args.boundary), "output": str(boundary_out)})
    else:
        shutil.copy2(args.boundary, boundary_out)
        for sc in layer_sidecars(args.boundary):
            shutil.copy2(sc, boundary_out.with_suffix(Path(sc).suffix))
        commands.append({"action": "copy_same_crs", "role": "boundary", "source": str(args.boundary), "output": str(boundary_out)})
    boundary_crs_ref = str(boundary_out) if not args.target_crs else args.target_crs

    outputs: dict[str, str] = {"boundary": str(boundary_out)}
    stage_outputs: dict[str, Any] = {"boundary": {"reprojected": str(boundary_out), "clipped": str(boundary_out)}}

    dem_out = out_dir / "dem.tif"
    dem_epsg = get_layer_epsg(args.dem)
    dem_skip = bool(target_epsg and dem_epsg and dem_epsg == target_epsg)
    stage_outputs["dem"] = normalize_raster_layer(
        qgis_process=args.qgis_process,
        source=args.dem,
        boundary=boundary_out,
        target_crs=boundary_crs_ref,
        out_path=dem_out,
        work_dir=work_dir,
        resampling=args.dem_resampling,
        target_resolution=args.target_resolution,
        env=env,
        audit=commands,
        skip_reproject=dem_skip,
    )
    outputs["dem"] = str(dem_out)

    for role, source in [("landuse", args.landuse), ("soil", args.soil)]:
        src_epsg = get_layer_epsg(source)
        skip = bool(target_epsg and src_epsg and src_epsg == target_epsg)
        if is_raster(source):
            out_path = out_dir / f"{role}.tif"
            stage_outputs[role] = normalize_raster_layer(
                qgis_process=args.qgis_process,
                source=source,
                boundary=boundary_out,
                target_crs=boundary_crs_ref,
                out_path=out_path,
                work_dir=work_dir,
                resampling=args.categorical_resampling,
                target_resolution=args.target_resolution,
                env=env,
                audit=commands,
                skip_reproject=skip,
            )
        elif is_vector(source):
            out_path = out_dir / f"{role}.shp"
            stage_outputs[role] = normalize_vector_layer(
                qgis_process=args.qgis_process,
                source=source,
                boundary=boundary_out,
                target_crs=boundary_crs_ref,
                out_path=out_path,
                work_dir=work_dir,
                env=env,
                audit=commands,
                skip_reproject=skip,
            )
        else:
            raise ValueError(f"Unsupported {role} layer type: {source}")
        outputs[role] = str(out_path)

    manifest = {
        "ok": True,
        "skill": "swmm-gis",
        "adapter": "qgis_normalize_layers",
        "qgis_process": str(args.qgis_process),
        "target_crs": target_crs,
        "target_crs_policy": "explicit target CRS" if args.target_crs else "boundary CRS",
        "target_resolution": args.target_resolution,
        "resampling": {
            "dem": args.dem_resampling,
            "categorical": args.categorical_resampling,
        },
        "inputs": {role: str(path) for role, path in sources.items()},
        "outputs": outputs,
        "stage_outputs": stage_outputs,
        "processing_commands": commands,
        "evidence_boundary": (
            "This step standardizes CRS and clips input GIS layers to the study boundary. "
            "It does not delineate streams, choose pour points, or run SWMM."
        ),
    }
    write_json(out_dir / "qgis_normalized_layers_manifest.json", manifest)
    return manifest


def extract_overlay_tables(
    subcatchments_geojson: Path,
    *,
    landuse_field: str,
    soil_field: str,
    id_field: str,
    out_landuse_csv: Path,
    out_soil_csv: Path,
) -> dict[str, Any]:
    obj = load_json(subcatchments_geojson)
    if obj.get("type") != "FeatureCollection":
        raise ValueError(f"Expected FeatureCollection: {subcatchments_geojson}")

    landuse_rows: list[dict[str, str]] = []
    soil_rows: list[dict[str, str]] = []
    issues = []
    seen: set[str] = set()
    for idx, feature in enumerate(obj.get("features") or [], start=1):
        props = feature.get("properties") or {}
        sid = str(props.get(id_field) or feature.get("id") or "").strip()
        if not sid:
            raise ValueError(f"Feature {idx} missing subcatchment id field '{id_field}'")
        if sid in seen:
            raise ValueError(f"Duplicate subcatchment id: {sid}")
        seen.add(sid)

        landuse = str(props.get(landuse_field) or "DEFAULT").strip() or "DEFAULT"
        soil = str(props.get(soil_field) or "-").strip() or "-"
        if landuse == "DEFAULT":
            issues.append({"severity": "warning", "id": sid, "message": f"Missing land use field '{landuse_field}', using DEFAULT."})
        if soil == "-":
            issues.append({"severity": "warning", "id": sid, "message": f"Missing soil field '{soil_field}', using fallback '-'."})
        landuse_rows.append({"subcatchment_id": sid, "landuse_class": landuse})
        soil_rows.append({"subcatchment_id": sid, "soil_texture": soil})

    write_csv(out_landuse_csv, landuse_rows, ["subcatchment_id", "landuse_class"])
    write_csv(out_soil_csv, soil_rows, ["subcatchment_id", "soil_texture"])

    return {
        "ok": True,
        "landuse_csv": str(out_landuse_csv),
        "soil_csv": str(out_soil_csv),
        "subcatchment_count": len(seen),
        "issue_count": len(issues),
        "issues": issues,
    }


def copy_network_and_qa(
    network_json: Path,
    out_network_json: Path,
    out_qa_json: Path,
    *,
    network_skill_dir: Path = DEFAULT_SKILLS_ROOT / "swmm-network",
) -> dict[str, Any]:
    out_network_json.parent.mkdir(parents=True, exist_ok=True)
    shutil.copyfile(network_json, out_network_json)
    qa = run_python(
        [
            str(network_skill_dir / "scripts/network_qa.py"),
            str(out_network_json),
            "--report-json",
            str(out_qa_json),
        ]
    )
    return {"network_json": str(out_network_json), "network_qa_json": str(out_qa_json), "qa": qa}


def export_swmm_intermediates(args: argparse.Namespace) -> dict[str, Any]:
    run_dir: Path = args.run_dir
    raw_dir = run_dir / "00_raw"
    gis_dir = run_dir / "01_gis"
    params_dir = run_dir / "02_params"
    network_dir = run_dir / "04_network"

    skills_root = resolve_skills_root(args.skills_root)
    params_skill_dir = skills_root / "swmm-params"
    network_skill_dir = skills_root / "swmm-network"

    layers = {
        "dem": str(args.dem) if args.dem else None,
        "subcatchments": str(args.subcatchments_geojson),
        "landuse": str(args.landuse_layer) if args.landuse_layer else None,
        "soil": str(args.soil_layer) if args.soil_layer else None,
        # outlet is a node ID string, not a file path — exclude from file validation
        "rainfall": str(args.rainfall) if args.rainfall else None,
        "network": str(args.network_json),
    }
    layers_manifest = build_layers_manifest(layers, raw_dir / "qgis_layers_manifest.json")
    crs_report = validate_crs(layers_manifest, raw_dir / "qgis_crs_report.json")
    if not layers_manifest["ok"]:
        raise ValueError(f"Layer validation failed. See {raw_dir / 'qgis_layers_manifest.json'}")
    if args.strict_crs and not crs_report["ok"]:
        raise ValueError(f"CRS validation failed. See {raw_dir / 'qgis_crs_report.json'}")

    subcatchments_work = gis_dir / "subcatchments.geojson"
    subcatchments_work.parent.mkdir(parents=True, exist_ok=True)
    shutil.copyfile(args.subcatchments_geojson, subcatchments_work)

    preprocess = run_python(
        [
            str(GIS_DIR / "scripts/preprocess_subcatchments.py"),
            "--subcatchments-geojson",
            str(subcatchments_work),
            "--network-json",
            str(args.network_json),
            "--out-csv",
            str(gis_dir / "subcatchments.csv"),
            "--out-json",
            str(gis_dir / "subcatchments.json"),
            "--id-field",
            args.id_field,
            "--outlet-hint-field",
            args.outlet_hint_field,
            "--default-rain-gage",
            args.default_rain_gage,
            "--default-slope-pct",
            str(args.default_slope_pct),
            "--min-slope-pct",
            str(args.min_slope_pct),
            "--min-width-m",
            str(args.min_width_m),
        ]
    )

    overlay = extract_overlay_tables(
        subcatchments_work,
        landuse_field=args.landuse_field,
        soil_field=args.soil_field,
        id_field=args.id_field,
        out_landuse_csv=params_dir / "landuse.csv",
        out_soil_csv=params_dir / "soil.csv",
    )

    landuse = run_python(
        [
            str(params_skill_dir / "scripts/landuse_to_swmm_params.py"),
            "--input",
            str(params_dir / "landuse.csv"),
            "--output",
            str(params_dir / "landuse.json"),
        ]
    )
    soil = run_python(
        [
            str(params_skill_dir / "scripts/soil_to_greenampt.py"),
            "--input",
            str(params_dir / "soil.csv"),
            "--output",
            str(params_dir / "soil.json"),
        ]
    )
    merged = run_python(
        [
            str(params_skill_dir / "scripts/merge_swmm_params.py"),
            "--landuse-json",
            str(params_dir / "landuse.json"),
            "--soil-json",
            str(params_dir / "soil.json"),
            "--output",
            str(params_dir / "merged_params.json"),
        ]
    )

    network = copy_network_and_qa(
        args.network_json,
        network_dir / "network.json",
        network_dir / "network_qa.json",
        network_skill_dir=network_skill_dir,
    )

    report = {
        "ok": True,
        "skill": "swmm-gis",
        "adapter": "qgis_data_prep",
        "case_id": args.case_id,
        "run_dir": str(run_dir),
        "workflow_boundary": (
            "MVP expects QGIS to provide delineated subcatchment polygons and overlay attributes; "
            "this script validates sources and exports Agentic SWMM-ready intermediates."
        ),
        "outputs": {
            "layers_manifest": str(raw_dir / "qgis_layers_manifest.json"),
            "crs_report": str(raw_dir / "qgis_crs_report.json"),
            "subcatchments_geojson": str(subcatchments_work),
            "subcatchments_csv": str(gis_dir / "subcatchments.csv"),
            "subcatchments_json": str(gis_dir / "subcatchments.json"),
            "landuse_csv": str(params_dir / "landuse.csv"),
            "soil_csv": str(params_dir / "soil.csv"),
            "landuse_json": str(params_dir / "landuse.json"),
            "soil_json": str(params_dir / "soil.json"),
            "merged_params_json": str(params_dir / "merged_params.json"),
            "network_json": str(network_dir / "network.json"),
            "network_qa_json": str(network_dir / "network_qa.json"),
            "qgis_export_manifest": str(run_dir / "qgis_export_manifest.json"),
        },
        "stage_summaries": {
            "layers": layers_manifest,
            "crs": crs_report,
            "preprocess": preprocess,
            "overlay": overlay,
            "landuse": landuse,
            "soil": soil,
            "merged_params": merged,
            "network": network,
        },
    }
    write_json(run_dir / "qgis_export_manifest.json", report)
    return report


def add_common_export_args(ap: argparse.ArgumentParser) -> None:
    ap.add_argument("--run-dir", type=Path, required=True)
    ap.add_argument("--case-id", default="qgis-case")
    ap.add_argument("--subcatchments-geojson", type=Path, required=True)
    ap.add_argument("--network-json", type=Path, required=True)
    ap.add_argument("--dem", type=Path, default=None)
    ap.add_argument("--landuse-layer", type=Path, default=None)
    ap.add_argument("--soil-layer", type=Path, default=None)
    ap.add_argument("--outlet", type=Path, default=None)
    ap.add_argument("--rainfall", type=Path, default=None)
    ap.add_argument("--id-field", default="subcatchment_id")
    ap.add_argument("--outlet-hint-field", default="outlet_hint")
    ap.add_argument("--landuse-field", default="landuse_class")
    ap.add_argument("--soil-field", default="soil_texture")
    ap.add_argument("--default-rain-gage", default="RG1")
    ap.add_argument("--default-slope-pct", type=float, default=1.0)
    ap.add_argument("--min-slope-pct", type=float, default=0.1)
    ap.add_argument("--min-width-m", type=float, default=10.0)
    ap.add_argument("--strict-crs", action="store_true")
    ap.add_argument(
        "--skills-root",
        type=Path,
        default=None,
        help=(
            "Root directory containing sibling skills swmm-params/swmm-network. "
            f"Overrides: flag > {SKILLS_ROOT_ENV} env var > default '<repo>/skills'."
        ),
    )


def main() -> None:
    parser = argparse.ArgumentParser(description="QGIS-oriented raw GIS to Agentic SWMM intermediate exporter.")
    sub = parser.add_subparsers(dest="command", required=True)

    load_ap = sub.add_parser("load-layers", help="Validate raw/QGIS-exported layer paths and sidecars.")
    load_ap.add_argument("--out", type=Path, required=True)
    load_ap.add_argument("--dem", type=Path)
    load_ap.add_argument("--boundary", type=Path)
    load_ap.add_argument("--subcatchments", type=Path)
    load_ap.add_argument("--landuse", type=Path)
    load_ap.add_argument("--soil", type=Path)
    load_ap.add_argument("--outlet", type=Path)
    load_ap.add_argument("--rainfall", type=Path)
    load_ap.add_argument("--network", type=Path)

    crs_ap = sub.add_parser("validate-crs", help="Check source CRS hints from a layer manifest.")
    crs_ap.add_argument("--layers-manifest", type=Path, required=True)
    crs_ap.add_argument("--out", type=Path, required=True)

    normalize_ap = sub.add_parser("normalize-layers", help="Reproject DEM/boundary/landuse/soil to one CRS and clip them by the boundary.")
    normalize_ap.add_argument("--dem", type=Path, required=True)
    normalize_ap.add_argument("--boundary", type=Path, required=True)
    normalize_ap.add_argument("--landuse", type=Path, required=True)
    normalize_ap.add_argument("--soil", type=Path, required=True)
    normalize_ap.add_argument("--out-dir", type=Path, required=True)
    normalize_ap.add_argument("--target-crs", help="Target CRS auth id, WKT/PROJ string, or layer path. Defaults to the boundary CRS.")
    normalize_ap.add_argument("--target-resolution", type=float, help="Optional target raster resolution in target CRS units.")
    normalize_ap.add_argument("--dem-resampling", type=int, default=1, help="GDAL resampling enum for DEM reprojection; default 1 is bilinear.")
    normalize_ap.add_argument("--categorical-resampling", type=int, default=0, help="GDAL resampling enum for categorical rasters; default 0 is nearest.")
    normalize_ap.add_argument("--qgis-process", type=Path, default=Path(DEFAULT_QGIS_PROCESS))
    normalize_ap.add_argument("--proj-lib", type=Path, default=Path(DEFAULT_PROJ_LIB))
    normalize_ap.add_argument("--gisbase", type=Path, default=Path(DEFAULT_GISBASE))

    overlay_ap = sub.add_parser("overlay-landuse-soil", help="Extract subcatchment_id, landuse_class, and soil_texture from a QGIS overlay export.")
    overlay_ap.add_argument("--subcatchments-geojson", type=Path, required=True)
    overlay_ap.add_argument("--out-landuse-csv", type=Path, required=True)
    overlay_ap.add_argument("--out-soil-csv", type=Path, required=True)
    overlay_ap.add_argument("--id-field", default="subcatchment_id")
    overlay_ap.add_argument("--landuse-field", default="landuse_class")
    overlay_ap.add_argument("--soil-field", default="soil_texture")

    network_ap = sub.add_parser("import-drainage-assets", help="Copy a prepared network JSON into 04_network and run network QA.")
    network_ap.add_argument("--network-json", type=Path, required=True)
    network_ap.add_argument("--out-network-json", type=Path, required=True)
    network_ap.add_argument("--out-qa-json", type=Path, required=True)
    network_ap.add_argument(
        "--skills-root",
        type=Path,
        default=None,
        help=(
            "Root directory containing the sibling swmm-network skill. "
            f"Overrides: flag > {SKILLS_ROOT_ENV} env var > default '<repo>/skills'."
        ),
    )

    export_ap = sub.add_parser("export-swmm-intermediates", help="Export 01_gis, 02_params, and 04_network artifacts from QGIS-prepared sources.")
    add_common_export_args(export_ap)

    args = parser.parse_args()
    if args.command == "load-layers":
        layers = {
            "dem": str(args.dem) if args.dem else None,
            "boundary": str(args.boundary) if args.boundary else None,
            "subcatchments": str(args.subcatchments) if args.subcatchments else None,
            "landuse": str(args.landuse) if args.landuse else None,
            "soil": str(args.soil) if args.soil else None,
            "outlet": str(args.outlet) if args.outlet else None,
            "rainfall": str(args.rainfall) if args.rainfall else None,
            "network": str(args.network) if args.network else None,
        }
        result = build_layers_manifest(layers, args.out)
    elif args.command == "validate-crs":
        result = validate_crs(load_json(args.layers_manifest), args.out)
    elif args.command == "normalize-layers":
        result = normalize_layers(args)
    elif args.command == "overlay-landuse-soil":
        result = extract_overlay_tables(
            args.subcatchments_geojson,
            landuse_field=args.landuse_field,
            soil_field=args.soil_field,
            id_field=args.id_field,
            out_landuse_csv=args.out_landuse_csv,
            out_soil_csv=args.out_soil_csv,
        )
    elif args.command == "import-drainage-assets":
        skills_root = resolve_skills_root(args.skills_root)
        result = copy_network_and_qa(
            args.network_json,
            args.out_network_json,
            args.out_qa_json,
            network_skill_dir=skills_root / "swmm-network",
        )
    elif args.command == "export-swmm-intermediates":
        result = export_swmm_intermediates(args)
    else:
        raise ValueError(f"Unknown command: {args.command}")

    print(json.dumps(result, indent=2))


if __name__ == "__main__":
    main()