# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
validation/geometry.py
----------------------
S-58 spatial operators over OGR geometries.

S-58 1.4.2 says "for all spatial operators a default tolerance of 1/COMF should
be applied". COMF is the coordinate multiplication factor from the cell's DSPM
record (typically 1e7), so the tolerance is normally 1e-7 degrees. We snap that
tolerance into the predicates by buffering by a tiny epsilon where the standard
calls for tolerant comparison.

The five mutually exclusive DE-9IM relations (WITHIN, CROSSES, TOUCHES,
DISJOINT, OVERLAPS) and the helpers (INTERSECTS, CONTAINS, EQUALS) map straight
onto GEOS via OGR. COVERED_BY and COINCIDENT are NOT standard ISO 19125-1
operators -- S-58 Section 2.4 defines them specially -- so they are implemented
here per those definitions.

This module imports osgeo lazily inside the constructor, so importing the engine
without GDAL does not fail; only constructing GeomOps requires GDAL.
"""


class GeomOps:
    """Tolerance-aware spatial predicates. One instance is shared by the run."""

    def __init__(self, tolerance=1e-7):
        # Imported here so the module is importable without GDAL present.
        from osgeo import ogr  # noqa: F401
        self._ogr = ogr
        self.tol = tolerance

    # -- robustness helpers ---------------------------------------------------
    def _repair(self, g):
        """Return a valid version of a possibly-invalid geometry, or None.
        Unclosed rings / self-intersections make GEOS predicates raise;
        MakeValid (GDAL 3+) fixes most, Buffer(0) is the fallback."""
        if g is None:
            return None
        try:
            if hasattr(g, "MakeValid"):
                fixed = g.MakeValid()
                if fixed is not None and not fixed.IsEmpty():
                    return fixed
        except Exception:                       # noqa: BLE001
            pass
        try:
            return g.Buffer(0)
        except Exception:                       # noqa: BLE001
            return None

    def _rel(self, a, b, op, default=False):
        """Evaluate a GEOS predicate op(a, b) defensively. If invalid geometry
        makes GEOS raise, repair both operands and retry once; fall back to
        ``default`` if it still cannot be evaluated, so one bad geometry never
        aborts a whole check."""
        if not a or not b:
            return default
        try:
            return bool(op(a, b))
        except Exception:                       # noqa: BLE001
            a2, b2 = self._repair(a), self._repair(b)
            if a2 is not None and b2 is not None:
                try:
                    return bool(op(a2, b2))
                except Exception:               # noqa: BLE001
                    pass
            return default

    # -- standard DE-9IM relations (delegated to GEOS) ------------------------
    def equals(self, a, b):
        return self._rel(a, b, lambda x, y: x.Equals(y))

    def disjoint(self, a, b):
        return self._rel(a, b, lambda x, y: x.Disjoint(y), default=True)

    def intersects(self, a, b):
        return self._rel(a, b, lambda x, y: x.Intersects(y))

    def touches(self, a, b):
        return self._rel(a, b, lambda x, y: x.Touches(y))

    def crosses(self, a, b):
        return self._rel(a, b, lambda x, y: x.Crosses(y))

    def overlaps(self, a, b):
        return self._rel(a, b, lambda x, y: x.Overlaps(y))

    def within(self, a, b):
        return self._rel(a, b, lambda x, y: x.Within(y))

    def contains(self, a, b):
        return self._rel(a, b, lambda x, y: x.Contains(y))

    # -- S-58 special operators (Section 2.4) ---------------------------------
    def covered_by(self, a, b):
        """a is COVERED_BY b  <=>  (a INTERSECTION b) == a, i.e. no point of a
        lies outside b. Implemented tolerantly: a.Difference(b) is empty within
        tolerance. This is the operator most relational ENC checks actually
        want ("object sits on land/in a depth area") -- WITHIN alone wrongly
        excludes objects that touch the boundary. Invalid geometry is repaired
        and retried rather than allowed to abort the check.
        """
        if not a or not b:
            return False
        try:
            return self._covered_by_raw(a, b)
        except Exception:                       # noqa: BLE001
            a2, b2 = self._repair(a), self._repair(b)
            if a2 is None or b2 is None:
                return False
            try:
                return self._covered_by_raw(a2, b2)
            except Exception:                   # noqa: BLE001
                return False

    def _covered_by_raw(self, a, b):
        if a.Within(b) or a.Equals(b):
            return True
        # Tolerant containment: buffer b by tol so boundary-coincident parts of
        # a are absorbed, then test that nothing of a is left outside.
        b_pad = b.Buffer(self.tol) if self.tol else b
        diff = a.Difference(b_pad)
        if diff is None:
            return False
        return diff.IsEmpty() or diff.Area() <= self._area_eps(a)

    def coincident(self, a, b):
        """Line/line (or boundary-as-line) COINCIDENCE: the two share a linear
        portion -- they OVERLAP (as lines) or one is WITHIN the other. Polygon
        inputs are reduced to their boundary rings first, per S-58 2.4 ("for
        this operation the boundary of a Polygon ... is treated as Lines").

        Returns True when the linear intersection of a and b has length > tol.
        """
        la = self._as_lines(a)
        lb = self._as_lines(b)
        if la is None or lb is None:
            return False
        try:
            inter = la.Intersection(lb)
        except Exception:                       # noqa: BLE001
            la2, lb2 = self._repair(la), self._repair(lb)
            if la2 is None or lb2 is None:
                return False
            try:
                inter = la2.Intersection(lb2)
            except Exception:                   # noqa: BLE001
                return False
        if inter is None or inter.IsEmpty():
            return False
        return self._max_length(inter) > self.tol

    # -- helpers --------------------------------------------------------------
    def _as_lines(self, geom):
        """Return a linestring/multilinestring view of a geometry: the geometry
        itself if it is already linear, or its boundary if it is a polygon."""
        if geom is None:
            return None
        ogr = self._ogr
        gt = ogr.GT_Flatten(geom.GetGeometryType())
        if gt in (ogr.wkbLineString, ogr.wkbMultiLineString):
            return geom
        if gt in (ogr.wkbPolygon, ogr.wkbMultiPolygon):
            try:
                return geom.GetBoundary()
            except Exception:                   # noqa: BLE001
                g2 = self._repair(geom)
                try:
                    return g2.GetBoundary() if g2 is not None else None
                except Exception:               # noqa: BLE001
                    return None
        return None

    @staticmethod
    def _area_eps(geom):
        """A negligible-area threshold scaled to the feature, so sliver
        differences from floating-point boundary math don't read as 'outside'."""
        try:
            return max(geom.Area() * 1e-9, 1e-18)
        except Exception:       # noqa: BLE001
            return 1e-18

    def _max_length(self, geom):
        try:
            return geom.Length()
        except Exception:       # noqa: BLE001
            return 0.0
