{
  "meta": {
    "name": "basic-demo-network",
    "flow_units": "CMS"
  },
  "junctions": [
    {
      "id": "J1",
      "invert_elev": 100.0,
      "max_depth": 3.0,
      "init_depth": 0.0,
      "sur_depth": 0.0,
      "aponded": 0.0,
      "coordinates": { "x": 0.0, "y": 0.0 }
    },
    {
      "id": "J2",
      "invert_elev": 99.2,
      "max_depth": 3.0,
      "init_depth": 0.0,
      "sur_depth": 0.0,
      "aponded": 0.0,
      "coordinates": { "x": 100.0, "y": 0.0 }
    }
  ],
  "outfalls": [
    {
      "id": "OF1",
      "invert_elev": 98.5,
      "type": "FREE",
      "gated": false,
      "route_to": null,
      "coordinates": { "x": 200.0, "y": 0.0 }
    }
  ],
  "conduits": [
    {
      "id": "C1",
      "from_node": "J1",
      "to_node": "J2",
      "length": 100.0,
      "roughness": 0.013,
      "in_offset": 0.0,
      "out_offset": 0.0,
      "init_flow": 0.0,
      "max_flow": null,
      "xsection": {
        "shape": "CIRCULAR",
        "geom1": 0.6,
        "geom2": 0.0,
        "geom3": 0.0,
        "geom4": 0.0,
        "barrels": 1
      },
      "vertices": [
        { "x": 50.0, "y": 10.0 }
      ]
    },
    {
      "id": "C2",
      "from_node": "J2",
      "to_node": "OF1",
      "length": 100.0,
      "roughness": 0.013,
      "in_offset": 0.0,
      "out_offset": 0.0,
      "init_flow": 0.0,
      "max_flow": null,
      "xsection": {
        "shape": "CIRCULAR",
        "geom1": 0.75,
        "geom2": 0.0,
        "geom3": 0.0,
        "geom4": 0.0,
        "barrels": 1
      }
    }
  ]
}

{
  "meta": {
    "name": "import-demo-network",
    "source": "example-geojson"
  },
  "junctions": {
    "format": "geojson",
    "fields": {
      "id": "node_id",
      "invert_elev": "inv_el",
      "max_depth": "max_d",
      "init_depth": null,
      "sur_depth": null,
      "aponded": null
    }
  },
  "outfalls": {
    "format": "geojson",
    "fields": {
      "id": "node_id",
      "invert_elev": "inv_el",
      "type": "out_type",
      "stage_data": null,
      "gated": null,
      "route_to": null
    },
    "defaults": {
      "type": "FREE",
      "gated": false
    }
  },
  "conduits": {
    "format": "geojson",
    "fields": {
      "id": "link_id",
      "from_node": "from_id",
      "to_node": "to_id",
      "length": "len_m",
      "roughness": "n_val",
      "diameter": "diam_m",
      "in_offset": null,
      "out_offset": null,
      "init_flow": null,
      "max_flow": null,
      "shape": null,
      "geom2": null,
      "geom3": null,
      "geom4": null,
      "barrels": null
    },
    "defaults": {
      "shape": "CIRCULAR",
      "roughness": 0.013,
      "geom1": 0.5,
      "geom2": 0.0,
      "geom3": 0.0,
      "geom4": 0.0,
      "barrels": 1,
      "in_offset": 0.0,
      "out_offset": 0.0,
      "init_flow": 0.0,
      "max_flow": null
    }
  }
}

#!/usr/bin/env python3
"""Assign each SWMM subcatchment to a real network node as its outlet.

Without this step the subcatchments.csv produced by
``basin_shp_to_subcatchments`` carries the literal outfall as every
subcatchment's outlet, which means surface runoff dumps straight to the
outfall and the pipe network sits idle in the SWMM model. This tool
rewrites the ``outlet`` column so runoff actually enters the pipe
network at a sensible upstream node.

Three modes:

- ``nearest_junction`` (default): centroid of each subcatchment polygon
  is matched to the closest node listed in network.json (junctions
  + outfalls). Use when no manhole layer is available.
- ``nearest_catch_basin``: same, but the candidate node set is read
  from a separate manholes / catch-basin GeoJSON. Use when a richer
  catch-basin layer is available than what was inferred during
  network import.
- ``manual_lookup``: read a CSV with two columns
  ``subcatchment_id,outlet_node_id`` and apply the mapping verbatim.
  Use when the agent (or a human) wants to override.
"""
from __future__ import annotations

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

import geopandas as gpd
from shapely.geometry import Point, shape


MODES = ("nearest_junction", "nearest_catch_basin", "manual_lookup")


def _load_subcatchment_centroids(geojson_path: Path) -> list[tuple[str, Point]]:
    obj = json.loads(geojson_path.read_text(encoding="utf-8"))
    out: list[tuple[str, Point]] = []
    for feat in obj.get("features") or []:
        props = feat.get("properties") or {}
        sid = props.get("subcatchment_id")
        if sid is None:
            continue
        geom = shape(feat["geometry"])
        out.append((str(sid), geom.centroid))
    if not out:
        raise ValueError(f"no subcatchments with subcatchment_id field in {geojson_path}")
    return out


def _node_xy(node: dict) -> tuple[float, float] | None:
    """Read xy from either {coordinates: {x, y}} or top-level {x, y}."""
    coords = node.get("coordinates")
    if isinstance(coords, dict) and coords.get("x") is not None and coords.get("y") is not None:
        return float(coords["x"]), float(coords["y"])
    if node.get("x") is not None and node.get("y") is not None:
        return float(node["x"]), float(node["y"])
    return None


def _candidate_nodes_from_network(network_json_path: Path, include_outfalls: bool) -> list[tuple[str, Point]]:
    obj = json.loads(network_json_path.read_text(encoding="utf-8"))
    out: list[tuple[str, Point]] = []
    for j in obj.get("junctions") or []:
        nid = j.get("id")
        xy = _node_xy(j)
        if nid is None or xy is None:
            continue
        out.append((str(nid), Point(*xy)))
    if include_outfalls:
        for o in obj.get("outfalls") or []:
            nid = o.get("id")
            xy = _node_xy(o)
            if nid is None or xy is None:
                continue
            out.append((str(nid), Point(*xy)))
    if not out:
        raise ValueError(f"no candidate nodes in {network_json_path} (junctions+outfalls)")
    return out


def _candidate_nodes_from_geojson(path: Path, id_field: str) -> list[tuple[str, Point]]:
    gdf = gpd.read_file(path)
    if id_field not in gdf.columns:
        raise ValueError(f"id_field '{id_field}' not in {path} columns: {list(gdf.columns)}")
    out: list[tuple[str, Point]] = []
    for _, row in gdf.iterrows():
        nid = row[id_field]
        if nid is None:
            continue
        geom = row.geometry
        if geom is None:
            continue
        if geom.geom_type != "Point":
            geom = geom.centroid
        out.append((str(nid), geom))
    if not out:
        raise ValueError(f"no usable point candidates in {path}")
    return out


def _read_manual_lookup(path: Path) -> dict[str, str]:
    mapping: dict[str, str] = {}
    with path.open("r", encoding="utf-8", newline="") as f:
        reader = csv.DictReader(f)
        if not reader.fieldnames or "subcatchment_id" not in reader.fieldnames or "outlet_node_id" not in reader.fieldnames:
            raise ValueError(
                f"{path} must have headers 'subcatchment_id' and 'outlet_node_id'; "
                f"got {reader.fieldnames}"
            )
        for row in reader:
            mapping[str(row["subcatchment_id"]).strip()] = str(row["outlet_node_id"]).strip()
    if not mapping:
        raise ValueError(f"{path} has no rows")
    return mapping


def _nearest(point: Point, candidates: list[tuple[str, Point]]) -> tuple[str, float]:
    best_id = None
    best_dist = float("inf")
    for nid, pt in candidates:
        d = point.distance(pt)
        if d < best_dist:
            best_dist = d
            best_id = nid
    if best_id is None:
        raise RuntimeError("no candidates supplied")
    return best_id, float(best_dist)


def parse_args() -> argparse.Namespace:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--subcatchments-csv-in", required=True)
    ap.add_argument("--subcatchments-geojson", required=True)
    ap.add_argument("--out-csv", required=True)
    ap.add_argument("--mode", choices=MODES, default="nearest_junction")
    ap.add_argument("--network-json", default=None,
                    help="Required for mode=nearest_junction (provides node coordinates).")
    ap.add_argument("--include-outfalls-as-candidates", action="store_true",
                    help="When mode=nearest_junction, include outfall nodes alongside junctions.")
    ap.add_argument("--candidates-geojson", default=None,
                    help="Required for mode=nearest_catch_basin (Point or polygon layer).")
    ap.add_argument("--candidates-id-field", default="node_id",
                    help="Field name in --candidates-geojson that holds the node id.")
    ap.add_argument("--lookup-csv", default=None,
                    help="Required for mode=manual_lookup; columns: subcatchment_id,outlet_node_id.")
    return ap.parse_args()


def main() -> None:
    args = parse_args()
    csv_in = Path(args.subcatchments_csv_in)
    geojson_in = Path(args.subcatchments_geojson)
    csv_out = Path(args.out_csv)
    for p in (csv_in, geojson_in):
        if not p.exists():
            raise FileNotFoundError(p)

    centroids = _load_subcatchment_centroids(geojson_in)

    assignments: dict[str, str] = {}
    distances: dict[str, float] = {}

    if args.mode == "nearest_junction":
        if not args.network_json:
            raise ValueError("mode=nearest_junction requires --network-json")
        candidates = _candidate_nodes_from_network(Path(args.network_json), include_outfalls=args.include_outfalls_as_candidates)
        for sid, c in centroids:
            nid, d = _nearest(c, candidates)
            assignments[sid] = nid
            distances[sid] = d
    elif args.mode == "nearest_catch_basin":
        if not args.candidates_geojson:
            raise ValueError("mode=nearest_catch_basin requires --candidates-geojson")
        candidates = _candidate_nodes_from_geojson(Path(args.candidates_geojson), args.candidates_id_field)
        for sid, c in centroids:
            nid, d = _nearest(c, candidates)
            assignments[sid] = nid
            distances[sid] = d
    elif args.mode == "manual_lookup":
        if not args.lookup_csv:
            raise ValueError("mode=manual_lookup requires --lookup-csv")
        lookup = _read_manual_lookup(Path(args.lookup_csv))
        for sid, _ in centroids:
            if sid not in lookup:
                raise ValueError(f"subcatchment {sid} missing from lookup CSV")
            assignments[sid] = lookup[sid]
    else:
        raise ValueError(f"unknown mode: {args.mode}")

    # Rewrite the CSV with the new outlet column.
    rows_in: list[dict[str, str]] = []
    with csv_in.open("r", encoding="utf-8", newline="") as f:
        reader = csv.DictReader(f)
        if not reader.fieldnames or "subcatchment_id" not in reader.fieldnames or "outlet" not in reader.fieldnames:
            raise ValueError(
                f"{csv_in} must have 'subcatchment_id' and 'outlet' columns; got {reader.fieldnames}"
            )
        fieldnames = list(reader.fieldnames)
        for row in reader:
            sid = row["subcatchment_id"]
            if sid in assignments:
                row["outlet"] = assignments[sid]
            rows_in.append(row)

    csv_out.parent.mkdir(parents=True, exist_ok=True)
    with csv_out.open("w", encoding="utf-8", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=fieldnames)
        writer.writeheader()
        writer.writerows(rows_in)

    report = {
        "ok": True,
        "skill": "swmm-network",
        "tool": "assign_subcatchment_outlets",
        "mode": args.mode,
        "counts": {
            "subcatchments_assigned": len(assignments),
            "subcatchments_in_csv": len(rows_in),
        },
        "assignments": [
            {
                "subcatchment_id": sid,
                "outlet_node_id": nid,
                "centroid_to_node_distance_m": distances.get(sid),
            }
            for sid, nid in assignments.items()
        ],
        "outputs": {"subcatchments_csv": str(csv_out)},
        "inputs": {
            "subcatchments_csv_in": str(csv_in),
            "subcatchments_geojson": str(geojson_in),
            "network_json": str(args.network_json) if args.network_json else None,
            "candidates_geojson": str(args.candidates_geojson) if args.candidates_geojson else None,
            "lookup_csv": str(args.lookup_csv) if args.lookup_csv else None,
        },
    }
    print(json.dumps(report, indent=2))


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

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

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

