# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
s57/geomops.py
--------------
Spatial **predicates and measures** over native geometry values — the layer the
``s57.geometry`` docstring reserves as ``GeomOps``.

The engine *assembles* geometry values (``Point``/``MultiPoint``/``LineString``/
``MultiLineString``/``Polygon``/``MultiPolygon`` in ``s57.geometry``). Anything
that asks a *question* about those values — how big, how long, does a point fall
inside, do two shapes overlap — belongs here, behind one seam, so consumers ask
``GeomOps`` rather than reaching for a geometry library directly.

Two backends, by design:

* **Pure-Python core (default, GDAL-free).** Measures (bbox, area, length,
  perimeter, centroid) and the predicates that are robust without a full
  geometry engine (point-in-polygon, bbox intersection, ring closure). This
  keeps the engine importable and these operations usable with no osgeo present.

* **GEOS backend (optional, lazy).** The genuinely set-theoretic predicates
  (``intersects``/``overlaps``/``contains``/``touches``/``is_valid``/
  ``is_simple``) are delegated to GEOS via ``osgeo`` — fed the value's own WKT,
  ``ogr.CreateGeometryFromWkt(geom.wkt())``. So when GDAL is present it is only a
  *calculator*: the native engine remains the reader, and the verified native
  assembly is what GEOS operates on. ``has_geos()`` reports availability; the
  set predicates raise ``GeomOpsUnavailable`` if asked without it.

