"""Build the route-binding network resource from the Rijkswaterstaat FIS-VNDS
``route`` layer — and derive its connectivity ourselves.

The ``route`` layer is the source of truth for the route MANAGER and vessel
BINDING: full-stretch, properly-named fairway routes (e.g. id 20940 =
'Westerschelde', ~105 km) carrying linear-referencing km. The fine-grained
``vaarwegen`` passages stay a display-only chart layer.

We deliberately do NOT use the service's own ``routejunctie`` layer. Instead we
pass over the route geometries and flag a connection only where two routes
actually meet:

  * SHARED VERTEX : a location where vertices of >=2 distinct routes coincide
                    within ``vertex_tol_m``. Routes are digitised to share exact
                    vertices at connections, so ~2 m is ample and stable.
  * T-BRANCH      : a route endpoint with no shared vertex that lands within
                    ``tbranch_tol_m`` of another route's *interior* — projected
                    onto that route and flagged (the few real T-junctions).

The result is baked into a bundled GeoPackage (``data/routes.gpkg``: a ``route``
lines table + a ``junctions`` points table) so binding runs fully offline.
Re-run ``build_resource()`` to refresh from the live service.
"""
import json
import math
import os
import urllib.request
from collections import defaultdict

SERVICE_URL = (
    "https://geo.rijkswaterstaat.nl/services/ogc/gdr/fis_vnds/ows"
    "?service=WFS&version=2.0.0&request=GetFeature"
    "&typeNames=fis_vnds:route&outputFormat=application/json&count=50000"
)

_MDEG = 111320.0


# --------------------------------------------------------------------------- #
#  Fetch + parse
# --------------------------------------------------------------------------- #
def fetch_routes_geojson(url=SERVICE_URL, timeout=180):
    """GET the full ``route`` layer as a GeoJSON FeatureCollection (dict)."""
    if not str(url).lower().startswith("https://"):
        raise ValueError("refusing non-HTTPS URL: %r" % (url,))
    req = urllib.request.Request(url, headers={"User-Agent": "ais-watchkeeper"})
    with urllib.request.urlopen(req, timeout=timeout) as r:   # nosec B310 - scheme checked (https only)
        return json.loads(r.read().decode("utf-8"))


def parse_routes(geojson):
    """GeoJSON -> list of route dicts.

    Each dict: ``{routeid, name, watername, code, kmbegin, kmend, coords}`` with
    ``coords`` a list of (lon, lat). A MultiLineString becomes one entry per
    part (same attributes), matching how the network engine indexes segments.
    """
    out = []
    for f in geojson.get("features", []):
        g = f.get("geometry") or {}
        p = f.get("properties") or {}
        t = g.get("type")
        if t == "LineString":
            parts = [g["coordinates"]]
        elif t == "MultiLineString":
            parts = g["coordinates"]
        else:
            continue
        for part in parts:
            if len(part) >= 2:
                out.append({
                    "routeid": p.get("id"),
                    "name": p.get("name") or "",
                    "watername": p.get("watername") or "",
                    "code": p.get("code") or "",
                    "kmbegin": p.get("routekmbegin"),
                    "kmend": p.get("routekmend"),
                    "coords": [(float(c[0]), float(c[1])) for c in part],
                })
    return out


# --------------------------------------------------------------------------- #
#  Geometry helpers (local equirectangular metric frame)
# --------------------------------------------------------------------------- #
def _frame(lines):
    """Mean-latitude metric frame: returns (to_m, to_lonlat) closures."""
    n = sum(len(ln["coords"]) for ln in lines) or 1
    lat0 = sum(la for ln in lines for _, la in ln["coords"]) / n
    cos0 = math.cos(math.radians(lat0))
    return (lambda lo, la: (lo * _MDEG * cos0, la * _MDEG),
            lambda x, y: (x / (_MDEG * cos0), y / _MDEG))