from _hash_util import sha256_file


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


def save_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 read_csv_rows(path: Path | None) -> list[dict[str, str]]:
    if path is None:
        return []
    with path.open("r", encoding="utf-8", newline="") as f:
        return list(csv.DictReader(f))


def read_geojson_features(path: Path | None) -> list[dict[str, Any]]:
    if path is None:
        return []
    obj = load_json(path)
    if obj.get("type") != "FeatureCollection":
        raise ValueError(f"Expected GeoJSON FeatureCollection: {path}")
    return list(obj.get("features") or [])


def as_float(value: Any, default: float | None = None) -> float | None:
    if value is None or str(value).strip() == "":
        return default
    return float(str(value).strip())


def as_int(value: Any, default: int = 1) -> int:
    if value is None or str(value).strip() == "":
        return default
    return int(float(str(value).strip()))


def as_bool(value: Any, default: bool = False) -> bool:
    if value is None or str(value).strip() == "":
        return default
    if isinstance(value, bool):
        return value
    return str(value).strip().lower() in {"1", "true", "yes", "y"}


def text(value: Any, default: str = "") -> str:
    if value is None:
        return default
    return str(value).strip()


def field_value(props: dict[str, Any], fields: dict[str, str | None], key: str, default: Any = None) -> Any:
    field = fields.get(key)
    if not field:
        return default
    return props.get(field, default)


def coord_key(x: float, y: float, precision: int) -> str:
    return f"{round(x, precision):.{precision}f},{round(y, precision):.{precision}f}"


def default_node_id(prefix: str, x: float, y: float, precision: int) -> str:
    token = coord_key(x, y, precision).replace("-", "m").replace(".", "p").replace(",", "_")
    return f"{prefix}_{token}"


def line_length(coords: list[list[float]]) -> float:
    total = 0.0
    for a, b in zip(coords, coords[1:]):
        total += math.hypot(float(b[0]) - float(a[0]), float(b[1]) - float(a[1]))
    return total


def csv_pipe_geometry(row: dict[str, Any], fields: dict[str, str | None]) -> list[list[float]]:
    required = ["from_x", "from_y", "to_x", "to_y"]
    if not all(fields.get(k) and text(row.get(fields[k])) for k in required):
        return []
    return [
        [float(row[str(fields["from_x"])]), float(row[str(fields["from_y"])])],
        [float(row[str(fields["to_x"])]), float(row[str(fields["to_y"])])],
    ]


def geojson_pipe_geometry(feature: dict[str, Any]) -> list[list[float]]:
    geom = feature.get("geometry") or {}
    if geom.get("type") != "LineString":
        raise ValueError("Pipe GeoJSON features must use LineString geometry")
    coords = geom.get("coordinates") or []
    if len(coords) < 2:
        raise ValueError("Pipe LineString must contain at least two coordinates")
    return [[float(xy[0]), float(xy[1])] for xy in coords]


def point_geometry(feature: dict[str, Any]) -> tuple[float, float]:
    geom = feature.get("geometry") or {}
    if geom.get("type") != "Point":
        raise ValueError("Node/outfall GeoJSON features must use Point geometry")
    coords = geom.get("coordinates") or []
    return float(coords[0]), float(coords[1])


def csv_point_geometry(row: dict[str, Any], fields: dict[str, str | None]) -> tuple[float, float]:
    x_field = fields.get("x")
    y_field = fields.get("y")
    if not x_field or not y_field:
        raise ValueError("CSV node/outfall mapping requires x and y fields")
    return float(row[x_field]), float(row[y_field])


def load_point_assets(
    *,
    csv_path: Path | None,
    geojson_path: Path | None,
    cfg: dict[str, Any],
    asset_kind: str,
) -> list[dict[str, Any]]:
    fields = cfg.get("fields", {})
    defaults = cfg.get("defaults", {})
    rows = [{"props": row, "xy": csv_point_geometry(row, fields), "source": str(csv_path)} for row in read_csv_rows(csv_path)]
    rows.extend(
        {"props": feat.get("properties") or {}, "xy": point_geometry(feat), "source": str(geojson_path)}
        for feat in read_geojson_features(geojson_path)
    )

    out = []
    for idx, rec in enumerate(rows, start=1):
        props = rec["props"]
        x, y = rec["xy"]
        node_id = text(field_value(props, fields, "id"), f"{asset_kind.upper()}_{idx}")
        common = {
            "id": node_id,
            "invert_elev": as_float(field_value(props, fields, "invert_elev"), defaults.get("invert_elev", 0.0)),
            "coordinates": {"x": x, "y": y},
            "source_asset": {
                "kind": asset_kind,
                "source": rec["source"],
                "source_id": text(field_value(props, fields, "source_id"), node_id),
            },
        }
        if asset_kind == "outfall":
            out.append(
                {
                    **common,
                    "type": text(field_value(props, fields, "type"), defaults.get("type", "FREE")).upper(),
                    "stage_data": field_value(props, fields, "stage_data", defaults.get("stage_data")),
                    "gated": as_bool(field_value(props, fields, "gated"), defaults.get("gated", False)),
                    "route_to": field_value(props, fields, "route_to", defaults.get("route_to")),
                    "asset_type": text(field_value(props, fields, "asset_type"), defaults.get("asset_type", "storm")),
                    "system_layer": text(field_value(props, fields, "system_layer"), defaults.get("system_layer", "minor_pipe")),
                }
            )
        else:
            out.append(
                {
                    **common,
                    "max_depth": as_float(field_value(props, fields, "max_depth"), defaults.get("max_depth", 2.0)),
                    "init_depth": as_float(field_value(props, fields, "init_depth"), defaults.get("init_depth", 0.0)),
                    "sur_depth": as_float(field_value(props, fields, "sur_depth"), defaults.get("sur_depth", 0.0)),
                    "aponded": as_float(field_value(props, fields, "aponded"), defaults.get("aponded", 0.0)),
                    "asset_type": text(field_value(props, fields, "asset_type"), defaults.get("asset_type", "storm")),
                    "system_layer": text(field_value(props, fields, "system_layer"), defaults.get("system_layer", "minor_pipe")),
                }
            )
    return out


def collect_pipe_records(csv_path: Path | None, geojson_path: Path | None, fields: dict[str, str | None]) -> list[dict[str, Any]]:
    records = []
    for row in read_csv_rows(csv_path):
        records.append({"props": row, "coords": csv_pipe_geometry(row, fields), "source": str(csv_path)})
    for feat in read_geojson_features(geojson_path):
        records.append({"props": feat.get("properties") or {}, "coords": geojson_pipe_geometry(feat), "source": str(geojson_path)})
    return records


def add_inferred_node(
    *,
    nodes_by_id: dict[str, dict[str, Any]],
    coord_to_id: dict[str, str],
    node_id: str,
    x: float,
    y: float,
    invert_elev: float,
    max_depth: float,
    precision: int,
    source_pipe: str,
    system_layer: str,
    asset_type: str,
) -> None:
    key = coord_key(x, y, precision)
    if node_id in nodes_by_id:
        coord_to_id.setdefault(key, node_id)
        return
    nodes_by_id[node_id] = {
        "id": node_id,
        "invert_elev": invert_elev,
        "max_depth": max_depth,
        "init_depth": 0.0,
        "sur_depth": 0.0,
        "aponded": 0.0,
        "coordinates": {"x": x, "y": y},
        "asset_type": asset_type,
        "system_layer": system_layer,
        "inferred": True,
        "source_asset": {
            "kind": "inferred_junction",
            "source_pipe": source_pipe,
            "coordinate_key": key,
        },
    }
    coord_to_id[key] = node_id


