# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
enc_coverage.py
---------------
ENC data-coverage geometry (M_COVR, CATCOV=1) and the cross-producer overlap
conflict detector.

Design (locked with Fabian 2026-07-02): the true coverage of a cell is its
M_COVR **CATCOV=1** polygon — NOT the rectangular file extent (rectangles overlap
freely and mean nothing). Per IHO / the NL+BE joint Westerschelde practice,
CATCOV=1 must not overlap CATCOV=1; overlap is only legal in CATCOV=2 (no-data)
zones, and along the jagged UNCLOS border the two producers' CATCOV=1 merely
TOUCH (overlapping rectangles + complementary CATCOV=1/2).

So a *conflict* between two producers is a CATCOV=1 overlap that survives eroding
each coverage inward by a tolerance (default 100 m) — this absorbs "just touching"
and the border slivers, and only flags genuine deep double-coverage. The erosion
+ intersection are done in a LOCAL METRIC CRS (auto UTM by centroid) so the
tolerance is true ground metres (EPSG:3857 would inflate ~1.6x at 51 N).

QGIS-FREE: the geometry MATH uses osgeo ogr/osr (GEOS) on WKT the caller has
already read with the in-house ISO 8211 engine — this module never reads an ENC
file itself, so both the dock (P1 scan-time alert) and the QGIS-free tile bake
(P2 CATCOV=1 masking) can use it. Input WKT is expected in EPSG:4326 (lon/lat),
as the native class layers and the decode buckets produce.
"""

CONFLICT_BUFFER_M = 100.0        # erosion tolerance for the overlap test (metres)
_MIN_CONFLICT_AREA_M2 = 1.0      # ignore degenerate slivers after erosion


def _ogr_osr():
    from osgeo import ogr, osr
    return ogr, osr


def _srs(osr, epsg):
    s = osr.SpatialReference()
    s.ImportFromEPSG(epsg)
    try:                                                # GDAL 3: force X=lon,Y=lat
        s.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    except Exception:                                   # noqa: BLE001
        pass
    return s


def _utm_epsg(lon, lat):
    """EPSG code of the UTM zone containing (lon, lat)."""
    zone = int((lon + 180.0) // 6) % 60 + 1
    return (32600 if lat >= 0 else 32700) + zone


def _geom(ogr, wkt):
    try:
        g = ogr.CreateGeometryFromWkt(wkt) if wkt else None
    except Exception:                                   # noqa: BLE001
        return None
    if g is not None and not g.IsValid():
        try:
            g = g.Buffer(0)                             # heal self-touching rings
        except Exception:                               # noqa: BLE001
            pass
    return g


def catcov1_union_wkt(rows):
    """Union of the CATCOV=1 area geometries in ``rows`` = iterable of
    (wkt, catcov) → a single (Multi)Polygon WKT in EPSG:4326, or None. A missing
    or blank CATCOV counts as 1 (coverage available), matching S-57 practice and
    the loader's coverage builder."""
    ogr, _osr = _ogr_osr()
    union = None
    for wkt, catcov in rows:
        try:
            if catcov not in (None, "") and int(float(str(catcov))) != 1:
                continue
        except (TypeError, ValueError):
            pass                                        # unparseable → treat as 1
        g = _geom(ogr, wkt)
        if g is None or g.IsEmpty():
            continue
        try:
            union = g.Clone() if union is None else union.Union(g)
        except Exception:                               # noqa: BLE001
            continue
    return union.ExportToWkt() if union is not None else None


def union_wkt(wkts):
    """Union of several 4326 WKT geometries → one WKT, or None."""
    ogr, _osr = _ogr_osr()
    u = None
    for w in wkts:
        g = _geom(ogr, w)
        if g is None or g.IsEmpty():
            continue
        try:
            u = g.Clone() if u is None else u.Union(g)
        except Exception:                               # noqa: BLE001
            continue
    return u.ExportToWkt() if u is not None else None


def difference_wkt(base_wkt, subtract_wkt):
    """``base`` minus ``subtract`` → WKT, or None if the result is empty. A
    None/blank ``subtract`` returns ``base`` unchanged (nothing to remove)."""
    ogr, _osr = _ogr_osr()
    b = _geom(ogr, base_wkt)
    if b is None or b.IsEmpty():
        return None
    if not subtract_wkt:
        return base_wkt
    s = _geom(ogr, subtract_wkt)
    if s is None or s.IsEmpty():
        return base_wkt
    try:
        d = b.Difference(s)
    except Exception:                                   # noqa: BLE001
        return base_wkt
    return d.ExportToWkt() if (d is not None and not d.IsEmpty()) else None


def clip_rows(rows, mask_wkt, kind):
    """Clip decode rows ``[(rcid, wkt, named), …]`` to ``mask_wkt`` (a coverage
    polygon in EPSG:4326). ``kind`` is the geometry kind ('point' | 'line' |
    'area'): points are kept when they fall inside the mask; lines/areas are
    intersected with it (empties dropped, the geometry replaced by the clipped
    WKT). A None/blank mask leaves the rows unchanged. Returns a NEW list.

    A bbox reject and a Contains fast-path keep this cheap: features fully inside
    the mask are kept untouched (no expensive intersection); features fully
    outside are dropped without one."""
    if not mask_wkt:
        return rows
    ogr, _osr = _ogr_osr()
    mask = _geom(ogr, mask_wkt)
    if mask is None or mask.IsEmpty():
        return rows
    me = mask.GetEnvelope()                             # (minx,maxx,miny,maxy)
    is_point = (kind == "point")
    out = []
    for rcid, wkt, named in rows:
        g = _geom(ogr, wkt)
        if g is None or g.IsEmpty():
            continue
        ge = g.GetEnvelope()
        if ge[1] < me[0] or me[1] < ge[0] or ge[3] < me[2] or me[3] < ge[2]:
            continue                                    # bbox-disjoint → drop
        if is_point:
            if mask.Intersects(g):
                out.append((rcid, wkt, named))
            continue
        try:
            if mask.Contains(g):                        # fully inside → keep as-is
                out.append((rcid, wkt, named))
                continue
        except Exception:                               # noqa: BLE001
            pass
        try:
            clipped = g.Intersection(mask)
        except Exception:                               # noqa: BLE001
            clipped = None
        if clipped is not None and not clipped.IsEmpty():
            out.append((rcid, clipped.ExportToWkt(), named))
    return out