Coordinate frame: GeomOps is frame-agnostic and works in whatever units the
values carry (the engine emits lon/lat degrees). Planar formulae are used for
area/length — adequate for the cell-local, predicate-style questions the checks
ask; a geodesic mode can be added behind the same seam if a measure ever needs
true metres.
"""

import math


class GeomOpsUnavailable(RuntimeError):
    """Raised when a set-theoretic predicate is requested but no GEOS backend is
    available (osgeo not importable in this environment)."""


# ---------------------------------------------------------------------------
# Coordinate-level helpers (pure)
# ---------------------------------------------------------------------------
def _rings_of(geom):
    """Yield each ring (list of (x, y)) of a Polygon/MultiPolygon."""
    if getattr(geom, "kind", None) == "POLYGON":
        for r in geom.rings:
            yield r
    elif getattr(geom, "kind", None) == "MULTIPOLYGON":
        for poly in geom.polygons:
            for r in poly:
                yield r


def _polys_of(geom):
    """Yield each polygon as ``[exterior, *holes]`` for Polygon/MultiPolygon."""
    if getattr(geom, "kind", None) == "POLYGON":
        yield geom.rings
    elif getattr(geom, "kind", None) == "MULTIPOLYGON":
        for poly in geom.polygons:
            yield poly


def _lines_of(geom):
    """Yield each line part (list of (x, y)) of a LineString/MultiLineString."""
    if getattr(geom, "kind", None) == "LINESTRING":
        yield geom.points
    elif getattr(geom, "kind", None) == "MULTILINESTRING":
        for p in geom.parts:
            yield p


def _ring_area(ring):
    """Signed shoelace area of a ring (positive = counter-clockwise)."""
    n = len(ring)
    if n < 3:
        return 0.0
    s = 0.0
    for i in range(n):
        x1, y1 = ring[i][0], ring[i][1]
        x2, y2 = ring[(i + 1) % n][0], ring[(i + 1) % n][1]
        s += x1 * y2 - x2 * y1
    return s / 2.0


def _seg_len(a, b):
    return math.hypot(b[0] - a[0], b[1] - a[1])


def _point_in_ring(x, y, ring):
    """Ray-cast: is (x, y) inside this ring? Boundary result is not guaranteed
    (callers treat on-boundary as a separate concern)."""
    inside = False
    n = len(ring)
    if n < 3:
        return False
    j = n - 1
    for i in range(n):
        xi, yi = ring[i][0], ring[i][1]
        xj, yj = ring[j][0], ring[j][1]
        if ((yi > y) != (yj > y)) and \
           (x < (xj - xi) * (y - yi) / ((yj - yi) or 1e-30) + xi):
            inside = not inside
        j = i
    return inside


# ---------------------------------------------------------------------------
# GeomOps
# ---------------------------------------------------------------------------
class GeomOps:
    """Predicate/measure operations over native geometry values.

    Construct once and reuse; the GEOS handle (if any) is resolved lazily on
    first set-predicate use. Measures and point/bbox predicates never need GEOS.
    """

    def __init__(self, prefer_geos=True):
        self._prefer_geos = prefer_geos
        self._ogr = None
        self._geos_resolved = False

    # -- measures (pure) ------------------------------------------------------
    @staticmethod
    def bbox(geom):
        """(minx, miny, maxx, maxy) over all coordinates, or None if empty."""
        if geom is None:
            return None
        cs = geom.coords()
        if not cs:
            return None
        xs = [c[0] for c in cs]
        ys = [c[1] for c in cs]
        return (min(xs), min(ys), max(xs), max(ys))

    @staticmethod
    def area(geom):
        """Absolute planar area of a Polygon/MultiPolygon (exterior minus holes);
        0.0 for non-areal geometry. Units are the square of the value's units."""
        total = 0.0
        for poly in _polys_of(geom):
            if not poly:
                continue
            ext = abs(_ring_area(poly[0]))
            holes = sum(abs(_ring_area(h)) for h in poly[1:])
            total += max(0.0, ext - holes)
        return total

    @staticmethod
    def length(geom):
        """Planar length of a LineString/MultiLineString, or the perimeter of a
        Polygon/MultiPolygon (all rings). 0.0 otherwise."""
        total = 0.0
        for part in _lines_of(geom):
            for i in range(len(part) - 1):
                total += _seg_len(part[i], part[i + 1])
        for ring in _rings_of(geom):
            for i in range(len(ring) - 1):
                total += _seg_len(ring[i], ring[i + 1])
        return total

    @staticmethod
    def centroid(geom):
        """Area-weighted centroid for areal geometry, else the mean of vertices.
        Returns (x, y) or None for empty geometry."""
        polys = list(_polys_of(geom))
        if polys:
            ax = ay = asum = 0.0
            for poly in polys:
                if not poly:
                    continue
                ring = poly[0]
                a = _ring_area(ring)
                if a == 0.0:
                    continue
                cx = cy = 0.0
                n = len(ring)
                for i in range(n):
                    x1, y1 = ring[i][0], ring[i][1]
                    x2, y2 = ring[(i + 1) % n][0], ring[(i + 1) % n][1]
                    cross = x1 * y2 - x2 * y1
                    cx += (x1 + x2) * cross
                    cy += (y1 + y2) * cross
                cx /= (6 * a)
                cy /= (6 * a)
                ax += cx * a
                ay += cy * a
                asum += a
            if asum:
                return (ax / asum, ay / asum)
        cs = geom.coords() if geom is not None else []
        if not cs:
            return None
        return (sum(c[0] for c in cs) / len(cs),
                sum(c[1] for c in cs) / len(cs))

    # -- point / bbox predicates (pure) ---------------------------------------
    @staticmethod
    def point_in_polygon(x, y, geom):
        """True if (x, y) is inside a Polygon/MultiPolygon: within an exterior
        ring and not within any of that polygon's holes."""
        for poly in _polys_of(geom):
            if not poly:
                continue
            if _point_in_ring(x, y, poly[0]) and \
               not any(_point_in_ring(x, y, h) for h in poly[1:]):
                return True
        return False

    @classmethod
    def bbox_intersects(cls, a, b):
        """Cheap rejection test: do the bounding boxes overlap? Useful as a
        pre-filter before an exact (GEOS) predicate."""
        ba, bb = cls.bbox(a), cls.bbox(b)
        if ba is None or bb is None:
            return False
        return not (ba[2] < bb[0] or bb[2] < ba[0] or
                    ba[3] < bb[1] or bb[3] < ba[1])

    @staticmethod
    def ring_is_closed(ring):
        """A ring is closed when its first and last vertex coincide."""
        return len(ring) >= 4 and ring[0] == ring[-1]

    # -- set predicates (GEOS-backed) -----------------------------------------
    def has_geos(self):
        """Whether a GEOS backend (osgeo) is available in this environment."""
        return self._geos() is not None

    def _geos(self):
        if not self._geos_resolved:
            self._geos_resolved = True
            if self._prefer_geos:
                try:
                    from osgeo import ogr
                    self._ogr = ogr
                except Exception:               # noqa: BLE001
                    self._ogr = None
        return self._ogr

    def _to_geos(self, geom):
        ogr = self._geos()
        if ogr is None or geom is None:
            return None
        try:
            return ogr.CreateGeometryFromWkt(geom.wkt())
        except Exception:                       # noqa: BLE001
            return None

    def _predicate(self, name, a, b):
        ogr = self._geos()
        if ogr is None:
            raise GeomOpsUnavailable(
                f"{name} needs a GEOS backend (osgeo not available)")
        # bbox pre-filter for the binary predicates that imply spatial overlap
        if name in ("intersects", "overlaps", "contains", "within") \
                and not self.bbox_intersects(a, b):
            return False
        ga, gb = self._to_geos(a), self._to_geos(b)
        if ga is None or gb is None:
            raise GeomOpsUnavailable(f"{name}: geometry not representable")
        return bool(getattr(ga, name.capitalize())(gb))

    def intersects(self, a, b):
        """True if a and b share any point (GEOS)."""
        return self._predicate("intersects", a, b)

    def overlaps(self, a, b):
        """True if a and b overlap in 2-D interior (GEOS)."""
        return self._predicate("overlaps", a, b)

    def contains(self, a, b):
        """True if a contains b (GEOS)."""
        return self._predicate("contains", a, b)

    def touches(self, a, b):
        """True if a and b touch only at boundaries (GEOS)."""
        return self._predicate("touches", a, b)

    def is_valid(self, geom):
        """OGC validity of a single geometry (GEOS)."""
        ogr = self._geos()
        if ogr is None:
            raise GeomOpsUnavailable("is_valid needs a GEOS backend")
        g = self._to_geos(geom)
        if g is None:
            return False
        try:
            return bool(g.IsValid())
        except Exception:                       # noqa: BLE001
            return False

    def is_simple(self, geom):
        """Whether a line geometry is simple (no self-intersection) (GEOS)."""
        ogr = self._geos()
        if ogr is None:
            raise GeomOpsUnavailable("is_simple needs a GEOS backend")
        g = self._to_geos(geom)
        if g is None:
            return False
        try:
            return bool(g.IsSimple())
        except Exception:                       # noqa: BLE001
            return False