def build_network(args: argparse.Namespace, mapping: dict[str, Any]) -> dict[str, Any]:
    precision = int(mapping.get("coordinate_precision", 3))
    pipes_cfg = mapping.get("pipes", {})
    node_cfg = mapping.get("junctions", {})
    outfall_cfg = mapping.get("outfalls", {})
    pipe_fields = pipes_cfg.get("fields", {})
    pipe_defaults = pipes_cfg.get("defaults", {})
    infer_cfg = mapping.get("inference", {})

    explicit_junctions = load_point_assets(
        csv_path=args.junctions_csv,
        geojson_path=args.junctions_geojson,
        cfg=node_cfg,
        asset_kind="junction",
    )
    explicit_outfalls = load_point_assets(
        csv_path=args.outfalls_csv,
        geojson_path=args.outfalls_geojson,
        cfg=outfall_cfg,
        asset_kind="outfall",
    )

    nodes_by_id = {str(j["id"]): j for j in explicit_junctions}
    outfalls_by_id = {str(o["id"]): o for o in explicit_outfalls}
    coord_to_id: dict[str, str] = {}
    for rec in explicit_junctions + explicit_outfalls:
        xy = rec["coordinates"]
        coord_to_id[coord_key(float(xy["x"]), float(xy["y"]), precision)] = str(rec["id"])

    pipe_records = collect_pipe_records(args.pipes_csv, args.pipes_geojson, pipe_fields)
    conduits = []
    inferred_count = 0
    for idx, rec in enumerate(pipe_records, start=1):
        props = rec["props"]
        coords = rec["coords"]
        if not coords:
            raise ValueError(f"Pipe record {idx} has no geometry or endpoint coordinate fields")
        start = coords[0]
        end = coords[-1]
        start_key = coord_key(float(start[0]), float(start[1]), precision)
        end_key = coord_key(float(end[0]), float(end[1]), precision)
        pipe_id = text(field_value(props, pipe_fields, "id"), f"P{idx}")
        system_layer = text(field_value(props, pipe_fields, "system_layer"), pipe_defaults.get("system_layer", "minor_pipe"))
        asset_type = text(field_value(props, pipe_fields, "asset_type"), pipe_defaults.get("asset_type", "storm"))
        from_node = text(field_value(props, pipe_fields, "from_node"), coord_to_id.get(start_key, ""))
        to_node = text(field_value(props, pipe_fields, "to_node"), coord_to_id.get(end_key, ""))

        if not from_node:
            from_node = default_node_id(str(infer_cfg.get("junction_prefix", "J")), float(start[0]), float(start[1]), precision)
        if not to_node:
            to_node = default_node_id(str(infer_cfg.get("junction_prefix", "J")), float(end[0]), float(end[1]), precision)

        from_invert = as_float(field_value(props, pipe_fields, "from_invert_elev"), pipe_defaults.get("from_invert_elev", 0.0))
        to_invert = as_float(field_value(props, pipe_fields, "to_invert_elev"), pipe_defaults.get("to_invert_elev", from_invert))
        max_depth = float(infer_cfg.get("max_depth", 2.0))
        if from_node not in outfalls_by_id and from_node not in nodes_by_id:
            inferred_count += 1
            add_inferred_node(
                nodes_by_id=nodes_by_id,
                coord_to_id=coord_to_id,
                node_id=from_node,
                x=float(start[0]),
                y=float(start[1]),
                invert_elev=float(from_invert or 0.0),
                max_depth=max_depth,
                precision=precision,
                source_pipe=pipe_id,
                system_layer=system_layer,
                asset_type=asset_type,
            )
        if to_node not in outfalls_by_id and to_node not in nodes_by_id:
            inferred_count += 1
            add_inferred_node(
                nodes_by_id=nodes_by_id,
                coord_to_id=coord_to_id,
                node_id=to_node,
                x=float(end[0]),
                y=float(end[1]),
                invert_elev=float(to_invert or 0.0),
                max_depth=max_depth,
                precision=precision,
                source_pipe=pipe_id,
                system_layer=system_layer,
                asset_type=asset_type,
            )

        vertices = [{"x": float(x), "y": float(y)} for x, y in coords[1:-1]]
        length = as_float(field_value(props, pipe_fields, "length"), None)
        if length is None:
            length = max(line_length(coords), float(pipe_defaults.get("minimum_length", 1.0)))
        geom1 = as_float(
            field_value(props, pipe_fields, "geom1", field_value(props, pipe_fields, "diameter")),
            pipe_defaults.get("geom1", 0.5),
        )
        conduits.append(
            {
                "id": pipe_id,
                "from_node": from_node,
                "to_node": to_node,
                "length": length,
                "roughness": as_float(field_value(props, pipe_fields, "roughness"), pipe_defaults.get("roughness", 0.013)),
                "in_offset": as_float(field_value(props, pipe_fields, "in_offset"), pipe_defaults.get("in_offset", 0.0)),
                "out_offset": as_float(field_value(props, pipe_fields, "out_offset"), pipe_defaults.get("out_offset", 0.0)),
                "init_flow": as_float(field_value(props, pipe_fields, "init_flow"), pipe_defaults.get("init_flow", 0.0)),
                "max_flow": as_float(field_value(props, pipe_fields, "max_flow"), pipe_defaults.get("max_flow")),
                "xsection": {
                    "shape": text(field_value(props, pipe_fields, "shape"), pipe_defaults.get("shape", "CIRCULAR")).upper(),
                    "geom1": geom1,
                    "geom2": as_float(field_value(props, pipe_fields, "geom2"), pipe_defaults.get("geom2", 0.0)),
                    "geom3": as_float(field_value(props, pipe_fields, "geom3"), pipe_defaults.get("geom3", 0.0)),
                    "geom4": as_float(field_value(props, pipe_fields, "geom4"), pipe_defaults.get("geom4", 0.0)),
                    "barrels": as_int(field_value(props, pipe_fields, "barrels"), int(pipe_defaults.get("barrels", 1))),
                },
                "vertices": vertices,
                "asset_type": asset_type,
                "system_layer": system_layer,
                "material": text(field_value(props, pipe_fields, "material"), pipe_defaults.get("material", "")),
                "source_asset": {
                    "kind": "pipe",
                    "source": rec["source"],
                    "source_id": text(field_value(props, pipe_fields, "source_id"), pipe_id),
                },
            }
        )

    layers = sorted({c.get("system_layer", "") for c in conduits if c.get("system_layer")})
    network = {
        "meta": {
            **(mapping.get("meta") or {}),
            "adapter": "city_network_adapter",
            "dual_system_ready": bool(mapping.get("dual_system_ready", True)),
            "system_layers": layers,
            "sources": {
                "pipes_csv": str(args.pipes_csv) if args.pipes_csv else None,
                "pipes_geojson": str(args.pipes_geojson) if args.pipes_geojson else None,
                "junctions_csv": str(args.junctions_csv) if args.junctions_csv else None,
                "junctions_geojson": str(args.junctions_geojson) if args.junctions_geojson else None,
                "outfalls_csv": str(args.outfalls_csv) if args.outfalls_csv else None,
                "outfalls_geojson": str(args.outfalls_geojson) if args.outfalls_geojson else None,
            },
            "source_hashes": {
                "pipes_csv": sha256_file(args.pipes_csv),
                "pipes_geojson": sha256_file(args.pipes_geojson),
                "junctions_csv": sha256_file(args.junctions_csv),
                "junctions_geojson": sha256_file(args.junctions_geojson),
                "outfalls_csv": sha256_file(args.outfalls_csv),
                "outfalls_geojson": sha256_file(args.outfalls_geojson),
            },
            "counts": {
                "pipes": len(conduits),
                "explicit_junctions": len(explicit_junctions),
                "explicit_outfalls": len(explicit_outfalls),
                "inferred_junctions": inferred_count,
            },
        },
        "junctions": sorted(nodes_by_id.values(), key=lambda x: str(x["id"])),
        "outfalls": sorted(outfalls_by_id.values(), key=lambda x: str(x["id"])),
        "conduits": conduits,
    }
    return network


def main() -> None:
    ap = argparse.ArgumentParser(
        description="Convert structured urban pipe/surface asset exports into Agentic SWMM network.json."
    )
    ap.add_argument("--pipes-csv", type=Path, default=None)
    ap.add_argument("--pipes-geojson", type=Path, default=None)
    ap.add_argument("--junctions-csv", type=Path, default=None)
    ap.add_argument("--junctions-geojson", type=Path, default=None)
    ap.add_argument("--outfalls-csv", type=Path, default=None)
    ap.add_argument("--outfalls-geojson", type=Path, default=None)
    ap.add_argument("--mapping-json", type=Path, required=True)
    ap.add_argument("--out", type=Path, required=True)
    args = ap.parse_args()

    if args.pipes_csv is None and args.pipes_geojson is None:
        raise ValueError("At least one pipe source is required: --pipes-csv or --pipes-geojson")

    mapping = load_json(args.mapping_json)
    network = build_network(args, mapping)
    save_json(args.out, network)
    print(
        json.dumps(
            {
                "ok": True,
                "out": str(args.out),
                "junction_count": len(network["junctions"]),
                "outfall_count": len(network["outfalls"]),
                "conduit_count": len(network["conduits"]),
                "system_layers": network["meta"]["system_layers"],
                "inferred_junctions": network["meta"]["counts"]["inferred_junctions"],
            },
            indent=2,
        )
    )


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Pick a SWMM outfall point from raw pipe and watercourse layers.

Two pluggable modes:

- ``endpoint_nearest_watercourse`` (default): scan every pipe endpoint,
  measure its distance to the nearest watercourse geometry, and pick the
  endpoint with the smallest distance. Best when a digitised watercourse
  layer covers the basin's receiving water.

- ``lowest_endpoint``: pick the pipe endpoint with the smallest y
  coordinate. Useful as a fallback when no watercourse layer is
  available; assumes a projected, north-positive CRS so that the
  topographic low is south.

The tool always emits a single outfall (``node_id = OUT1``). Multi-outfall
networks need a follow-up tool; this one is intentionally minimal.
"""
from __future__ import annotations

import argparse
import json
import sys
from pathlib import Path

from shapely.geometry import Point, shape


MODES = ("endpoint_nearest_watercourse", "lowest_endpoint")


def _collect_endpoints(pipes_geojson: dict) -> list[dict]:
    """Return list of {pipe_index, pipe_id, position, point} entries."""
    out: list[dict] = []
    for idx, feat in enumerate(pipes_geojson.get("features") or []):
        geom = feat.get("geometry") or {}
        if geom.get("type") != "LineString":
            continue
        coords = geom.get("coordinates") or []
        if len(coords) < 2:
            continue
        pid = (
            (feat.get("properties") or {}).get("FACILITYID")
            or (feat.get("properties") or {}).get("id")
            or f"pipe_{idx}"
        )
        out.append({
            "pipe_index": idx,
            "pipe_id": pid,
            "position": "start",
            "point": Point(float(coords[0][0]), float(coords[0][1])),
        })
        out.append({
            "pipe_index": idx,
            "pipe_id": pid,
            "position": "end",
            "point": Point(float(coords[-1][0]), float(coords[-1][1])),
        })
    if not out:
        raise ValueError("no LineString features with endpoints in pipes geojson")
    return out


def _watercourse_geoms(watercourse_geojson: dict) -> list:
    geoms = []
    for feat in watercourse_geojson.get("features") or []:
        geom = feat.get("geometry")
        if geom is None:
            continue
        try:
            geoms.append(shape(geom))
        except Exception:
            continue
    if not geoms:
        raise ValueError("watercourse geojson has no usable geometries")
    return geoms


def _pick_endpoint_nearest_watercourse(endpoints: list[dict], watercourse: list) -> dict:
    best = None
    best_dist = float("inf")
    for ep in endpoints:
        dist = min(ep["point"].distance(g) for g in watercourse)
        if dist < best_dist:
            best_dist = dist
            best = {**ep, "distance_to_watercourse_m": float(dist)}
    if best is None:
        raise RuntimeError("no endpoint chosen — pipes empty?")
    return best


def _pick_lowest_endpoint(endpoints: list[dict]) -> dict:
    best = min(endpoints, key=lambda ep: ep["point"].y)
    return {**best, "distance_to_watercourse_m": None}


def build_outfalls_geojson(chosen: dict, source_crs: dict | None) -> dict:
    pt = chosen["point"]
    properties = {
        "node_id": "OUT1",
        "type": "FREE",
        "invert_elev": 0.0,
        "asset_type": "storm",
        "system_layer": "minor_pipe",
        "source_pipe": chosen["pipe_id"],
        "source_position": chosen["position"],
    }
    if chosen.get("distance_to_watercourse_m") is not None:
        properties["dist_to_watercourse_m"] = chosen["distance_to_watercourse_m"]
    out = {
        "type": "FeatureCollection",
        "name": "outfalls",
        "features": [
            {
                "type": "Feature",
                "properties": properties,
                "geometry": {"type": "Point", "coordinates": [pt.x, pt.y]},
            }
        ],
    }
    if source_crs is not None:
        out["crs"] = source_crs
    return out


def parse_args() -> argparse.Namespace:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--pipes-geojson", required=True)
    ap.add_argument(
        "--watercourse-geojson",
        default=None,
        help="Required for mode=endpoint_nearest_watercourse; ignored for mode=lowest_endpoint.",
    )
    ap.add_argument("--mode", choices=MODES, default="endpoint_nearest_watercourse")
    ap.add_argument("--out", required=True)
    return ap.parse_args()


def main() -> None:
    args = parse_args()
    pipes_path = Path(args.pipes_geojson)
    out_path = Path(args.out)
    if not pipes_path.exists():
        raise FileNotFoundError(pipes_path)

    pipes_geojson = json.loads(pipes_path.read_text(encoding="utf-8"))
    endpoints = _collect_endpoints(pipes_geojson)

    if args.mode == "endpoint_nearest_watercourse":
        if not args.watercourse_geojson:
            raise ValueError("--watercourse-geojson required for mode=endpoint_nearest_watercourse")
        wc_path = Path(args.watercourse_geojson)
        if not wc_path.exists():
            raise FileNotFoundError(wc_path)
        wc_geojson = json.loads(wc_path.read_text(encoding="utf-8"))
        watercourse = _watercourse_geoms(wc_geojson)
        chosen = _pick_endpoint_nearest_watercourse(endpoints, watercourse)
    elif args.mode == "lowest_endpoint":
        chosen = _pick_lowest_endpoint(endpoints)
    else:
        raise ValueError(f"unknown mode: {args.mode}")

    source_crs = pipes_geojson.get("crs")
    outfalls = build_outfalls_geojson(chosen, source_crs)
    out_path.parent.mkdir(parents=True, exist_ok=True)
    out_path.write_text(json.dumps(outfalls, indent=2), encoding="utf-8")

    report = {
        "ok": True,
        "skill": "swmm-network",
        "tool": "infer_outfall",
        "mode": args.mode,
        "chosen": {
            "pipe_id": chosen["pipe_id"],
            "position": chosen["position"],
            "x": chosen["point"].x,
            "y": chosen["point"].y,
            "distance_to_watercourse_m": chosen.get("distance_to_watercourse_m"),
        },
        "counts": {
            "endpoints_considered": len(endpoints),
            "outfalls_emitted": 1,
        },
        "outputs": {"outfalls_geojson": str(out_path)},
        "inputs": {
            "pipes_geojson": str(pipes_path),
            "watercourse_geojson": str(args.watercourse_geojson) if args.watercourse_geojson else None,
        },
    }
    print(json.dumps(report, indent=2))


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

#!/usr/bin/env python3
"""Convert a SWMM ``.inp`` into the swmm-network ``network.json`` shape.

