# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
sourceval/areamatch.py — area-geometry conflation with a generalization dead-band.
================================================================================

The AREA companion to geomatch (lines): compare a charted polygon (basin,
land area, restricted area) against its source counterpart on SHAPE, inside a
dead-band that stays silent for acceptable generalization, and flag only an
OBVIOUS disagreement.

The two traps a naive surface (m²/ha) comparison falls into, and the design
that dodges them:

  * CHART POLYGONS ARE CLIPPED AT THE CELL BORDER, source polygons are not —
    a raw area comparison false-flags every border feature. Therefore the
    source polygon is FIRST CLIPPED TO THE CHART'S COVERAGE (M_COVR): both
    sides then describe the same footprint and the border is a non-issue by
    construction.
  * EQUAL AREA CAN HIDE DIFFERENT SHAPES (a circle and a star). Therefore the
    verdict is based on the coverage-clipped SYMMETRIC DIFFERENCE (localizable
    disagreement lobes, reported in m²) and IoU — the scalar surface check is
    subsumed, not skipped.

Belt-and-braces: the symmetric difference is ERODED by half the tolerance
before judging, which removes the thin slivers produced by the coverage
polygon's own digitization wobble and by within-tolerance generalization.

Verdict (the dead-band): pass silently when IoU ≥ AREA_PASS_IOU and no eroded
lobe survives above the lobe floor; flag when IoU < AREA_FLAG_IOU or a lobe
exceeds the floor. Between the bands: silence (no nagging middle tier).

