# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
s101/geometry.py
----------------
Native S-101 geometry assembly (S-100 Part 7 spatial schema) over the same
``iso8211`` byte reader — GDAL-free. The S-101 sibling of ``s57.geometry``.

It turns the cell's spatial records into geometry values, reusing the S-57
geometry value types (``Point``/``LineString``/``Polygon``/… with WKT output) so
downstream code and the WKT comparator are shared across both standards.

S-101 spatial model (as observed in the Belgian Ed 2.0.0 cells and S-100 Part 7):

    Point          PRID (RCNM 110)  + C2IT  (one YCOO/XCOO)
    Multipoint     MRID (RCNM 115)  + C3IL  (3-D YCOO/XCOO/ZCOO list -> MULTIPOINT Z)
    Curve          CRID (RCNM 120)  + 2x PTAS (begin TOPI=1 / end TOPI=2
                                       -> PRID) + C2IL (interior control pts)
    Composite curve CCID (RCNM 125) + CUCO (component CRIDs, each with ORNT)
    Surface        SRID (RCNM 130)  + RIAS (ring -> CCID/CRID, ORNT, USAG
                                       1=exterior / 2=interior)

A feature points at its geometry through SPAS / PTAS (RRNM = target record name,
RRID = target id, ORNT = orientation). Curve coordinates are ordered
begin-point, interior control points (C2IL), end-point. Composite-curve and ring
coordinates honour ORNT (2 = reverse) and drop the duplicated junction vertex.

Coordinates are integers scaled by the coordinate multiplication factor (COMF,
read from DSSI CMFX/CMFY, default 1e7 → 7-dp WGS84 lon/lat); soundings depth by
CMFZ (default 10). XCOO is longitude, YCOO latitude (S-57 axis order).

VERIFICATION STATUS: there is no mature GDAL S-101 oracle, so per project
discipline these geometric reads are **flagged unverified** — they are checked by
*self-consistency* (every spatial pointer resolves; curve endpoints match their
PTAS-referenced points; surface rings close; coordinates fall inside the dataset
bounding box) rather than diffed against a reference engine. See test_s101_geometry
and the headless self-consistency harness.