This is the input bridge that lets ``network_qa.py`` run its structural
checks (isolated_node / no_outfall_path / counts / xsection validity) on
ANY model that emitted an ``.inp`` — including a SWMManywhere-synthesized
network, which produces an ``.inp`` but no ``network.json``. With this,
the same structural QA serves the real-data paths (which already build a
``network.json``) and the synth path uniformly.

Standalone: stdlib only, no ``agentic_swmm`` / sibling-skill imports, so it
runs identically whether invoked by the MCP server, the CLI, or a test.

Only the sections the QA needs are parsed: ``[COORDINATES]``, ``[JUNCTIONS]``,
``[OUTFALLS]``, ``[CONDUITS]``, ``[XSECTIONS]``. Other sections are ignored.
A conduit with no matching ``[XSECTIONS]`` row is emitted without an
``xsection`` key on purpose — that is exactly what network_qa flags as
``missing_xsection``, so the gap stays visible rather than being papered over.
"""
from __future__ import annotations

import re
from pathlib import Path
from typing import Any

_SECTION_RE = re.compile(r"^\s*\[([A-Z_]+)\]\s*$")


def parse_inp_sections(text: str) -> dict[str, list[str]]:
    """Split INP text into ``{SECTION: [non-comment data rows]}``."""
    sections: dict[str, list[str]] = {}
    current: str | None = None
    for raw in text.splitlines():
        stripped = raw.strip()
        if not stripped or stripped.startswith(";"):
            continue
        m = _SECTION_RE.match(raw)
        if m:
            current = m.group(1).upper()
            sections.setdefault(current, [])
            continue
        if current is None:
            continue
        sections[current].append(stripped)
    return sections


def inp_to_network(inp_path: str | Path) -> dict[str, Any]:
    """Build the ``network.json`` dict that ``network_qa.run_qa`` consumes."""
    text = Path(inp_path).read_text(encoding="utf-8", errors="replace")
    sections = parse_inp_sections(text)

    # [COORDINATES]: Node Xcoord Ycoord
    coords: dict[str, dict[str, float]] = {}
    for row in sections.get("COORDINATES", []):
        cols = row.split()
        if len(cols) >= 3:
            try:
                coords[cols[0]] = {"x": float(cols[1]), "y": float(cols[2])}
            except ValueError:
                continue

    def _node(name: str) -> dict[str, Any]:
        node: dict[str, Any] = {"id": name}
        if name in coords:
            node["coordinates"] = coords[name]
        return node

    junctions = [_node(r.split()[0]) for r in sections.get("JUNCTIONS", []) if r.split()]
    outfalls = [_node(r.split()[0]) for r in sections.get("OUTFALLS", []) if r.split()]

    # [XSECTIONS]: Link Shape Geom1 ...
    xsections: dict[str, dict[str, Any]] = {}
    for row in sections.get("XSECTIONS", []):
        cols = row.split()
        if len(cols) >= 3:
            try:
                xsections[cols[0]] = {"shape": cols[1].upper(), "geom1": float(cols[2])}
            except ValueError:
                continue

    # [CONDUITS]: Name FromNode ToNode Length Roughness ...
    conduits: list[dict[str, Any]] = []
    for row in sections.get("CONDUITS", []):
        cols = row.split()
        if len(cols) < 5:
            continue
        try:
            conduit: dict[str, Any] = {
                "id": cols[0],
                "from_node": cols[1],
                "to_node": cols[2],
                "length": float(cols[3]),
                "roughness": float(cols[4]),
            }
        except ValueError:
            continue
        if cols[0] in xsections:
            conduit["xsection"] = xsections[cols[0]]
        conduits.append(conduit)

    return {
        "junctions": junctions,
        "outfalls": outfalls,
        "conduits": conduits,
        "meta": {"source": "inp", "inp_path": str(inp_path)},
    }

#!/usr/bin/env python3
"""Emit a minimal 1-junction + 1-outfall stub network.json.