def _proj_point_on_seg(px, py, ax, ay, bx, by):
    """(distance, t, cx, cy): nearest point on segment a->b to p; t in [0,1]."""
    dx, dy = bx - ax, by - ay
    L2 = dx * dx + dy * dy
    if L2 == 0.0:
        return math.hypot(px - ax, py - ay), 0.0, ax, ay
    t = max(0.0, min(1.0, ((px - ax) * dx + (py - ay) * dy) / L2))
    cx, cy = ax + t * dx, ay + t * dy
    return math.hypot(px - cx, py - cy), t, cx, cy


# --------------------------------------------------------------------------- #
#  Junction derivation
# --------------------------------------------------------------------------- #
def derive_junctions(lines, vertex_tol_m=2.0, tbranch_tol_m=10.0):
    """Flag connection points across the route network.

    Returns a list of (lon, lat). A point is flagged when >=2 distinct routes
    share a vertex within ``vertex_tol_m`` (touch / branch / multi-way junction),
    plus the few dangling endpoints that sit within ``tbranch_tol_m`` of another
    route's interior (real T-branches), projected onto that route.
    """
    if not lines:
        return []
    to_m, to_ll = _frame(lines)

    # ---- shared-vertex clusters ----
    cell = defaultdict(list)                 # gridcell -> [(routeid, x, y), ...]
    for ln in lines:
        rid = ln["routeid"]
        for lo, la in ln["coords"]:
            x, y = to_m(lo, la)
            cell[(round(x / vertex_tol_m), round(y / vertex_tol_m))].append(
                (rid, x, y))
    shared_cells = set()
    junctions = []
    for k, items in cell.items():
        if len({r for r, _, _ in items}) >= 2:
            shared_cells.add(k)
            mx = sum(x for _, x, _ in items) / len(items)
            my = sum(y for _, _, y in items) / len(items)
            junctions.append(to_ll(mx, my))

    # ---- T-branches: dangling endpoints near another route's interior ----
    scell = 200.0
    seg = defaultdict(list)
    for ln in lines:
        m = [to_m(lo, la) for lo, la in ln["coords"]]
        for i in range(len(m) - 1):
            ax, ay = m[i]; bx, by = m[i + 1]
            c0 = (int(ax // scell), int(ay // scell))
            c1 = (int(bx // scell), int(by // scell))
            for cx in range(min(c0[0], c1[0]), max(c0[0], c1[0]) + 1):
                for cy in range(min(c0[1], c1[1]), max(c0[1], c1[1]) + 1):
                    seg[(cx, cy)].append((ax, ay, bx, by, ln["routeid"]))

    def is_shared(x, y):
        """True if (x, y) already lies in a shared-vertex cluster cell."""
        return (round(x / vertex_tol_m), round(y / vertex_tol_m)) in shared_cells

    n_tbranch = 0
    for ln in lines:
        for lo, la in (ln["coords"][0], ln["coords"][-1]):
            x, y = to_m(lo, la)
            if is_shared(x, y):
                continue
            cx, cy = int(x // scell), int(y // scell)
            best = None
            for ddx in (-1, 0, 1):
                for ddy in (-1, 0, 1):
                    for ax, ay, bx, by, orid in seg.get((cx + ddx, cy + ddy), ()):
                        if orid == ln["routeid"]:
                            continue
                        d, t, qx, qy = _proj_point_on_seg(x, y, ax, ay, bx, by)
                        if 0.0 < t < 1.0 and d <= tbranch_tol_m \
                                and (best is None or d < best[0]):
                            best = (d, qx, qy)
            if best is not None:
                junctions.append(to_ll(best[1], best[2]))
                n_tbranch += 1

    junctions = _dedup(junctions, to_m, vertex_tol_m)
    derive_junctions.last_stats = {
        "shared": len(shared_cells), "tbranch": n_tbranch,
        "total": len(junctions)}
    return junctions


def _dedup(points, to_m, tol_m):
    """Drop near-coincident junction points (within tol_m), keep one each."""
    seen = set()
    out = []
    for lo, la in points:
        x, y = to_m(lo, la)
        key = (round(x / tol_m), round(y / tol_m))
        if key in seen:
            continue
        seen.add(key)
        out.append((lo, la))
    return out


# --------------------------------------------------------------------------- #
#  GeoPackage writer (via OGR/GDAL)
# --------------------------------------------------------------------------- #
def write_route_gpkg(lines, junctions, path):
    """Write ``lines`` + ``junctions`` to a GeoPackage at ``path``.

    Tables: ``route`` (LineString, with routeid/name/watername/code/km) and
    ``junctions`` (Point). Columns match what ``load_network`` reads.
    """
    from osgeo import ogr, osr
    if os.path.exists(path):
        os.remove(path)
    drv = ogr.GetDriverByName("GPKG")
    ds = drv.CreateDataSource(path)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)

    rl = ds.CreateLayer("route", srs, ogr.wkbLineString)
    for nm, typ in (("routeid", ogr.OFTInteger), ("name", ogr.OFTString),
                    ("watername", ogr.OFTString), ("code", ogr.OFTString),
                    ("routekmbegin", ogr.OFTReal), ("routekmend", ogr.OFTReal)):
        rl.CreateField(ogr.FieldDefn(nm, typ))
    rd = rl.GetLayerDefn()
    for ln in lines:
        f = ogr.Feature(rd)
        rid = ln["routeid"]
        f.SetField("routeid", int(rid) if rid is not None else 0)
        f.SetField("name", ln["name"])
        f.SetField("watername", ln["watername"])
        f.SetField("code", str(ln["code"]))
        if ln["kmbegin"] is not None:
            f.SetField("routekmbegin", float(ln["kmbegin"]))
        if ln["kmend"] is not None:
            f.SetField("routekmend", float(ln["kmend"]))
        g = ogr.Geometry(ogr.wkbLineString)
        for lo, la in ln["coords"]:
            g.AddPoint_2D(lo, la)
        f.SetGeometry(g)
        rl.CreateFeature(f)
        f = None

    jl = ds.CreateLayer("junctions", srs, ogr.wkbPoint)
    jd = jl.GetLayerDefn()
    for lo, la in junctions:
        f = ogr.Feature(jd)
        g = ogr.Geometry(ogr.wkbPoint)
        g.AddPoint_2D(lo, la)
        f.SetGeometry(g)
        jl.CreateFeature(f)
        f = None
    ds = None


# --------------------------------------------------------------------------- #
#  Orchestration
# --------------------------------------------------------------------------- #
def default_gpkg_path():
    """Bundled resource path: ``<plugin>/data/routes.gpkg``."""
    return os.path.join(os.path.dirname(os.path.abspath(__file__)),
                        "data", "routes.gpkg")


def build_resource(out_path=None, geojson_path=None, vertex_tol_m=2.0,
                   tbranch_tol_m=10.0):
    """Build ``data/routes.gpkg`` from the live service (or a saved GeoJSON file).

    Returns a stats dict. Pass ``geojson_path`` to build from an already-saved
    response (offline / reproducible); otherwise the live service is fetched.
    """
    out_path = out_path or default_gpkg_path()
    if geojson_path:
        with open(geojson_path, encoding="utf-8") as fh:
            gj = json.load(fh)
    else:
        gj = fetch_routes_geojson()
    lines = parse_routes(gj)
    junctions = derive_junctions(lines, vertex_tol_m, tbranch_tol_m)
    os.makedirs(os.path.dirname(out_path), exist_ok=True)
    write_route_gpkg(lines, junctions, out_path)
    stats = {"routes": len({ln["routeid"] for ln in lines}),
             "line_parts": len(lines),
             "junctions": len(junctions),
             "out": out_path}
    stats.update(getattr(derive_junctions, "last_stats", {}))
    return stats


if __name__ == "__main__":
    import sys
    gj = sys.argv[1] if len(sys.argv) > 1 else None
    print(build_resource(geojson_path=gj))