3-D MULTIPOINT soundings (MRID + C3IL) are assembled as of v1.26.1: the ISO 8211
``\\*`` cartesian sub-group (``VCID`` header + repeating YCOO/XCOO/ZCOO) is now
expanded by the shared reader (FieldDef.array_start), so MRID emits MULTIPOINT Z
with depth = ZCOO / SOMF (CMFZ from DSSI).
"""

try:                                                    # plugin context
    from ..s57.geometry import (Point, MultiPoint, LineString,
                                MultiLineString, Polygon, MultiPolygon)
except (ImportError, ValueError):                       # headless / top-level
    from s57.geometry import (Point, MultiPoint, LineString,
                              MultiLineString, Polygon, MultiPolygon)

RCNM_PRID = 110
RCNM_MRID = 115
RCNM_CRID = 120
RCNM_CCID = 125
RCNM_SRID = 130

ORNT_REVERSE = 2
USAG_INTERIOR = 2
TOPI_BEGIN = 1
TOPI_END = 2


def _i(v):
    if v is None:
        return None
    if isinstance(v, int):
        return v
    try:
        return int(str(v).strip())
    except (ValueError, TypeError):
        return None


def _dedup(coords):
    """Drop consecutive duplicate vertices. A curve's begin/end node (from PTAS)
    can coincide with the first/last C2IL control point, and composite-curve
    junctions repeat a vertex; either way a doubled consecutive vertex is
    redundant. Ring closure (first == last) is preserved."""
    out = []
    for c in coords:
        if not out or out[-1] != c:
            out.append(c)
    return out


class S101GeometryBuilder:
    """Assembles geometry for S-101 features from a cell's decoded ``iso8211``
    records. Construct once per cell, then call :meth:`geometry` per feature.

    Resolution is defensive: an unresolvable pointer yields ``None`` for that
    feature rather than raising, so one broken reference cannot abort a sweep.
    """

    def __init__(self, records):
        self.points = {}
        self.curves = {}
        self.composites = {}
        self.surfaces = {}
        self.multipoints = {}
        self.comf = 1e7
        self.somf = 10.0
        for rec in records:
            dssi = rec.first("DSSI")
            if dssi:
                cmfx = _i(dssi.get("CMFX"))
                cmfz = _i(dssi.get("CMFZ"))
                if cmfx:
                    self.comf = float(cmfx)
                if cmfz:
                    self.somf = float(cmfz)
            for idtag, store in (("PRID", self.points), ("CRID", self.curves),
                                 ("CCID", self.composites), ("SRID", self.surfaces),
                                 ("MRID", self.multipoints)):
                o = rec.first(idtag)
                if o is not None:
                    store[(_i(o.get("RCNM")), _i(o.get("RCID")))] = rec
                    break

    # -- coordinate helpers ------------------------------------------------
    def _xy(self, occ):
        x, y = _i(occ.get("XCOO")), _i(occ.get("YCOO"))
        if x is None or y is None:
            return None
        return (x / self.comf, y / self.comf)

    def _point_xy(self, rrid):
        rec = self.points.get((RCNM_PRID, _i(rrid)))
        if not rec:
            return None
        c = rec.first("C2IT")
        return self._xy(c) if c else None

    def _multipoint_xyz(self, rrid):
        """[(lon, lat, depth)] for an MRID multipoint (3-D soundings) from its
        C3IL tuples; depth = ZCOO / SOMF (CMFZ). Returns [] if unresolved."""
        rec = self.multipoints.get((RCNM_MRID, _i(rrid)))
        if not rec:
            return []
        out = []
        for t in rec.fields.get("C3IL", []):
            x, y, z = _i(t.get("XCOO")), _i(t.get("YCOO")), _i(t.get("ZCOO"))
            if x is None or y is None:
                continue
            depth = (z / self.somf) if z is not None else None
            out.append((x / self.comf, y / self.comf, depth))
        return out

    def _curve_coords(self, rrid):
        """Ordered coords of a CRID: begin point, interior C2IL, end point."""
        rec = self.curves.get((RCNM_CRID, _i(rrid)))
        if not rec:
            return None
        begin = end = None
        for o in rec.fields.get("PTAS", []):
            topi = _i(o.get("TOPI"))
            if topi == TOPI_BEGIN:
                begin = o.get("RRID")
            elif topi == TOPI_END:
                end = o.get("RRID")
        coords = []
        if begin is not None:
            b = self._point_xy(begin)
            if b:
                coords.append(b)
        for c in rec.fields.get("C2IL", []):
            xy = self._xy(c)
            if xy:
                coords.append(xy)
        if end is not None:
            e = self._point_xy(end)
            if e:
                coords.append(e)
        return _dedup(coords)

    def _composite_coords(self, rrid):
        """Ordered coords of a CCID: its CUCO component curves, each honouring
        ORNT, with the shared junction vertex de-duplicated."""
        rec = self.composites.get((RCNM_CCID, _i(rrid)))
        if not rec:
            return None
        coords = []
        for cu in rec.fields.get("CUCO", []):
            cc = self._curve_coords(cu.get("RRID"))
            if not cc:
                continue
            if _i(cu.get("ORNT")) == ORNT_REVERSE:
                cc = list(reversed(cc))
            if coords and cc and coords[-1] == cc[0]:
                cc = cc[1:]
            coords.extend(cc)
        return _dedup(coords)

    def _ref_coords(self, rrnm, rrid, ornt=None):
        """Coords for a curve-or-composite reference, honouring ORNT."""
        rrnm = _i(rrnm)
        if rrnm == RCNM_CCID:
            coords = self._composite_coords(rrid)
        elif rrnm == RCNM_CRID:
            coords = self._curve_coords(rrid)
        else:
            return None
        if coords and _i(ornt) == ORNT_REVERSE:
            coords = list(reversed(coords))
        return coords

    def _surface_rings(self, rrid):
        """(exterior_ring, [interior_rings]) for a SRID, or None."""
        rec = self.surfaces.get((RCNM_SRID, _i(rrid)))
        if not rec:
            return None
        exterior, interiors = None, []
        for ri in rec.fields.get("RIAS", []):
            ring = self._ref_coords(ri.get("RRNM"), ri.get("RRID"), ri.get("ORNT"))
            if not ring:
                continue
            if _i(ri.get("USAG")) == USAG_INTERIOR:
                interiors.append(ring)
            elif exterior is None:
                exterior = ring
            else:
                interiors.append(ring)   # extra exterior -> treat as separate part
        if exterior is None:
            return None
        return exterior, interiors

    # -- public ------------------------------------------------------------
    def geometry(self, feature):
        """Assemble and return a geometry value for an ``S101Feature`` from its
        spatial references, or ``None`` if it has none / they don't resolve.

        Dispatches on the referenced record type: points → Point/MultiPoint,
        curves/composites → LineString/MultiLineString, surfaces →
        Polygon/MultiPolygon.
        """
        refs = [occ for tag, occ in feature.spatial_refs if tag in ("SPAS", "PTAS")]
        if not refs:
            return None

        pts, mpts, lines, polys = [], [], [], []
        for occ in refs:
            rrnm = _i(occ.get("RRNM"))
            rrid = occ.get("RRID")
            if rrnm == RCNM_PRID:
                xy = self._point_xy(rrid)
                if xy:
                    pts.append(xy)
            elif rrnm == RCNM_MRID:
                mpts.extend(self._multipoint_xyz(rrid))
            elif rrnm in (RCNM_CRID, RCNM_CCID):
                coords = self._ref_coords(rrnm, rrid, occ.get("ORNT"))
                if coords and len(coords) >= 2:
                    lines.append(coords)
            elif rrnm == RCNM_SRID:
                rings = self._surface_rings(rrid)
                if rings:
                    polys.append([rings[0]] + rings[1])

        if polys:
            if len(polys) == 1:
                return Polygon(polys[0])
            return MultiPolygon(polys)
        if lines:
            if len(lines) == 1:
                return LineString(lines[0])
            return MultiLineString(lines)
        if mpts:
            return MultiPoint(mpts, z_present=True)
        if pts:
            if len(pts) == 1:
                return Point(pts[0][0], pts[0][1])
            return MultiPoint(pts)
        return None