Use when the basin has a known outlet point + name but no pipe-network
geometry yet. Drives the swmm-builder builder MCP for smoke runs and
real-data fallbacks without inventing a multi-pipe drainage system.
"""
from __future__ import annotations

import argparse
import json
from pathlib import Path

import shapefile


def _resolve_field_index(reader: shapefile.Reader, name: str) -> int:
    fields = [field[0] for field in reader.fields[1:]]
    if name not in fields:
        raise SystemExit(
            f"Field '{name}' not found in {reader.shp.name}. Available: {fields}"
        )
    return fields.index(name)


def main() -> None:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--shp", type=Path, required=True, help="Subcatchment shapefile carrying outlet attrs")
    ap.add_argument("--outlet-id-field", required=True)
    ap.add_argument("--outlet-x-field", required=True)
    ap.add_argument("--outlet-y-field", required=True)
    ap.add_argument("--invert-elev", type=float, required=True)
    ap.add_argument("--outfall-invert-elev", type=float, required=True)
    ap.add_argument("--outfall-offset-m", type=float, default=500.0)
    ap.add_argument("--conduit-roughness", type=float, default=0.013)
    ap.add_argument("--conduit-diameter", type=float, default=1.5)
    ap.add_argument("--outfall-id", default=None)
    ap.add_argument("--conduit-id", default="C_OUT")
    ap.add_argument("--max-depth", type=float, default=3.0)
    ap.add_argument("--flow-units", default="CMS")
    ap.add_argument("--out-json", type=Path, required=True)
    args = ap.parse_args()

    reader = shapefile.Reader(str(args.shp))
    idx_id = _resolve_field_index(reader, args.outlet_id_field)
    idx_x = _resolve_field_index(reader, args.outlet_x_field)
    idx_y = _resolve_field_index(reader, args.outlet_y_field)

    records = list(reader.records())
    if not records:
        raise SystemExit(f"No records in {args.shp}")
    first = records[0]
    outlet_id = str(first[idx_id]).strip()
    if not outlet_id:
        raise SystemExit(f"Blank outlet id in field '{args.outlet_id_field}' of {args.shp}")
    outlet_x = float(first[idx_x])
    outlet_y = float(first[idx_y])

    outfall_id = args.outfall_id or f"{outlet_id}_OUT"
    if outfall_id == outlet_id:
        raise SystemExit(f"--outfall-id must differ from junction id '{outlet_id}'")

    network = {
        "meta": {
            "name": f"{outlet_id}-minimal-stub",
            "flow_units": args.flow_units,
            "note": "Stub network produced by swmm-network/minimal_stub_network.py — single junction + free outfall.",
        },
        "junctions": [
            {
                "id": outlet_id,
                "invert_elev": args.invert_elev,
                "max_depth": args.max_depth,
                "init_depth": 0.0,
                "sur_depth": 0.0,
                "aponded": 0.0,
                "coordinates": {"x": outlet_x, "y": outlet_y},
            }
        ],
        "outfalls": [
            {
                "id": outfall_id,
                "invert_elev": args.outfall_invert_elev,
                "type": "FREE",
                "gated": False,
                "route_to": None,
                "coordinates": {"x": outlet_x + args.outfall_offset_m, "y": outlet_y},
            }
        ],
        "conduits": [
            {
                "id": args.conduit_id,
                "from_node": outlet_id,
                "to_node": outfall_id,
                "length": args.outfall_offset_m,
                "roughness": args.conduit_roughness,
                "in_offset": 0.0,
                "out_offset": 0.0,
                "init_flow": 0.0,
                "max_flow": None,
                "xsection": {
                    "shape": "CIRCULAR",
                    "geom1": args.conduit_diameter,
                    "geom2": 0.0,
                    "geom3": 0.0,
                    "geom4": 0.0,
                    "barrels": 1,
                },
            }
        ],
    }

    args.out_json.parent.mkdir(parents=True, exist_ok=True)
    args.out_json.write_text(json.dumps(network, indent=2), encoding="utf-8")
    print(json.dumps({"ok": True, "out_json": str(args.out_json), "junction": outlet_id, "outfall": outfall_id}, indent=2))


if __name__ == "__main__":
    main()

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

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


def load_json(path: str | Path) -> Any:
    return json.loads(Path(path).read_text())


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


def _as_float(value: Any, default: float | None = None) -> float | None:
    if value in (None, "", "null"):
        return default
    return float(value)


def _as_bool(value: Any, default: bool = False) -> bool:
    if value in (None, ""):
        return default
    if isinstance(value, bool):
        return value
    s = str(value).strip().lower()
    return s in {"1", "true", "yes", "y"}


def read_geojson_features(path: Path) -> list[dict]:
    obj = load_json(path)
    if obj.get("type") != "FeatureCollection":
        raise ValueError(f"Expected GeoJSON FeatureCollection: {path}")
    return obj.get("features", [])


def read_csv_rows(path: Path) -> list[dict]:
    with path.open(newline="", encoding="utf-8") as f:
        return list(csv.DictReader(f))


def get_value(props: dict, field: str | None, default: Any = None) -> Any:
    if not field:
        return default
    return props.get(field, default)


def import_junctions(path: Path, cfg: dict) -> list[dict]:
    fmt = cfg.get("format", "geojson")
    fields = cfg.get("fields", {})
    out = []
    if fmt == "geojson":
        for feat in read_geojson_features(path):
            props = feat.get("properties", {})
            geom = feat.get("geometry", {})
            coords = geom.get("coordinates", [None, None])
            out.append({
                "id": str(get_value(props, fields.get("id"))),
                "invert_elev": _as_float(get_value(props, fields.get("invert_elev")), 0.0),
                "max_depth": _as_float(get_value(props, fields.get("max_depth")), 0.0),
                "init_depth": _as_float(get_value(props, fields.get("init_depth")), 0.0),
                "sur_depth": _as_float(get_value(props, fields.get("sur_depth")), 0.0),
                "aponded": _as_float(get_value(props, fields.get("aponded")), 0.0),
                "coordinates": {"x": float(coords[0]), "y": float(coords[1])},
            })
    elif fmt == "csv":
        for row in read_csv_rows(path):
            out.append({
                "id": str(get_value(row, fields.get("id"))),
                "invert_elev": _as_float(get_value(row, fields.get("invert_elev")), 0.0),
                "max_depth": _as_float(get_value(row, fields.get("max_depth")), 0.0),
                "init_depth": _as_float(get_value(row, fields.get("init_depth")), 0.0),
                "sur_depth": _as_float(get_value(row, fields.get("sur_depth")), 0.0),
                "aponded": _as_float(get_value(row, fields.get("aponded")), 0.0),
                "coordinates": {
                    "x": float(get_value(row, fields.get("x"))),
                    "y": float(get_value(row, fields.get("y"))),
                },
            })
    else:
        raise ValueError(f"Unsupported junction format: {fmt}")
    return out


def import_outfalls(path: Path, cfg: dict) -> list[dict]:
    fmt = cfg.get("format", "geojson")
    fields = cfg.get("fields", {})
    out = []
    if fmt == "geojson":
        for feat in read_geojson_features(path):
            props = feat.get("properties", {})
            geom = feat.get("geometry", {})
            coords = geom.get("coordinates", [None, None])
            out.append({
                "id": str(get_value(props, fields.get("id"))),
                "invert_elev": _as_float(get_value(props, fields.get("invert_elev")), 0.0),
                "type": str(get_value(props, fields.get("type"), cfg.get("defaults", {}).get("type", "FREE"))),
                "stage_data": get_value(props, fields.get("stage_data")),
                "gated": _as_bool(get_value(props, fields.get("gated")), cfg.get("defaults", {}).get("gated", False)),
                "route_to": get_value(props, fields.get("route_to")),
                "coordinates": {"x": float(coords[0]), "y": float(coords[1])},
            })
    elif fmt == "csv":
        for row in read_csv_rows(path):
            out.append({
                "id": str(get_value(row, fields.get("id"))),
                "invert_elev": _as_float(get_value(row, fields.get("invert_elev")), 0.0),
                "type": str(get_value(row, fields.get("type"), cfg.get("defaults", {}).get("type", "FREE"))),
                "stage_data": get_value(row, fields.get("stage_data")),
                "gated": _as_bool(get_value(row, fields.get("gated")), cfg.get("defaults", {}).get("gated", False)),
                "route_to": get_value(row, fields.get("route_to")),
                "coordinates": {
                    "x": float(get_value(row, fields.get("x"))),
                    "y": float(get_value(row, fields.get("y"))),
                },
            })
    else:
        raise ValueError(f"Unsupported outfall format: {fmt}")
    return out


def import_conduits(path: Path, cfg: dict) -> list[dict]:
    fmt = cfg.get("format", "geojson")
    if fmt != "geojson":
        raise ValueError("MVP conduit import currently supports geojson only")
    fields = cfg.get("fields", {})
    defaults = cfg.get("defaults", {})
    out = []
    for feat in read_geojson_features(path):
        props = feat.get("properties", {})
        geom = feat.get("geometry", {})
        coords = geom.get("coordinates", [])
        vertices = []
        if geom.get("type") == "LineString" and len(coords) > 2:
            for xy in coords[1:-1]:
                vertices.append({"x": float(xy[0]), "y": float(xy[1])})
        out.append({
            "id": str(get_value(props, fields.get("id"))),
            "from_node": str(get_value(props, fields.get("from_node"))),
            "to_node": str(get_value(props, fields.get("to_node"))),
            "length": _as_float(get_value(props, fields.get("length")), defaults.get("length", 1.0)),
            "roughness": _as_float(get_value(props, fields.get("roughness")), defaults.get("roughness", 0.013)),
            "in_offset": _as_float(get_value(props, fields.get("in_offset")), defaults.get("in_offset", 0.0)),
            "out_offset": _as_float(get_value(props, fields.get("out_offset")), defaults.get("out_offset", 0.0)),
            "init_flow": _as_float(get_value(props, fields.get("init_flow")), defaults.get("init_flow", 0.0)),
            "max_flow": _as_float(get_value(props, fields.get("max_flow")), defaults.get("max_flow")),
            "xsection": {
                "shape": str(get_value(props, fields.get("shape"), defaults.get("shape", "CIRCULAR"))),
                "geom1": _as_float(get_value(props, fields.get("geom1") or fields.get("diameter")), defaults.get("geom1", 0.5)),
                "geom2": _as_float(get_value(props, fields.get("geom2")), defaults.get("geom2", 0.0)),
                "geom3": _as_float(get_value(props, fields.get("geom3")), defaults.get("geom3", 0.0)),
                "geom4": _as_float(get_value(props, fields.get("geom4")), defaults.get("geom4", 0.0)),
                "barrels": int(get_value(props, fields.get("barrels"), defaults.get("barrels", 1))),
            },
            **({"vertices": vertices} if vertices else {}),
        })
    return out


def main() -> None:
    ap = argparse.ArgumentParser()
    ap.add_argument("--conduits", type=Path, required=True)
    ap.add_argument("--junctions", type=Path, required=True)
    ap.add_argument("--outfalls", type=Path, required=True)
    ap.add_argument("--mapping", type=Path, required=True)
    ap.add_argument("--out", type=Path, required=True)
    args = ap.parse_args()

    mapping = load_json(args.mapping)
    network = {
        "meta": mapping.get("meta", {}),
        "junctions": import_junctions(args.junctions, mapping["junctions"]),
        "outfalls": import_outfalls(args.outfalls, mapping["outfalls"]),
        "conduits": import_conduits(args.conduits, mapping["conduits"]),
    }
    save_json(args.out, network)
    print(json.dumps({
        "ok": True,
        "out": str(args.out),
        "summary": {
            "junction_count": len(network["junctions"]),
            "outfall_count": len(network["outfalls"]),
            "conduit_count": len(network["conduits"]),
        }
    }, indent=2))


if __name__ == "__main__":
    main()

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

import argparse
import json
from collections import defaultdict, deque
from pathlib import Path
from typing import Any


def load_json(path: str | Path) -> Any:
    return json.loads(Path(path).read_text())


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


def add_issue(issues: list[dict], severity: str, code: str, message: str, obj_id: str | None = None) -> None:
    rec = {"severity": severity, "code": code, "message": message}
    if obj_id is not None:
        rec["id"] = obj_id
    issues.append(rec)


def summarize(network: dict) -> dict:
    conduits = network.get("conduits", [])
    junctions = network.get("junctions", [])
    outfalls = network.get("outfalls", [])
    system_layers = sorted({str(c.get("system_layer")) for c in conduits if c.get("system_layer")})
    asset_types = sorted({str(c.get("asset_type")) for c in conduits if c.get("asset_type")})
    return {
        "junction_count": len(junctions),
        "outfall_count": len(outfalls),
        "conduit_count": len(conduits),
        "total_conduit_length": float(sum(c.get("length", 0) for c in conduits)),
        "inferred_junction_count": sum(1 for j in junctions if j.get("inferred")),
        "system_layers": system_layers,
        "asset_types": asset_types,
        "dual_system_ready": bool((network.get("meta") or {}).get("dual_system_ready") or len(system_layers) > 1),
    }


def run_qa(network: dict) -> dict:
    """Run structural QA on a network dict, return the report dict.

    Pure: no I/O. ``main`` builds ``network`` from either a network.json or
    (via ``inp_to_network``) a SWMM ``.inp`` and hands it here, so every
    entrypoint shares one QA implementation.
    """
    issues: list[dict] = []

    node_ids = [j["id"] for j in network.get("junctions", [])] + [o["id"] for o in network.get("outfalls", [])]
    counts = defaultdict(int)
    for nid in node_ids:
        counts[nid] += 1
    for cid in [c["id"] for c in network.get("conduits", [])]:
        counts[cid] += 1
    for obj_id, count in counts.items():
        if count > 1:
            add_issue(issues, "error", "duplicate_id", f"Duplicate ID detected: {obj_id}", obj_id)

    all_nodes = {j["id"]: j for j in network.get("junctions", [])}
    all_nodes.update({o["id"]: o for o in network.get("outfalls", [])})

    for j in network.get("junctions", []):
        if "coordinates" not in j or "x" not in j["coordinates"] or "y" not in j["coordinates"]:
            add_issue(issues, "error", "missing_coordinates", "Junction missing valid coordinates", j["id"])
    for o in network.get("outfalls", []):
        if "coordinates" not in o or "x" not in o["coordinates"] or "y" not in o["coordinates"]:
            add_issue(issues, "error", "missing_coordinates", "Outfall missing valid coordinates", o["id"])

    incoming = defaultdict(int)
    outgoing = defaultdict(int)
    graph = defaultdict(list)
    outfalls = {o["id"] for o in network.get("outfalls", [])}

    for c in network.get("conduits", []):
        cid = c["id"]
        if c.get("from_node") not in all_nodes:
            add_issue(issues, "error", "missing_from_node", f"Conduit from_node not found: {c.get('from_node')}", cid)
        if c.get("to_node") not in all_nodes:
            add_issue(issues, "error", "missing_to_node", f"Conduit to_node not found: {c.get('to_node')}", cid)
        if c.get("length", 0) <= 0:
            add_issue(issues, "error", "non_positive_length", "Conduit length must be > 0", cid)
        if c.get("roughness", 0) <= 0:
            add_issue(issues, "error", "non_positive_roughness", "Conduit roughness must be > 0", cid)
        if not c.get("xsection"):
            add_issue(issues, "error", "missing_xsection", "Conduit missing xsection", cid)
        else:
            xs = c["xsection"]
            if xs.get("geom1", 0) <= 0:
                add_issue(issues, "error", "invalid_xsection", "Conduit xsection geom1 must be > 0", cid)
        fn = c.get("from_node")
        tn = c.get("to_node")
        if fn in all_nodes and tn in all_nodes:
            outgoing[fn] += 1
            incoming[tn] += 1
            graph[fn].append(tn)
        layer = c.get("system_layer")
        if layer and layer not in {"minor_pipe", "major_surface", "dual_link", "storm", "sanitary", "combined"}:
            add_issue(issues, "warning", "unknown_system_layer", f"Unrecognized system_layer: {layer}", cid)

    for j in network.get("junctions", []):
        jid = j["id"]
        if incoming[jid] == 0 and outgoing[jid] == 0:
            add_issue(issues, "warning", "isolated_node", "Junction is isolated", jid)

    # SWMM allows exactly one inlet link and zero outlet links per outfall
    # (engine ERROR 141; a system with no acceptable outlet then dies with
    # ERROR 145). This structural rule was missing, which let an unrunnable
    # network sail through QA with ok=true (found 2026-08-08 by mining the
    # city-dual-system benchmark's rpt against its manifest).
    for o in network.get("outfalls", []):
        oid = o["id"]
        if incoming[oid] > 1:
            add_issue(
                issues,
                "error",
                "outfall_multiple_inlets",
                f"Outfall has {incoming[oid]} inlet links; SWMM allows exactly one (engine ERROR 141)",
                oid,
            )
        if outgoing[oid] > 0:
            add_issue(
                issues,
                "error",
                "outfall_has_outlet",
                "Outfall has an outgoing link; SWMM allows none (engine ERROR 141)",
                oid,
            )

    for start in [j["id"] for j in network.get("junctions", [])]:
        q = deque([start])
        seen = {start}
        reaches_outfall = False
        while q:
            cur = q.popleft()
            if cur in outfalls:
                reaches_outfall = True
                break
            for nxt in graph[cur]:
                if nxt not in seen:
                    seen.add(nxt)
                    q.append(nxt)
        if not reaches_outfall:
            add_issue(issues, "warning", "no_outfall_path", "No downstream path from junction to any outfall", start)

    return {
        "ok": not any(x["severity"] == "error" for x in issues),
        "summary": summarize(network),
        "issue_count": len(issues),
        "issues": issues,
    }


def main() -> None:
    ap = argparse.ArgumentParser(
        description="Structural QA for a SWMM network — a network.json or a .inp.",
    )
    ap.add_argument("network_json", type=Path, nargs="?", default=None,
                    help="Path to a network.json (omit when using --inp).")
    ap.add_argument("--inp", type=Path, default=None,
                    help="Run QA on a SWMM .inp instead of a network.json "
                         "(e.g. a SWMManywhere-synthesized model).")
    ap.add_argument("--report-json", default=None, type=Path)
    args = ap.parse_args()

    if args.inp is not None:
        # Lazy sibling import: the script's own dir is on sys.path when run
        # directly, so this resolves without packaging the scripts dir.
        from inp_to_network import inp_to_network
        network = inp_to_network(args.inp)
    elif args.network_json is not None:
        network = load_json(args.network_json)
    else:
        ap.error("provide a network_json path or --inp <model.inp>")

    report = run_qa(network)
    if args.report_json:
        save_json(args.report_json, report)
    print(json.dumps(report, indent=2))


if __name__ == "__main__":
    main()

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

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


def load_json(path: str | Path) -> Any:
    return json.loads(Path(path).read_text())


def format_num(x: Any) -> str:
    if x is None:
        return ""
    if isinstance(x, bool):
        return "YES" if x else "NO"
    if isinstance(x, int):
        return str(x)
    if isinstance(x, float):
        return f"{x:.6f}".rstrip("0").rstrip(".")
    return str(x)


def emit_junctions(network: dict) -> list[str]:
    lines = ["[JUNCTIONS]", ";;Name           Elevation      MaxDepth       InitDepth      SurDepth       Aponded"]
    for j in network.get("junctions", []):
        lines.append(
            f"{j['id']:<16} {format_num(j['invert_elev']):<14} {format_num(j['max_depth']):<14} {format_num(j.get('init_depth', 0)):<14} {format_num(j.get('sur_depth', 0)):<14} {format_num(j.get('aponded', 0))}"
        )
    return lines


def emit_outfalls(network: dict) -> list[str]:
    lines = ["[OUTFALLS]", ";;Name           Elevation      Type           Stage Data      Gated          Route To"]
    for o in network.get("outfalls", []):
        lines.append(
            f"{o['id']:<16} {format_num(o['invert_elev']):<14} {o['type']:<14} {format_num(o.get('stage_data', '')):<15} {format_num(o.get('gated', False)):<14} {format_num(o.get('route_to', ''))}"
        )
    return lines


def emit_conduits(network: dict) -> list[str]:
    lines = ["[CONDUITS]", ";;Name           From Node       To Node         Length         Roughness      InOffset       OutOffset      InitFlow       MaxFlow"]
    for c in network.get("conduits", []):
        lines.append(
            f"{c['id']:<16} {c['from_node']:<15} {c['to_node']:<15} {format_num(c['length']):<14} {format_num(c['roughness']):<14} {format_num(c.get('in_offset', 0)):<14} {format_num(c.get('out_offset', 0)):<14} {format_num(c.get('init_flow', 0)):<14} {format_num(c.get('max_flow', ''))}"
        )
    return lines


def emit_xsections(network: dict) -> list[str]:
    lines = ["[XSECTIONS]", ";;Link           Shape           Geom1          Geom2          Geom3          Geom4          Barrels"]
    for c in network.get("conduits", []):
        xs = c["xsection"]
        lines.append(
            f"{c['id']:<16} {xs['shape']:<15} {format_num(xs['geom1']):<14} {format_num(xs.get('geom2', 0)):<14} {format_num(xs.get('geom3', 0)):<14} {format_num(xs.get('geom4', 0)):<14} {format_num(xs.get('barrels', 1))}"
        )
    return lines


def emit_coordinates(network: dict) -> list[str]:
    lines = ["[COORDINATES]", ";;Node           X-Coord         Y-Coord"]
    for j in network.get("junctions", []):
        xy = j["coordinates"]
        lines.append(f"{j['id']:<16} {format_num(xy['x']):<15} {format_num(xy['y'])}")
    for o in network.get("outfalls", []):
        xy = o["coordinates"]
        lines.append(f"{o['id']:<16} {format_num(xy['x']):<15} {format_num(xy['y'])}")
    return lines


def emit_vertices(network: dict) -> list[str]:
    lines = ["[VERTICES]", ";;Link           X-Coord         Y-Coord"]
    found = False
    for c in network.get("conduits", []):
        for v in c.get("vertices", []) or []:
            found = True
            lines.append(f"{c['id']:<16} {format_num(v['x']):<15} {format_num(v['y'])}")
    return lines if found else []


def render_inp(network: dict) -> str:
    blocks = [
        emit_junctions(network),
        emit_outfalls(network),
        emit_conduits(network),
        emit_xsections(network),
        emit_coordinates(network),
    ]
    vertices = emit_vertices(network)
    if vertices:
        blocks.append(vertices)
    return "\n\n".join("\n".join(block) for block in blocks) + "\n"


def main() -> None:
    ap = argparse.ArgumentParser()
    ap.add_argument("network_json", type=Path)
    ap.add_argument("--out", default=None, type=Path)
    args = ap.parse_args()

    network = load_json(args.network_json)
    text = render_inp(network)
    if args.out:
        args.out.parent.mkdir(parents=True, exist_ok=True)
        args.out.write_text(text, encoding="utf-8")
    print(text)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Clip raw municipal storm shapefiles to a basin and emit adapter-ready GeoJSON.

Bridges raw `StormGravityMain.shp` (+ optional `StormManhole.shp`) into the
GeoJSON inputs that `swmm-network-mcp.import_city_network` expects, while
also filling a `mapping.json` from a template.

This tool DOES NOT pick the outfall (see `infer_outfall.py`) and DOES NOT
reorient pipes by flow direction (see `reorient_pipes.py`). It performs
clip + reproject-check + template fill only.
"""
from __future__ import annotations

