# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
validation/s101_geom.py
-----------------------
The **spatial layer** for the S-158-101 geometry-topology checks (the ``101_00xx``
band).

The mechanical S-101 checks in :mod:`validation.s101` read only the native
feature/attribute model. The topology band, though, is about *geometry* — holes
and overlaps in the Skin of the Earth, coincident curves, a point feature sitting
on a same-class surface, a cell crossing the antimeridian. To evaluate those we
need assembled geometry and set-theoretic spatial predicates.

This module is the S-101 sibling of :mod:`validation.geometry` (the S-58 spatial
operators). For each feature it assembles geometry with the **native** builder
(:class:`s101.geometry.S101GeometryBuilder` over the cell's retained ISO 8211
records — GDAL-free assembly, exactly as the on-canvas loader does), turns that
native WKT into an OGR geometry, and pairs it with a shared
:class:`validation.geometry.GeomOps` for the tolerant DE-9IM / COINCIDENT /
COVERED_BY predicates.

GDAL is used only as the *predicate calculator* (GEOS via OGR) over geometry that
was **read natively** — the standing "S-101 is never read with GDAL" rule governs
the ISO 8211 decode, which stays native throughout. If GDAL is unavailable the
spatial layer simply does not build and the geometry checks no-op (the mechanical
checks are unaffected), so the engine still imports and runs headlessly.
"""

try:                                                    # plugin context
    from ..s101.geometry import (S101GeometryBuilder, _i, RCNM_SRID, RCNM_CRID,
                                 RCNM_CCID, ORNT_REVERSE)
except (ImportError, ValueError):                       # headless / top-level
    from s101.geometry import (S101GeometryBuilder, _i, RCNM_SRID, RCNM_CRID,
                               RCNM_CCID, ORNT_REVERSE)

from .geometry import GeomOps


# S-101 spatial record names (RCNM), grouped by geometric primitive (S-100 Part 7).
_POINT_RCNM = (110, 115)        # PRID (point), MRID (multipoint)
_CURVE_RCNM = (120, 125)        # CRID (curve), CCID (composite curve)
_SURFACE_RCNM = (130,)          # SRID (surface)


def primitive_of_rrnm(rrnm):
    """Map a spatial-record name (RRNM) to its geometric primitive class
    (``point`` / ``curve`` / ``surface``), or ``None`` if unrecognised."""
    try:
        rrnm = int(rrnm)
    except (TypeError, ValueError):
        return None
    if rrnm in _POINT_RCNM:
        return "point"
    if rrnm in _CURVE_RCNM:
        return "curve"
    if rrnm in _SURFACE_RCNM:
        return "surface"
    return None


def feature_primitive_classes(feature):
    """The set of geometric-primitive classes a feature references through its
    raw spatial pointers (SPAS/PTAS RRNM). A well-formed feature references a
    single primitive; a mix is the 101_0016 violation. Reads only the native
    model, so it needs no geometry assembly (and no GDAL)."""
    classes = set()
    for tag, occ in getattr(feature, "spatial_refs", None) or []:
        if tag in ("SPAS", "PTAS"):
            c = primitive_of_rrnm(occ.get("RRNM"))
            if c:
                classes.add(c)
    return classes


def _iter_feature_wkt(cell):
    """Yield ``(feature, wkt)`` for a cell's assemblable features.

    Two sources, in priority order:
      * a pre-attached ``cell.feature_geoms`` = ``[(feature, Geometry|wkt), ...]``
        — the seam the headless tests drive (geometry doubles, no ISO 8211 bytes);
      * otherwise the cell's retained raw records (``cell._records``, set by
        :func:`validation.s101.read_cells`), assembled with the native builder.

    A cell exposing neither yields nothing, so the spatial layer is simply empty
    for feature/attribute-only doubles and the geometry checks no-op on them.
    """
    fg = getattr(cell, "feature_geoms", None)
    if fg is not None:
        for feature, geom in fg:
            wkt = geom.wkt() if hasattr(geom, "wkt") else geom
            if wkt:
                yield feature, wkt
        return
    recs = getattr(cell, "_records", None)
    if not recs:
        return
    builder = S101GeometryBuilder(recs)
    for feature in getattr(cell, "features", None) or []:
        g = builder.geometry(feature)
        if g is not None:
            yield feature, g.wkt()


class GeomEntry:
    """One assembled feature geometry: the owning cell, the native feature, its
    OGR geometry, primitive class, and cached (minX, maxX, minY, maxY) envelope."""
    __slots__ = ("cell", "feature", "geom", "primitive", "envelope")

    def __init__(self, cell, feature, geom, primitive, envelope):
        self.cell = cell
        self.feature = feature
        self.geom = geom
        self.primitive = primitive
        self.envelope = envelope


def _primitive_of_ogr(ogr, geom):
    gt = ogr.GT_Flatten(geom.GetGeometryType())
    if gt in (ogr.wkbPoint, ogr.wkbMultiPoint):
        return "point"
    if gt in (ogr.wkbLineString, ogr.wkbMultiLineString):
        return "curve"
    if gt in (ogr.wkbPolygon, ogr.wkbMultiPolygon):
        return "surface"
    return None


def envelopes_overlap(a, b):
    """True if two GetEnvelope() boxes (minX, maxX, minY, maxY) intersect — the
    cheap pre-filter before an expensive GEOS predicate."""
    return not (a[1] < b[0] or b[1] < a[0] or a[3] < b[2] or b[3] < a[2])


def union_geoms(geoms):
    """The set-theoretic union of OGR geometries (for coverage-equality tests), or
    ``None`` for an empty input. Defensive: a geometry that fails to union is
    skipped rather than aborting the merge."""
    if not geoms:
        return None
    u = geoms[0].Clone()
    for g in geoms[1:]:
        try:
            nxt = u.Union(g)
        except Exception:                   # noqa: BLE001
            continue
        if nxt is not None:
            u = nxt
    return u


class S101Spatial:
    """A run's assembled S-101 geometry plus the shared spatial operators. Built
    once per audit (cached on the run ctx) and queried by the topology checks.

    Construction is defensive: without GDAL, or if nothing assembles, ``build``
    returns ``None`` and every geometry check no-ops."""

    def __init__(self, ops):
        self.ops = ops
        self.entries = []                   # [GeomEntry]

    @classmethod
    def build(cls, cells):
        """Assemble every feature's geometry across ``cells`` into OGR geometries,
        or return ``None`` if GDAL is unavailable / nothing assembles."""
        try:
            from osgeo import ogr
        except Exception:                   # noqa: BLE001 — GDAL absent
            return None
        try:
            ops = GeomOps()
        except Exception:                   # noqa: BLE001
            return None
        self = cls(ops)
        for cell in cells or []:
            for feature, wkt in _iter_feature_wkt(cell):
                try:
                    geom = ogr.CreateGeometryFromWkt(wkt)
                except Exception:           # noqa: BLE001 — malformed WKT
                    geom = None
                if geom is None:
                    continue
                self.entries.append(GeomEntry(
                    cell, feature, geom,
                    _primitive_of_ogr(ogr, geom), geom.GetEnvelope()))
        return self if self.entries else None

    # -- queries --------------------------------------------------------------
    def cells_seen(self):
        """The cells that produced at least one geometry, in first-seen order.
        Topology checks scope PER CELL — adjacent cells legitimately abut/overlap
        at their shared borders, so a dataset-wide sweep would cry wolf."""
        seen, out = set(), []
        for e in self.entries:
            if id(e.cell) not in seen:
                seen.add(id(e.cell))
                out.append(e.cell)
        return out

    def entries_for(self, cell):
        return [e for e in self.entries if e.cell is cell]


# ---------------------------------------------------------------------------
# Curve-usage graph — the topology of shared boundary curves (for 101_0001 /
# 101_0004 / 101_0014). This is GDAL-FREE: it reads only the native records, so
# these checks run headlessly. It walks each surface feature's ring references
# (SPAS→SRID → RIAS → CRID, or → CCID → CUCO → CRID) down to the SIMPLE curves,
# composing orientation along the way, and records which feature used each curve.
# ---------------------------------------------------------------------------
class CurveGraph:
    """Per-cell curve topology: ``usage`` maps a simple-curve id (CRID rcid) to the
    list of ``(feature, net_reverse)`` that reference it (net_reverse = the curve's
    orientation within that feature's ring, composed through any CCID). ``curves``
    is the raw ``{(RCNM, rcid): record}`` CRID store (for the SEGH.INTP check)."""
    __slots__ = ("usage", "curves")

    def __init__(self, usage, curves):
        self.usage = usage
        self.curves = curves


def _surface_curve_usages(builder, srid_rrid):
    """Yield ``(curve_rcid, net_reverse)`` for every simple curve in a surface's
    rings, composing the RIAS orientation with the CCID's CUCO orientation
    (net reverse = exactly one of the two is reversed)."""
    rec = builder.surfaces.get((RCNM_SRID, _i(srid_rrid)))
    if not rec:
        return
    for ri in rec.fields.get("RIAS", []):
        rrnm, rrid = _i(ri.get("RRNM")), _i(ri.get("RRID"))
        rev = (_i(ri.get("ORNT")) == ORNT_REVERSE)
        if rrnm == RCNM_CRID:
            yield rrid, rev
        elif rrnm == RCNM_CCID:
            crec = builder.composites.get((RCNM_CCID, rrid))
            if not crec:
                continue
            for cu in crec.fields.get("CUCO", []):
                cu_rev = (_i(cu.get("ORNT")) == ORNT_REVERSE)
                yield _i(cu.get("RRID")), (rev ^ cu_rev)


def build_curve_graph(cell):
    """Assemble a cell's :class:`CurveGraph` from its retained native records
    (``cell._records``). A test double may pre-attach ``cell.curve_graph`` (used
    verbatim). Returns ``None`` when neither is available (feature/attribute-only
    doubles), so the curve-topology checks no-op on them."""
    pre = getattr(cell, "curve_graph", None)
    if pre is not None:
        return pre
    recs = getattr(cell, "_records", None)
    if not recs:
        return None
    builder = S101GeometryBuilder(recs)
    usage = {}
    for feature in getattr(cell, "features", None) or []:
        for tag, occ in getattr(feature, "spatial_refs", None) or []:
            if tag != "SPAS" or _i(occ.get("RRNM")) != RCNM_SRID:
                continue
            for curid, rev in _surface_curve_usages(builder, _i(occ.get("RRID"))):
                usage.setdefault(curid, []).append((feature, rev))
    return CurveGraph(usage, builder.curves)


def _curve_wkt(builder, rrnm, rrid):
    """LINESTRING WKT for a simple curve (CRID) or composite curve (CCID) from the
    native builder, or ``None`` if it doesn't resolve to >= 2 points."""
    rrnm = _i(rrnm)
    if rrnm == RCNM_CRID:
        coords = builder._curve_coords(rrid)
    elif rrnm == RCNM_CCID:
        coords = builder._composite_coords(rrid)
    else:
        return None
    if not coords or len(coords) < 2:
        return None
    return "LINESTRING (" + ",".join("%r %r" % (x, y) for x, y in coords) + ")"


def data_limit_mask_violations(spatial, cell):
    """Yield ``(feature, curve_rrid, mind)`` for MASK-flagged curves that run along
    the DataCoverage limit but whose masking indicator MIND is not 1 (they should
    be masked = truncated by the coverage limit). Conservative: only curves the
    encoder explicitly MASK-ed are examined, and MIND == 1 is skipped before any
    geometry work, so the cost is one coincidence test per non-masked flag.

    A test double may inject ``cell.mask_geoms`` = ``[(feature, rrid, mind, wkt)]``
    instead of supplying native records."""
    from osgeo import ogr
    ops = spatial.ops
    covs = [e.geom for e in spatial.entries_for(cell)
            if e.primitive == "surface"
            and getattr(e.feature, "type_name", None) == "DataCoverage"]
    if not covs:
        return
    border = None
    for g in covs:
        try:
            b = g.GetBoundary()
        except Exception:                   # noqa: BLE001
            b = None
        if b is not None:
            border = b if border is None else border.Union(b)
    if border is None:
        return

    def _hit(wkt, mind):
        if mind is None or mind == 1 or not wkt:
            return False
        try:
            line = ogr.CreateGeometryFromWkt(wkt)
        except Exception:                   # noqa: BLE001
            return False
        return line is not None and ops.coincident(line, border)

    injected = getattr(cell, "mask_geoms", None)
    if injected is not None:
        for feature, rrid, mind, wkt in injected:
            if _hit(wkt, mind):
                yield feature, rrid, mind
        return
    recs = getattr(cell, "_records", None)
    if not recs:
        return
    builder = S101GeometryBuilder(recs)
    for feature in getattr(cell, "features", None) or []:
        if getattr(feature, "type_name", None) == "DataCoverage":
            continue
        for tag, occ in getattr(feature, "spatial_refs", None) or []:
            if tag != "MASK":
                continue
            mind = _i(occ.get("MIND"))
            if mind is None or mind == 1:
                continue
            wkt = _curve_wkt(builder, occ.get("RRNM"), occ.get("RRID"))
            if _hit(wkt, mind):
                yield feature, _i(occ.get("RRID")), mind