Polygon boolean operations use OGR/GEOS (available wherever the plugin runs;
the module degrades to no-op without it, mirroring the geometry checks).
Coordinates arrive as lon/lat rings and are projected to local metres first
(geomatch.project_local), so every figure below is metric.
"""

from . import geomatch

# ── the dead-band (tunable, mirrors geomatch/QC1 constant style) ─────────────
AREA_PASS_IOU = 0.90      # ≥ this intersection-over-union → silent
AREA_FLAG_IOU = 0.75      # < this → obvious disagreement
LOBE_FLOOR_M2 = 100.0     # an eroded disagreement lobe below this never flags
LOBE_FLOOR_X = 25.0       # … and the floor grows with tolerance: X × tol²
MATCH_MIN_IOU = 0.30      # pairing: best IoU must reach this to count as a pair


def _ogr():
    try:
        from osgeo import ogr
        ogr.UseExceptions()
        return ogr
    except Exception:                                   # noqa: BLE001
        return None


def has_geos():
    """True when polygon boolean operations are available."""
    o = _ogr()
    if o is None:
        return False
    try:
        g = o.CreateGeometryFromWkt("POLYGON ((0 0,1 0,1 1,0 1,0 0))")
        return g.Intersection(g) is not None
    except Exception:                                   # noqa: BLE001
        return False


def _poly_from_rings(ogr_mod, rings_xy):
    """[[(x, y), …] exterior, *holes] → OGR Polygon (closed rings enforced)."""
    poly = ogr_mod.Geometry(ogr_mod.wkbPolygon)
    for ring in rings_xy:
        if len(ring) < 3:
            continue
        r = ogr_mod.Geometry(ogr_mod.wkbLinearRing)
        for x, y in ring:
            r.AddPoint_2D(float(x), float(y))
        if ring[0] != ring[-1]:
            r.AddPoint_2D(float(ring[0][0]), float(ring[0][1]))
        poly.AddGeometry(r)
    return poly


def _multi_from_areas(ogr_mod, areas_xy):
    """[[rings], …] (several polygons) → one OGR (Multi)Polygon, made valid."""
    if len(areas_xy) == 1:
        g = _poly_from_rings(ogr_mod, areas_xy[0])
    else:
        g = ogr_mod.Geometry(ogr_mod.wkbMultiPolygon)
        for rings in areas_xy:
            g.AddGeometry(_poly_from_rings(ogr_mod, rings))
    try:                                    # heal self-touches/slivers cheaply
        if not g.IsValid():
            g = g.Buffer(0)
    except Exception:                                   # noqa: BLE001
        pass
    return g


def _project_areas(areas_ll, lon0, lat0):
    return [[geomatch.project_local(r, lon0, lat0) for r in rings]
            for rings in areas_ll]


def compare_areas(chart_areas_ll, source_areas_ll, tol_m, coverage_ll=None):
    """Compare charted polygons against source polygons (both
    ``[[exterior, *holes], …]`` in lon/lat) with tolerance ``tol_m``.

    ``coverage_ll`` (optional, [[exterior, *holes], …]) is the chart's real
    footprint (M_COVR): the SOURCE is clipped to it first, so a source polygon
    continuing beyond the cell border never reads as a disagreement.

    Returns a dict (or None when OGR/GEOS is unavailable):
      flag         True only for an OBVIOUS difference
      iou          intersection-over-union after coverage clipping
      symdiff_m2   surface disagreement (m², after erosion by tol/2)
      worst_lobe   (area_m2, (lon, lat)) of the largest disagreement lobe
      chart_m2     charted surface (m²)
    """
    o = _ogr()
    if o is None or not chart_areas_ll or not source_areas_ll:
        return None
    tol = max(0.01, float(tol_m))
    lon0, lat0 = chart_areas_ll[0][0][0]
    chart = _multi_from_areas(o, _project_areas(chart_areas_ll, lon0, lat0))
    source = _multi_from_areas(o, _project_areas(source_areas_ll, lon0, lat0))
    try:
        if coverage_ll:
            cov = _multi_from_areas(o, _project_areas(coverage_ll, lon0, lat0))
            source = source.Intersection(cov)
            if source is None or source.IsEmpty():
                return None                       # source lies outside the chart
        inter = chart.Intersection(source)
        union = chart.Union(source)
        ia = inter.GetArea() if inter is not None else 0.0
        ua = union.GetArea() if union is not None else 0.0
        iou = (ia / ua) if ua > 0 else 0.0
        sym = chart.SymDifference(source)
        # Erode: kills coverage-wobble slivers + within-tolerance fuzz.
        sym_er = sym.Buffer(-tol / 2.0) if sym is not None else None
        symdiff_m2 = 0.0
        worst_lobe = None
        if sym_er is not None and not sym_er.IsEmpty():
            parts = ([sym_er.GetGeometryRef(i)
                      for i in range(sym_er.GetGeometryCount())]
                     if sym_er.GetGeometryCount() else [sym_er])
            best_a = 0.0
            best_c = None
            for p in parts:
                a = p.GetArea()
                symdiff_m2 += a
                if a > best_a:
                    best_a = a
                    c = p.Centroid()
                    best_c = (c.GetX(), c.GetY()) if c is not None else None
            if best_c is not None:
                worst_lobe = (best_a,
                              geomatch.unproject_local(best_c[0], best_c[1],
                                                       lon0, lat0))
    except Exception:                                   # noqa: BLE001
        return None
    lobe_floor = max(LOBE_FLOOR_M2, LOBE_FLOOR_X * tol * tol)
    obvious_iou = iou < AREA_FLAG_IOU
    obvious_lobe = worst_lobe is not None and worst_lobe[0] >= lobe_floor
    ok_band = iou >= AREA_PASS_IOU and not obvious_lobe
    return {
        "flag": bool(not ok_band and (obvious_iou or obvious_lobe)),
        "iou": iou,
        "symdiff_m2": symdiff_m2,
        "worst_lobe": worst_lobe,
        "chart_m2": chart.GetArea(),
    }


def match_areas(chart_areas, source_areas, tol_m):
    """Pair charted polygons with source polygons by best IoU (bbox prefilter).
    Inputs are lists of ``[[exterior, *holes], …]`` per feature (each feature
    may be a multipolygon). Returns (pairs, unmatched_chart, unmatched_source).
    None when OGR/GEOS is unavailable."""
    o = _ogr()
    if o is None:
        return None
    if not chart_areas or not source_areas:
        return ([], list(range(len(chart_areas))),
                list(range(len(source_areas))))
    # chart_areas[0] = feature → [0][0] polygon → [0][0][0] ring →
    # [0][0][0][0] the first (lon, lat) point: unpack it as the local origin.
    lon0, lat0 = chart_areas[0][0][0][0]
    ch = [_multi_from_areas(o, _project_areas(a, lon0, lat0))
          for a in chart_areas]
    so = [_multi_from_areas(o, _project_areas(a, lon0, lat0))
          for a in source_areas]

    def env(g):
        x0, x1, y0, y1 = g.GetEnvelope()
        return (x0, y0, x1, y1)

    ch_e, so_e = [env(g) for g in ch], [env(g) for g in so]
    pairs = []
    used = set()
    for i, cg in enumerate(ch):
        best, best_j = 0.0, -1
        for j, sg in enumerate(so):
            if j in used:
                continue
            a, b = ch_e[i], so_e[j]
            if a[2] < b[0] or b[2] < a[0] or a[3] < b[1] or b[3] < a[1]:
                continue
            try:
                inter = cg.Intersection(sg)
                union = cg.Union(sg)
                ia = inter.GetArea() if inter is not None else 0.0
                ua = union.GetArea() if union is not None else 0.0
            except Exception:                           # noqa: BLE001
                continue
            iou = (ia / ua) if ua > 0 else 0.0
            if iou > best:
                best, best_j = iou, j
        if best_j >= 0 and best >= MATCH_MIN_IOU:
            pairs.append((i, best_j))
            used.add(best_j)
    matched = {i for i, _ in pairs}
    return (pairs,
            [i for i in range(len(ch)) if i not in matched],
            [j for j in range(len(so)) if j not in used])