import argparse
import json
import sys
from pathlib import Path

import geopandas as gpd
from shapely.geometry import shape

from _hash_util import sha256_file


def load_basin_polygon(path: Path) -> tuple[gpd.GeoDataFrame, "CRS"]:
    gdf = gpd.read_file(path)
    if len(gdf) == 0:
        raise ValueError(f"basin clip is empty: {path}")
    if gdf.crs is None:
        raise ValueError(f"basin clip has no CRS: {path}")
    return gdf, gdf.crs


def clip_to_basin(
    layer_path: Path,
    basin_gdf: gpd.GeoDataFrame,
    basin_crs,
    layer_label: str,
) -> gpd.GeoDataFrame:
    gdf = gpd.read_file(layer_path)
    if gdf.crs is None:
        raise ValueError(f"{layer_label} has no CRS: {layer_path}")
    if gdf.crs != basin_crs:
        gdf = gdf.to_crs(basin_crs)
    clipped = gpd.clip(gdf, basin_gdf)
    # gpd.clip can leave empty geometries when nothing intersects; drop them
    clipped = clipped[~clipped.geometry.is_empty & clipped.geometry.notna()].copy()
    clipped.reset_index(drop=True, inplace=True)
    return clipped


def fill_mapping_template(
    template_path: Path,
    case_name: str,
    source_description: str,
    diameter_policy: str | None,
) -> dict:
    mapping = json.loads(template_path.read_text(encoding="utf-8"))
    meta = mapping.setdefault("meta", {})
    meta["name"] = case_name
    meta["source"] = source_description
    if diameter_policy is not None:
        meta["diameter_policy"] = diameter_policy
    return mapping


def parse_args() -> argparse.Namespace:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--pipes-shp", required=True, help="Path to storm pipe LineString shapefile.")
    ap.add_argument("--manholes-shp", default=None, help="Optional Point shapefile of manholes / junctions.")
    ap.add_argument("--basin-clip", required=True, help="GeoJSON polygon defining the basin to clip to.")
    ap.add_argument("--mapping-template", required=True, help="city_mapping_raw_shapefile.template.json path.")
    ap.add_argument("--out-dir", required=True, help="Output directory; will be created if missing.")
    ap.add_argument("--case-name", required=True, help="Value to write into meta.name of mapping.json.")
    ap.add_argument(
        "--source-description",
        required=True,
        help="Value to write into meta.source of mapping.json (e.g. 'Saanich StormWaterSHP subset').",
    )
    ap.add_argument(
        "--diameter-policy",
        default=None,
        help="Optional value for meta.diameter_policy in mapping.json.",
    )
    return ap.parse_args()


def main() -> None:
    args = parse_args()
    pipes_path = Path(args.pipes_shp)
    basin_path = Path(args.basin_clip)
    template_path = Path(args.mapping_template)
    out_dir = Path(args.out_dir)
    manholes_path = Path(args.manholes_shp) if args.manholes_shp else None

    for required in (pipes_path, basin_path, template_path):
        if not required.exists():
            raise FileNotFoundError(required)
    if manholes_path is not None and not manholes_path.exists():
        raise FileNotFoundError(manholes_path)

    out_dir.mkdir(parents=True, exist_ok=True)
    pipes_out = out_dir / "pipes.geojson"
    manholes_out = out_dir / "manholes.geojson" if manholes_path else None
    mapping_out = out_dir / "mapping.json"

    basin_gdf, basin_crs = load_basin_polygon(basin_path)

    pipes_clipped = clip_to_basin(pipes_path, basin_gdf, basin_crs, "pipes")
    if len(pipes_clipped) == 0:
        raise ValueError(
            f"no pipe features intersect basin {basin_path} (clipped count = 0). "
            "Check CRS consistency and that the basin polygon overlaps the pipe layer."
        )
    pipes_clipped.to_file(pipes_out, driver="GeoJSON")

    manholes_count: int | None = None
    if manholes_path is not None and manholes_out is not None:
        manholes_clipped = clip_to_basin(manholes_path, basin_gdf, basin_crs, "manholes")
        manholes_count = len(manholes_clipped)
        manholes_clipped.to_file(manholes_out, driver="GeoJSON")

    mapping = fill_mapping_template(
        template_path,
        case_name=args.case_name,
        source_description=args.source_description,
        diameter_policy=args.diameter_policy,
    )
    mapping_out.write_text(json.dumps(mapping, indent=2), encoding="utf-8")

    summary = {
        "ok": True,
        "skill": "swmm-network",
        "tool": "prepare_storm_inputs",
        "case_name": args.case_name,
        "basin_crs": str(basin_crs),
        "counts": {
            "pipes_clipped": len(pipes_clipped),
            "manholes_clipped": manholes_count,
        },
        "outputs": {
            "pipes_geojson": str(pipes_out),
            "manholes_geojson": str(manholes_out) if manholes_out else None,
            "mapping_json": str(mapping_out),
        },
        "inputs": {
            "pipes_shp": str(pipes_path),
            "manholes_shp": str(manholes_path) if manholes_path else None,
            "basin_clip": str(basin_path),
            "mapping_template": str(template_path),
        },
        "input_hashes": {
            "basin_clip_sha256": sha256_file(basin_path),
            "mapping_template_sha256": sha256_file(template_path),
        },
    }
    print(json.dumps(summary, indent=2))


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