def clip_points_xyz(points, mask_wkt):
    """Keep the (x, y, z) tuples whose (x, y) falls inside ``mask_wkt`` (a 4326
    coverage polygon). A None/blank mask leaves them unchanged. Returns a NEW
    list. Used for the bake's SOUNDG list (soundings are (x, y, z))."""
    if not mask_wkt:
        return list(points)
    ogr, _osr = _ogr_osr()
    mask = _geom(ogr, mask_wkt)
    if mask is None or mask.IsEmpty():
        return list(points)
    me = mask.GetEnvelope()
    out = []
    for xyz in points:
        x, y = xyz[0], xyz[1]
        if x < me[0] or x > me[1] or y < me[2] or y > me[3]:
            continue
        pt = ogr.Geometry(ogr.wkbPoint)
        pt.AddPoint_2D(float(x), float(y))
        if mask.Intersects(pt):
            out.append(xyz)
    return out


def _to_metric(ogr, osr, wkt, epsg):
    """Reproject a 4326 WKT geometry into the given (metric) EPSG."""
    g = _geom(ogr, wkt)
    if g is None or g.IsEmpty():
        return None
    ct = osr.CoordinateTransformation(_srs(osr, 4326), _srs(osr, epsg))
    if g.Transform(ct) != 0:
        return None
    return g


def coverages_conflict(wkt_a, wkt_b, buffer_m=CONFLICT_BUFFER_M):
    """True when two CATCOV=1 coverages overlap by MORE than ``buffer_m`` metres:
    each is eroded inward by buffer_m in a local UTM CRS, then intersected. A
    touch, or a sliver overlap from the jagged UNCLOS border, erodes away → not a
    conflict; only genuine deep double-coverage survives. BOTH geometries are put
    in the SAME UTM zone (from A's centroid) so a zone boundary can't corrupt the
    intersection."""
    ogr, osr = _ogr_osr()
    try:
        ga0 = _geom(ogr, wkt_a)
        if ga0 is None or ga0.IsEmpty():
            return False
        c = ga0.Centroid()
        epsg = _utm_epsg(c.GetX(), c.GetY())
        ga = _to_metric(ogr, osr, wkt_a, epsg)
        gb = _to_metric(ogr, osr, wkt_b, epsg)
        if ga is None or gb is None or not ga.Intersects(gb):
            return False                                # cheap reject: no overlap
        ea, eb = ga.Buffer(-buffer_m), gb.Buffer(-buffer_m)
        if ea is None or eb is None or ea.IsEmpty() or eb.IsEmpty():
            return False                                # overlap within tolerance
        inter = ea.Intersection(eb)
        return (inter is not None and not inter.IsEmpty()
                and inter.GetArea() > _MIN_CONFLICT_AREA_M2)
    except Exception:                                   # noqa: BLE001
        return False


def _envelope(wkt):
    """(minx, maxx, miny, maxy) of a WKT geometry, or None."""
    ogr, _osr = _ogr_osr()
    g = _geom(ogr, wkt)
    return g.GetEnvelope() if g is not None and not g.IsEmpty() else None


def _bbox_overlap(ea, eb, pad_deg):
    """Do two (minx,maxx,miny,maxy) envelopes overlap within a padding (degrees)?"""
    return not (ea[1] + pad_deg < eb[0] or eb[1] + pad_deg < ea[0]
                or ea[3] + pad_deg < eb[2] or eb[3] + pad_deg < ea[2])


def find_conflicts(coverages, buffer_m=CONFLICT_BUFFER_M):
    """Detect cross-producer coverage conflicts.

    ``coverages`` = [{"cell","usage","agen","wkt"}]. Returns a list of
    (cellA, cellB, usage, agenA, agenB) for pairs of the SAME usage band but
    DIFFERENT producing agency whose CATCOV=1 overlap beyond ``buffer_m``.
    Same-agency pairs are skipped — staging already dedups those by cell name /
    edition (Rule 2), and a well-authored joint production (shared agency) is
    authored to TOUCH, not overlap. A cheap bounding-box pre-filter skips the
    (expensive) metric erode+intersect for pairs that can't possibly overlap."""
    covs = [dict(c, _env=_envelope(c["wkt"]))
            for c in coverages if c.get("wkt")]
    covs = [c for c in covs if c["_env"] is not None]
    pad = buffer_m / 111320.0 * 2.0          # metres → deg latitude, generous
    out = []
    for i in range(len(covs)):
        for j in range(i + 1, len(covs)):
            a, b = covs[i], covs[j]
            if a.get("usage") != b.get("usage"):
                continue
            if (a.get("agen") or "").strip() == (b.get("agen") or "").strip():
                continue
            if not _bbox_overlap(a["_env"], b["_env"], pad):
                continue                     # bboxes disjoint → cannot conflict
            if coverages_conflict(a["wkt"], b["wkt"], buffer_m):
                out.append((a.get("cell"), b.get("cell"), a.get("usage"),
                            (a.get("agen") or "?"), (b.get("agen") or "?")))
    return out