#!/usr/bin/env python3
"""Reorient LineString pipes so geometry direction matches flow direction.

Raw municipal storm shapefiles store pipes as LineStrings whose vertex
order reflects digitisation, not flow. The downstream `city_network_adapter`
treats the first endpoint as `from_node` and the last as `to_node`, so
misoriented pipes produce a network that the QA step flags with
`no_outfall_path` warnings.

Algorithm: a breadth-first walk starting from the outfall vertices.
For each pipe touching the current frontier vertex, ensure its `to_node`
end matches the frontier vertex. If not, reverse the LineString. Then
the other endpoint joins the frontier as a newly-known downstream vertex.

Pipes that have no path to any outfall are left with their original
orientation and reported in `unreached_pipes`.
"""
from __future__ import annotations

import argparse
import json
import sys
from collections import deque
from pathlib import Path
from typing import Iterable


def _round_key(x: float, y: float, precision: int) -> tuple[float, float]:
    return (round(x, precision), round(y, precision))


def _line_endpoints(coords: list[list[float]]) -> tuple[tuple[float, float], tuple[float, float]]:
    if len(coords) < 2:
        raise ValueError("LineString must have at least 2 vertices")
    head = (float(coords[0][0]), float(coords[0][1]))
    tail = (float(coords[-1][0]), float(coords[-1][1]))
    return head, tail


def _key(point: tuple[float, float], precision: int) -> tuple[float, float]:
    return _round_key(point[0], point[1], precision)


def _reverse_line(coords: list[list[float]]) -> list[list[float]]:
    return list(reversed(coords))


def _outfall_points(features: list[dict]) -> list[tuple[float, float]]:
    points: list[tuple[float, float]] = []
    for f in features:
        geom = f.get("geometry") or {}
        if geom.get("type") != "Point":
            continue
        c = geom.get("coordinates") or []
        if len(c) < 2:
            continue
        points.append((float(c[0]), float(c[1])))
    return points


def _nearest_pipe_vertex(
    target: tuple[float, float],
    pipe_endpoints: Iterable[tuple[tuple[float, float], tuple[float, float]]],
    precision: int,
) -> tuple[float, float] | None:
    """Pick the pipe vertex closest to an outfall point. Returns the rounded key."""
    best: tuple[float, float] | None = None
    best_sq = float("inf")
    for head, tail in pipe_endpoints:
        for endpoint in (head, tail):
            dx = endpoint[0] - target[0]
            dy = endpoint[1] - target[1]
            sq = dx * dx + dy * dy
            if sq < best_sq:
                best_sq = sq
                best = endpoint
    if best is None:
        return None
    return _key(best, precision)


def reorient(
    pipes_geojson: dict,
    outfalls_geojson: dict,
    precision: int = 3,
) -> tuple[dict, dict]:
    features = list(pipes_geojson.get("features") or [])
    outfalls = _outfall_points(list(outfalls_geojson.get("features") or []))
    if not features:
        raise ValueError("pipes geojson has no features")
    if not outfalls:
        raise ValueError("outfalls geojson has no Point features")

    # Map vertex key -> list of pipe indices connected at either endpoint
    vertex_to_pipes: dict[tuple[float, float], list[int]] = {}
    pipe_endpoints: list[tuple[tuple[float, float], tuple[float, float]]] = []
    for idx, feat in enumerate(features):
        geom = feat.get("geometry") or {}
        if geom.get("type") != "LineString":
            raise ValueError(f"pipe feature {idx} is not a LineString")
        head, tail = _line_endpoints(geom["coordinates"])
        pipe_endpoints.append((head, tail))
        head_key = _key(head, precision)
        tail_key = _key(tail, precision)
        vertex_to_pipes.setdefault(head_key, []).append(idx)
        vertex_to_pipes.setdefault(tail_key, []).append(idx)

    # Seed BFS frontier with the pipe vertex closest to each outfall
    frontier: deque[tuple[float, float]] = deque()
    seeded: set[tuple[float, float]] = set()
    for ofall in outfalls:
        seed = _nearest_pipe_vertex(ofall, pipe_endpoints, precision)
        if seed is not None and seed not in seeded:
            seeded.add(seed)
            frontier.append(seed)

    if not frontier:
        raise ValueError("could not seed BFS frontier — no pipe vertex near any outfall")

    visited_vertices: set[tuple[float, float]] = set(seeded)
    visited_pipes: set[int] = set()
    reversed_indices: list[int] = []

    while frontier:
        vertex = frontier.popleft()
        for pipe_idx in vertex_to_pipes.get(vertex, []):
            if pipe_idx in visited_pipes:
                continue
            visited_pipes.add(pipe_idx)
            head, tail = pipe_endpoints[pipe_idx]
            head_key = _key(head, precision)
            tail_key = _key(tail, precision)
            # The pipe must flow INTO the current vertex (to_node == vertex).
            if tail_key == vertex:
                upstream = head_key
            elif head_key == vertex:
                # Need to flip so that the line ends at the current vertex.
                coords = features[pipe_idx]["geometry"]["coordinates"]
                features[pipe_idx]["geometry"]["coordinates"] = _reverse_line(coords)
                pipe_endpoints[pipe_idx] = (tail, head)
                reversed_indices.append(pipe_idx)
                upstream = tail_key
            else:
                # Pipe touches neither endpoint of the visited vertex; defensive skip.
                continue
            if upstream not in visited_vertices:
                visited_vertices.add(upstream)
                frontier.append(upstream)

    unreached = [
        {
            "index": idx,
            "id": features[idx].get("properties", {}).get("FACILITYID")
            or features[idx].get("properties", {}).get("id")
            or f"pipe_{idx}",
        }
        for idx in range(len(features))
        if idx not in visited_pipes
    ]

    out = {
        "type": "FeatureCollection",
        "name": pipes_geojson.get("name", "pipes_oriented"),
        "features": features,
    }
    if "crs" in pipes_geojson:
        out["crs"] = pipes_geojson["crs"]

    report = {
        "ok": True,
        "skill": "swmm-network",
        "tool": "reorient_pipes",
        "counts": {
            "pipes_total": len(features),
            "pipes_reversed": len(reversed_indices),
            "pipes_reached": len(visited_pipes),
            "pipes_unreached": len(unreached),
            "outfalls_used": len(outfalls),
        },
        "reversed_pipe_indices": reversed_indices,
        "unreached_pipes": unreached,
        "coordinate_precision": precision,
    }
    return out, report


def parse_args() -> argparse.Namespace:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--pipes-geojson", required=True)
    ap.add_argument("--outfalls-geojson", required=True)
    ap.add_argument("--out", required=True, help="Output path for reoriented pipes geojson.")
    ap.add_argument(
        "--coordinate-precision",
        type=int,
        default=3,
        help="Decimal places used to match pipe endpoints to vertices.",
    )
    return ap.parse_args()


def main() -> None:
    args = parse_args()
    pipes_path = Path(args.pipes_geojson)
    outfalls_path = Path(args.outfalls_geojson)
    out_path = Path(args.out)
    for p in (pipes_path, outfalls_path):
        if not p.exists():
            raise FileNotFoundError(p)

    pipes = json.loads(pipes_path.read_text(encoding="utf-8"))
    outfalls = json.loads(outfalls_path.read_text(encoding="utf-8"))

    oriented, report = reorient(pipes, outfalls, precision=args.coordinate_precision)
    out_path.parent.mkdir(parents=True, exist_ok=True)
    out_path.write_text(json.dumps(oriented, indent=2), encoding="utf-8")
    report["outputs"] = {"pipes_oriented_geojson": str(out_path)}
    report["inputs"] = {
        "pipes_geojson": str(pipes_path),
        "outfalls_geojson": str(outfalls_path),
    }
    print(json.dumps(report, indent=2))


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

#!/usr/bin/env python3
"""Snap nearby pipe LineString endpoints together so the network is graph-connected.

Real municipal storm pipe shapefiles are usually hand-digitised. Adjacent
pipe segments that should share a manhole often have endpoints separated
by sub-millimetre to several-centimetre vertex drift. The
``city_network_adapter`` infers junctions by exact-coordinate equality
and so treats those drifting endpoints as separate nodes — leaving the
pipe network as a forest of disconnected fragments.

This tool clusters every pipe endpoint within a tolerance and snaps each
cluster to the cluster centroid, rewriting the LineStrings so that
shared endpoints land on identical coordinates.

Algorithm: union-find over endpoints. Two endpoints are merged when
their Euclidean distance is below ``tolerance_m``. Each connected
component's centroid becomes the snapped coordinate. The interior
vertices of each LineString are left untouched; only the first and last
vertices are rewritten.
"""
from __future__ import annotations

import argparse
import json
import math
import sys
from pathlib import Path


def _euclid(a: tuple[float, float], b: tuple[float, float]) -> float:
    return math.hypot(a[0] - b[0], a[1] - b[1])


class _UnionFind:
    def __init__(self, n: int) -> None:
        self.parent = list(range(n))

    def find(self, i: int) -> int:
        while self.parent[i] != i:
            self.parent[i] = self.parent[self.parent[i]]
            i = self.parent[i]
        return i

    def union(self, i: int, j: int) -> None:
        ri, rj = self.find(i), self.find(j)
        if ri != rj:
            self.parent[ri] = rj


def snap(pipes_geojson: dict, tolerance_m: float) -> tuple[dict, dict]:
    feats = list(pipes_geojson.get("features") or [])
    if not feats:
        raise ValueError("pipes geojson has no features")
    if tolerance_m < 0:
        raise ValueError("tolerance_m must be non-negative")

    # Collect every (feature_index, endpoint_position, point) record.
    endpoints: list[tuple[int, str, tuple[float, float]]] = []
    for fi, f in enumerate(feats):
        geom = f.get("geometry") or {}
        if geom.get("type") != "LineString":
            raise ValueError(f"feature {fi} is not a LineString")
        coords = geom.get("coordinates") or []
        if len(coords) < 2:
            raise ValueError(f"feature {fi} has < 2 vertices")
        endpoints.append((fi, "start", (float(coords[0][0]), float(coords[0][1]))))
        endpoints.append((fi, "end", (float(coords[-1][0]), float(coords[-1][1]))))

    n = len(endpoints)
    uf = _UnionFind(n)

    # Bucket endpoints into a grid of size tolerance_m for O(n) clustering.
    # When tolerance_m == 0 there is nothing to snap; just return the input
    # untouched.
    if tolerance_m == 0:
        report = {
            "ok": True,
            "skill": "swmm-network",
            "tool": "snap_pipe_endpoints",
            "tolerance_m": 0.0,
            "counts": {
                "endpoints_total": n,
                "clusters": n,
                "clusters_merged": 0,
                "max_snap_distance_m": 0.0,
            },
        }
        return pipes_geojson, report

    cell = float(tolerance_m)
    buckets: dict[tuple[int, int], list[int]] = {}
    for idx, (_, _, pt) in enumerate(endpoints):
        key = (int(math.floor(pt[0] / cell)), int(math.floor(pt[1] / cell)))
        buckets.setdefault(key, []).append(idx)

    # For each endpoint, scan the 3x3 neighborhood of buckets for partners.
    for idx, (_, _, pt) in enumerate(endpoints):
        kx = int(math.floor(pt[0] / cell))
        ky = int(math.floor(pt[1] / cell))
        for dx in (-1, 0, 1):
            for dy in (-1, 0, 1):
                for jdx in buckets.get((kx + dx, ky + dy), []):
                    if jdx <= idx:
                        continue
                    if _euclid(pt, endpoints[jdx][2]) <= tolerance_m:
                        uf.union(idx, jdx)

    # Collect cluster centroids and per-cluster max distance.
    clusters: dict[int, list[int]] = {}
    for idx in range(n):
        clusters.setdefault(uf.find(idx), []).append(idx)

    snapped_xy: dict[int, tuple[float, float]] = {}
    max_dist = 0.0
    clusters_merged = 0
    for _, members in clusters.items():
        if len(members) < 2:
            snapped_xy[members[0]] = endpoints[members[0]][2]
            continue
        clusters_merged += 1
        xs = [endpoints[m][2][0] for m in members]
        ys = [endpoints[m][2][1] for m in members]
        cx = sum(xs) / len(xs)
        cy = sum(ys) / len(ys)
        for m in members:
            d = _euclid(endpoints[m][2], (cx, cy))
            if d > max_dist:
                max_dist = d
            snapped_xy[m] = (cx, cy)

    # Rewrite the geojson features in place (deep copy not required —
    # caller doesn't reuse pipes_geojson).
    for idx, (fi, pos, _) in enumerate(endpoints):
        new_xy = snapped_xy[idx]
        coords = feats[fi]["geometry"]["coordinates"]
        if pos == "start":
            coords[0] = [new_xy[0], new_xy[1]]
        else:
            coords[-1] = [new_xy[0], new_xy[1]]

    # Drop pipes whose two endpoints landed in the same cluster — they
    # would be self-loop conduits that SWMM and import_city_network
    # both reject. Common causes: an already-tiny pipe in the source
    # data, or a snap tolerance large enough to collapse a short
    # pipe's two ends.
    kept: list[dict] = []
    dropped_self_loops: list[dict] = []
    for fi, f in enumerate(feats):
        coords = f["geometry"]["coordinates"]
        head = (coords[0][0], coords[0][1])
        tail = (coords[-1][0], coords[-1][1])
        if head == tail:
            dropped_self_loops.append({
                "index": fi,
                "id": (f.get("properties") or {}).get("FACILITYID")
                    or (f.get("properties") or {}).get("id")
                    or f"pipe_{fi}",
                "snapped_endpoint": [head[0], head[1]],
            })
        else:
            kept.append(f)

    out = {
        "type": "FeatureCollection",
        "name": pipes_geojson.get("name", "pipes_snapped"),
        "features": kept,
    }
    if "crs" in pipes_geojson:
        out["crs"] = pipes_geojson["crs"]

    report = {
        "ok": True,
        "skill": "swmm-network",
        "tool": "snap_pipe_endpoints",
        "tolerance_m": float(tolerance_m),
        "counts": {
            "endpoints_total": n,
            "clusters": len(clusters),
            "clusters_merged": clusters_merged,
            "max_snap_distance_m": float(max_dist),
            "pipes_in": len(feats),
            "pipes_out": len(kept),
            "pipes_dropped_as_self_loops": len(dropped_self_loops),
        },
        "dropped_self_loops": dropped_self_loops,
    }
    return out, report


def parse_args() -> argparse.Namespace:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--pipes-geojson", required=True)
    ap.add_argument("--tolerance-m", type=float, required=True,
                    help="Maximum distance (in CRS units, expected metres) to merge two endpoints.")
    ap.add_argument("--out", required=True)
    return ap.parse_args()


def main() -> None:
    args = parse_args()
    pipes_path = Path(args.pipes_geojson)
    out_path = Path(args.out)
    if not pipes_path.exists():
        raise FileNotFoundError(pipes_path)

    pipes = json.loads(pipes_path.read_text(encoding="utf-8"))
    snapped, report = snap(pipes, args.tolerance_m)
    out_path.parent.mkdir(parents=True, exist_ok=True)
    out_path.write_text(json.dumps(snapped, indent=2), encoding="utf-8")
    report["outputs"] = {"pipes_snapped_geojson": str(out_path)}
    report["inputs"] = {"pipes_geojson": str(pipes_path)}
    print(json.dumps(report, indent=2))


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

"""Shared SHA-256 file hashing for skills/swmm-network/scripts/.

``prepare_storm_inputs.py`` and ``city_network_adapter.py`` each
hand-rolled the same chunked hashlib loop. Converged here per
ADR-0006 D5. Deliberately agentic_swmm-import-free: these scripts run
as standalone subprocess entry points spawned by mcp/swmm-network's
server.js, not through the agentic_swmm package.
"""

from __future__ import annotations

import hashlib
from pathlib import Path


def sha256_file(path: Path | None) -> str | None:
    """Hex SHA-256 digest of ``path``'s bytes, or ``None`` if ``path`` is ``None``."""
    if path is None:
        return None
    digest = hashlib.sha256()
    with path.open("rb") as f:
        for chunk in iter(lambda: f.read(1024 * 1024), b""):
            digest.update(chunk)
    return digest.hexdigest()

{
  "meta": {
    "name": "<case-name>",
    "source": "<municipal storm shapefile origin, e.g. Saanich StormWaterSHP>",
    "evidence_boundary": "Imported from raw municipal LineString pipe layer + Point manhole layer. Junctions are inferred from pipe endpoints when not supplied; outfall must be supplied separately. Diameter / invert fallbacks declared in defaults below.",
    "diameter_policy": "<describe: e.g. ignore raw DIAMETER when nonnumeric; substitute defaults.geom1>"
  },
  "dual_system_ready": false,
  "coordinate_precision": 3,
  "inference": {
    "junction_prefix": "J_AUTO",
    "max_depth": 2.5
  },
  "pipes": {
    "fields": {
      "id": "FACILITYID",
      "shape": "CRSECSHAPE",
      "material": "MATERIAL",
      "source_id": "FACILITYID"
    },
    "defaults": {
      "shape": "CIRCULAR",
      "roughness": 0.013,
      "geom1": 0.3,
      "geom2": 0.0,
      "geom3": 0.0,
      "geom4": 0.0,
      "barrels": 1,
      "in_offset": 0.0,
      "out_offset": 0.0,
      "init_flow": 0.0,
      "max_flow": null,
      "minimum_length": 1.0,
      "asset_type": "storm",
      "system_layer": "minor_pipe",
      "from_invert_elev": 0.0,
      "to_invert_elev": 0.0
    }
  },
  "junctions": {
    "fields": {
      "id": "node_id",
      "x": "x",
      "y": "y",
      "invert_elev": "invert_elev",
      "max_depth": "max_depth",
      "asset_type": "asset_type",
      "system_layer": "system_layer"
    },
    "defaults": {
      "invert_elev": 0.0,
      "max_depth": 2.5,
      "asset_type": "storm",
      "system_layer": "minor_pipe"
    }
  },
  "outfalls": {
    "fields": {
      "id": "node_id",
      "x": "x",
      "y": "y",
      "invert_elev": "invert_elev",
      "type": "type",
      "asset_type": "asset_type",
      "system_layer": "system_layer"
    },
    "defaults": {
      "type": "FREE",
      "gated": false,
      "asset_type": "storm",
      "system_layer": "minor_pipe"
    }
  }
}

# `mapping.json` templates for `import_city_network`

The `swmm-network-mcp.import_city_network` tool (and the underlying
`city_network_adapter.py` script) require a `mapping.json` config that
describes how the raw input columns map to the SWMM network model.

There are two practical shapes of `mapping.json` because there are two
practical shapes of input data. **Pick the one that matches what you
actually have, then edit it.**

---

## Style 1 — structured / dual-system export (CAD or asset-DB origin)

Use when your pipes table has **explicit `from_node` / `to_node` columns**
and per-pipe `from_x / from_y / to_x / to_y` and invert elevations. This is
typical of structured CAD or asset-DB exports.

**Template:** `../examples/city-dual-system/mapping.json`
**Example input:** `../examples/city-dual-system/pipes.csv` and friends.

Salient features:

- `pipes.fields` enumerates `from_node`, `to_node`, `from_x`, `from_y`,
  `to_x`, `to_y`, `from_invert_elev`, `to_invert_elev`, `diameter`,
  `roughness`, `material` — the adapter trusts these directly.
- `dual_system_ready: true` is allowed and meaningful.
- Junction inference is rarely used (junctions usually come in as their
  own explicit CSV).

---

## Style 2 — raw municipal storm shapefile (LineString pipes only)

Use when your pipes geojson is a **bare LineString layer** (typical of a
municipal `StormGravityMain.shp` export). Pipes have no `from_node` /
`to_node` columns; coordinates only exist as geometry vertices; invert
elevations are usually absent; some attribute fields (e.g. `DIAMETER`)
may contain nonnumeric values.

**Template:** `city_mapping_raw_shapefile.template.json`
**Reference real-world use:** `docs/framework-validation/saanich-smoke-20260513/`
contains a working mapping derived from this template for Saanich data.

Salient features:

- `pipes.fields` is **sparse on purpose** — only the few fields the
  shapefile actually provides (`FACILITYID`, `CRSECSHAPE`, `MATERIAL`).
  The adapter then **infers junctions from pipe endpoints** rather than
  looking them up.
- `pipes.defaults.geom1` is the fallback diameter applied when the raw
  DIAMETER field is missing or nonnumeric. Document the policy in
  `meta.diameter_policy` so the manifest captures it.
- `pipes.defaults.from_invert_elev / to_invert_elev` default to 0.0 m.
  SWMM will warn about elevation drops; this is intentional for smoke
  runs and should be replaced by a DEM-based invert inference once
  available (see `BACKLOG.md F12`).
- `dual_system_ready: false`.
- `inference.junction_prefix: "J_AUTO"` causes inferred junctions to be
  named `J_AUTO_<x>p<y>_<...>`.

### What this template does NOT do

Two known gaps live outside the mapping config and need separate
remediation:

1. **Outfall identification.** The mapping config does not pick the
   outfall node; you must supply a separate `outfalls.geojson` to
   `import_city_network`. See `BACKLOG.md B3`.
2. **Pipe orientation.** The adapter treats pipe geometry direction as
   flow direction. For raw municipal shapefiles this is rarely true.
   The first `qa` call will flag `no_outfall_path` warnings. See
   `BACKLOG.md B5`.

When (B3) and (B5) are closed, those steps will be invoked by a separate
MCP tool and remain orthogonal to this mapping config.

---

## How to use a template

1. Copy the template into your run directory (e.g.
   `runs/<case>/04_network/mapping.json`).
2. Replace placeholder strings under `meta.*` with your case-specific
   strings.
3. If your shapefile uses different field names than the defaults
   (e.g. `OBJECTID` rather than `FACILITYID`), edit `pipes.fields.*`.
4. Adjust `defaults.geom1` to a sensible diameter for unknown pipes in
   your dataset.
5. Pass the file via `--arguments-json '{"mappingPath": "<path>"}'` when
   invoking the MCP tool.

The adapter does not validate the mapping config against a schema yet;
typos in `fields.*` will surface as "field not found" errors at import
time.