# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
s57/geometry.py
---------------
Native geometry **assembly** — the derived-view layer above the S-57 record
model. Points (PRIM=1), lines (PRIM=2), and areas (PRIM=3) are all assembled
natively; GDAL is used only as an *oracle* to verify the result.

This is the read-side mirror of the eventual writer. The node/edge records in
``s57.records`` are the *canonical, lossless* source of truth (raw integer
coordinates, COMF/SOMF, orientation, masking); the geometry values built here
are a **derived view** — georeferenced floats produced on demand, never stored
back. Nothing in this module mutates the record model.

What it emits is a small set of **native geometry values** (``Point``,
``MultiPoint``, ``LineString``, ``MultiLineString``, ``Polygon``,
``MultiPolygon``) — not OGR/GEOS objects. The engine stays GDAL-free; GDAL is used
only as an *oracle* to verify these values (see ``oracle_harness.py``) and,
behind ``s57.geomops.GeomOps``, as the optional calculator for set-theoretic
spatial predicates. The verification currency is WKT, so
each value can ``.wkt()`` itself, and this module also carries a small,
dependency-free WKT comparator (``parse_wkt`` / ``wkt_match``) used by the oracle
diff and by the headless tests.

Georeferencing (S-57 Part 3 / App B.1, via DSPM):
    lon = XCOO / COMF      lat = YCOO / COMF      depth = VE3D / SOMF
COMF is typically 1e7 (→ 7 decimal places), SOMF typically 10 (→ 1 dp).
Record coords are stored ``(y, x)`` or ``(y, x, z)`` (see ``records._vector``).

Build order: points → lines → areas — all three implemented. The line/area
seam is ``_edge_chain``, which a future writer inverts.
"""

import math
import re

try:                                    # package context
    from .records import _i, RCNM_VI, RCNM_VC, RCNM_VE
except (ImportError, ValueError):       # headless / top-level
    from records import _i, RCNM_VI, RCNM_VC, RCNM_VE   # noqa: F401

# S-57 pointer enumerations (Part 3). Centralised so read and the future writer
# share one vocabulary.
ORNT_FORWARD = 1
ORNT_REVERSE = 2
USAG_EXTERIOR = 1
USAG_INTERIOR = 2
USAG_EXTERIOR_TRUNCATED = 3
TOPI_BEGIN = 1
TOPI_END = 2

# S-57 typical multipliers, assumed only when DSPM omits them (builder flags it).
_DEFAULT_COMF = 1e7
_DEFAULT_SOMF = 10.0


# ---------------------------------------------------------------------------
# Formatting helpers
# ---------------------------------------------------------------------------
def _decimals_for(multiplier, fallback):
    """Decimal places implied by a power-of-ten multiplier (1e7 → 7).
    Non-power-of-ten or unusable input → ``fallback``."""
    try:
        m = float(multiplier)
    except (TypeError, ValueError):
        return fallback
    if m <= 0:
        return fallback
    d = math.log10(m)
    r = round(d)
    if abs(d - r) < 1e-9:
        return max(0, int(r))
    return fallback


def _fmt(value, decimals):
    """Deterministic fixed-point formatting; normalises a signed zero to '0...'."""
    s = f"{value:.{decimals}f}"
    if float(s) == 0.0:                 # kill '-0.0000000'
        s = f"{0.0:.{decimals}f}"
    return s


# ---------------------------------------------------------------------------
# Native geometry values (derived view)
# ---------------------------------------------------------------------------
class Geometry:
    """Base for native geometry values. Carries the formatting precision the
    builder derived from COMF/SOMF so ``.wkt()`` is deterministic and parameter
    free."""
    __slots__ = ()
    kind = "GEOMETRY"

    def wkt(self):
        raise NotImplementedError

    def coords(self):
        """Flat ordered list of coordinate tuples (x, y) or (x, y, z)."""
        raise NotImplementedError


class Point(Geometry):
    __slots__ = ("x", "y", "z", "xy_dp", "z_dp")
    kind = "POINT"

    def __init__(self, x, y, z=None, xy_dp=7, z_dp=1):
        self.x = x
        self.y = y
        self.z = z
        self.xy_dp = xy_dp
        self.z_dp = z_dp

    @property
    def has_z(self):
        return self.z is not None

    def coords(self):
        if self.z is None:
            return [(self.x, self.y)]
        return [(self.x, self.y, self.z)]

    def wkt(self):
        x, y = _fmt(self.x, self.xy_dp), _fmt(self.y, self.xy_dp)
        if self.z is None:
            return f"POINT ({x} {y})"
        return f"POINT Z ({x} {y} {_fmt(self.z, self.z_dp)})"

    def __repr__(self):
        return f"<Point {self.wkt()}>"


class MultiPoint(Geometry):
    __slots__ = ("points", "z_present", "xy_dp", "z_dp")
    kind = "MULTIPOINT"

    def __init__(self, points, z_present=False, xy_dp=7, z_dp=1):
        self.points = list(points)      # [(x, y)] or [(x, y, z)]
        self.z_present = z_present
        self.xy_dp = xy_dp
        self.z_dp = z_dp

    @property
    def has_z(self):
        return self.z_present

    def coords(self):
        return list(self.points)

    def wkt(self):
        xy, z = self.xy_dp, self.z_dp
        if self.z_present:
            inner = ",".join(
                f"({_fmt(p[0], xy)} {_fmt(p[1], xy)} {_fmt(p[2], z)})"
                for p in self.points)
            return f"MULTIPOINT Z ({inner})"
        inner = ",".join(f"({_fmt(p[0], xy)} {_fmt(p[1], xy)})"
                         for p in self.points)
        return f"MULTIPOINT ({inner})"

    def __repr__(self):
        return f"<MultiPoint n={len(self.points)} z={self.z_present}>"


def _pairs(seq, dp):
    return ",".join(f"{_fmt(p[0], dp)} {_fmt(p[1], dp)}" for p in seq)


def _ring(seq, dp):
    return "(" + _pairs(seq, dp) + ")"


class LineString(Geometry):
    __slots__ = ("points", "xy_dp")
    kind = "LINESTRING"

    def __init__(self, points, xy_dp=7):
        self.points = list(points)              # [(x, y)]
        self.xy_dp = xy_dp

    def coords(self):
        return list(self.points)

    def wkt(self):
        return f"LINESTRING ({_pairs(self.points, self.xy_dp)})"

    def __repr__(self):
        return f"<LineString n={len(self.points)}>"


class MultiLineString(Geometry):
    __slots__ = ("parts", "xy_dp")
    kind = "MULTILINESTRING"

    def __init__(self, parts, xy_dp=7):
        self.parts = [list(p) for p in parts]   # [[(x, y)], ...]
        self.xy_dp = xy_dp

    def coords(self):
        out = []
        for p in self.parts:
            out.extend(p)
        return out

    def wkt(self):
        inner = ",".join(_ring(p, self.xy_dp) for p in self.parts)
        return f"MULTILINESTRING ({inner})"

    def __repr__(self):
        return f"<MultiLineString parts={len(self.parts)}>"


class Polygon(Geometry):
    """Exterior ring first, holes after. Rings are closed (first == last)."""
    __slots__ = ("rings", "xy_dp")
    kind = "POLYGON"

    def __init__(self, rings, xy_dp=7):
        self.rings = [list(r) for r in rings]   # [exterior, *holes]
        self.xy_dp = xy_dp

    def coords(self):
        out = []
        for r in self.rings:
            out.extend(r)
        return out

    def wkt(self):
        inner = ",".join(_ring(r, self.xy_dp) for r in self.rings)
        return f"POLYGON ({inner})"

    def __repr__(self):
        return f"<Polygon rings={len(self.rings)}>"


class MultiPolygon(Geometry):
    __slots__ = ("polygons", "xy_dp")
    kind = "MULTIPOLYGON"

    def __init__(self, polygons, xy_dp=7):
        # [[exterior, *holes], ...]
        self.polygons = [[list(r) for r in poly] for poly in polygons]
        self.xy_dp = xy_dp

    def coords(self):
        out = []
        for poly in self.polygons:
            for r in poly:
                out.extend(r)
        return out

    def wkt(self):
        polys = ",".join(
            "(" + ",".join(_ring(r, self.xy_dp) for r in poly) + ")"
            for poly in self.polygons)
        return f"MULTIPOLYGON ({polys})"

    def __repr__(self):
        return f"<MultiPolygon n={len(self.polygons)}>"


def _raw_same(a, b):
    """Two raw (y, x[, z]) coordinate tuples share the same node position."""
    return (a is not None and b is not None
            and a[0] == b[0] and a[1] == b[1])


# ---------------------------------------------------------------------------
# Builder
# ---------------------------------------------------------------------------
class GeometryBuilder:
    """Assemble native geometry from one parsed ``S57Cell``.

    Implements ``point()`` (PRIM=1), ``line()`` (PRIM=2) and ``area()`` (PRIM=3),
    with ``geometry()`` dispatching by primitive. The builder reads COMF/SOMF
    from the cell's DSPM once and reuses them. It is read-only over the record
    model and robust to dirty data — a feature it cannot assemble yields ``None``
    rather than raising, so a single bad record never aborts a sweep.
    """

    def __init__(self, cell):
        self.cell = cell
        dspm = cell.dspm or {}
        self.comf = self._num(dspm.get("COMF"), _DEFAULT_COMF)
        self.somf = self._num(dspm.get("SOMF"), _DEFAULT_SOMF)
        self.coun = _i(dspm.get("COUN"))            # 1 = lat/long
        self.xy_dp = _decimals_for(self.comf, 7)
        self.z_dp = _decimals_for(self.somf, 1)
        self.notes = []                              # provisional / assumption flags
        if not dspm or dspm.get("COMF") in (None, ""):
            self.notes.append("DSPM COMF absent; assumed 1e7")
        if self.coun not in (None, 1):
            self.notes.append(
                f"DSPM COUN={self.coun} (not 1=lat/long); "
                "coordinates still divided by COMF")

    # -- helpers --------------------------------------------------------------
    @staticmethod
    def _num(value, default):
        if value is None or value == "":
            return default
        try:
            return float(value)
        except (TypeError, ValueError):
            return default

    def _resolve(self, rcnm, rcid):
        return self.cell.vectors.get((rcnm, rcid))

    def _georef(self, c):
        """Raw (y, x[, z]) ints → (lon, lat[, depth]) floats."""
        y, x = c[0], c[1]
        lon = x / self.comf
        lat = y / self.comf
        if len(c) >= 3 and c[2] is not None:
            return (lon, lat, c[2] / self.somf)
        return (lon, lat)

    # -- Phase 1: points ------------------------------------------------------
    def point(self, feature):
        """Assemble a point feature (PRIM=1) → ``Point`` / ``MultiPoint`` / None.

        Gathers coordinates from the feature's FSPT-referenced spatial records
        (an isolated node VI for ordinary points and soundings, or a connected
        node VC). Returns ``None`` for a non-point feature, a feature with no
        usable spatial reference, or dangling/empty pointers.

        Type rule: a point feature whose coordinates carry depth (Z) is a
        sounding set and is emitted as ``MULTIPOINT Z`` — including the
        degenerate single-sounding case — to match GDAL's S-57 driver, which
        always represents soundings as multipoints. A 2-D feature with one
        coordinate is a ``Point``; with several, a (rare) ``MULTIPOINT``.
        """
        if feature is None or feature.prim != 1:
            return None
        raw = []
        for ptr in feature.fspt:
            vec = self._resolve(ptr[0], ptr[1])
            if vec is None or not vec.coords:
                continue
            raw.extend(vec.coords)
        if not raw:
            return None

        has_z = any(len(c) >= 3 and c[2] is not None for c in raw)
        pts = [self._georef(c) for c in raw]
        if has_z:
            # Normalise every tuple to 3-D so MULTIPOINT Z is well-formed even if
            # a stray tuple lacked VE3D (defensive — SG3D normally always has it).
            pts = [(p[0], p[1], p[2] if len(p) >= 3 else 0.0) for p in pts]

        if has_z or len(pts) > 1:
            return MultiPoint(pts, z_present=has_z,
                              xy_dp=self.xy_dp, z_dp=self.z_dp)
        x, y = pts[0]
        return Point(x, y, None, xy_dp=self.xy_dp, z_dp=self.z_dp)

    def iter_points(self):
        """Yield ``(feature, geometry)`` for every PRIM=1 feature that assembles.
        Convenience for the oracle harness and batch verification."""
        for f in self.cell.features:
            if f.prim == 1:
                g = self.point(f)
                if g is not None:
                    yield f, g

    def iter_isolated_node_points(self):
        """Yield ``(feature, geometry)`` for PRIM=1 features whose spatial
        reference includes an isolated node (VI, RCNM 110) — the v1 GIS-bridge
        scope. Soundings sit on isolated nodes too, so they appear here as
        ``MultiPoint`` (Z); a consumer wanting only single points can filter on
        ``isinstance(geom, Point)``."""
        for f in self.cell.features:
            if f.prim != 1:
                continue
            if not any(p[0] == RCNM_VI for p in f.fspt):
                continue
            g = self.point(f)
            if g is not None:
                yield f, g

    # -- the read/write seam --------------------------------------------------
    def _edge_chain(self, rcnm, rcid, ornt):
        """Resolve one edge (VE) pointer to its full ordered raw-coordinate run
        ``[begin, *interior, end]`` as (y, x[, z]) ints, oriented per ``ornt``.

        This is the single seam line and area assembly share, and the exact
        point a future **writer** inverts: given a ring/line it decomposes back
        into a VE (interior SG2D) plus begin/end connected-node (VC) pointers and
        an ORNT. Begin/end come from the edge's VRPT (TOPI begin/end); the
        interior comes from the edge's own SG2D. Returns ``None`` if the edge or
        both end nodes are unresolvable.
        """
        ve = self._resolve(rcnm, rcid)
        if ve is None:
            return None
        begin = end = None
        for vp in ve.vrpt:
            topi = vp[4] if len(vp) > 4 else None
            node = self._resolve(vp[0], vp[1])
            if node is None or not node.coords:
                continue
            if topi == TOPI_BEGIN:
                begin = node.coords[0]
            elif topi == TOPI_END:
                end = node.coords[0]
        seq = []
        if begin is not None:
            seq.append(begin)
        seq.extend(ve.coords)
        if end is not None:
            seq.append(end)
        if not seq:
            return None
        if ornt == ORNT_REVERSE:
            seq = list(reversed(seq))
        return seq

    # -- Phase 2: lines -------------------------------------------------------
    def line(self, feature):
        """Assemble a line feature (PRIM=2) → ``LineString`` /
        ``MultiLineString`` / None. Chains the feature's edges head-to-tail
        (dropping the shared node at each join); a break in connectivity starts a
        new part, so a feature with disjoint edges becomes a ``MultiLineString``,
        matching GDAL."""
        if feature is None or feature.prim != 2:
            return None
        parts, cur = [], []
        for ptr in feature.fspt:
            if ptr[0] != RCNM_VE:
                continue
            seq = self._edge_chain(ptr[0], ptr[1], ptr[2])
            if not seq:
                continue
            if not cur:
                cur = list(seq)
            elif _raw_same(cur[-1], seq[0]):
                cur.extend(seq[1:])
            elif _raw_same(cur[-1], seq[-1]):
                cur.extend(list(reversed(seq))[1:])
            else:
                parts.append(cur)
                cur = list(seq)
        if cur:
            parts.append(cur)
        geo = [[self._georef(c)[:2] for c in p] for p in parts]
        geo = [p for p in geo if len(p) >= 2]
        if not geo:
            return None
        if len(geo) == 1:
            return LineString(geo[0], xy_dp=self.xy_dp)
        return MultiLineString(geo, xy_dp=self.xy_dp)

    # -- Phase 3: areas -------------------------------------------------------
    def area(self, feature):
        """Assemble an area feature (PRIM=3) → ``Polygon`` / ``MultiPolygon`` /
        None. Edges are chained into closed rings (respecting connectivity),
        split by USAG into exterior shells and interior holes; holes are assigned
        to their containing shell by point-in-polygon, so multiple shells become
        a ``MultiPolygon``."""
        if feature is None or feature.prim != 3:
            return None
        segs = []
        for ptr in feature.fspt:
            if ptr[0] != RCNM_VE:
                continue
            seq = self._edge_chain(ptr[0], ptr[1], ptr[2])
            if seq:
                usag = ptr[3] if len(ptr) > 3 else USAG_EXTERIOR
                segs.append((usag, seq))
        if not segs:
            return None
        rings = self._build_rings(segs)
        geo_rings = []
        for usag, ring in rings:
            g = [self._georef(c)[:2] for c in ring]
            if len(g) < 3:
                continue
            if g[0] != g[-1]:
                g.append(g[0])                  # close the ring
            geo_rings.append((usag, g))
        if not geo_rings:
            return None
        return self._rings_to_polygons(geo_rings)

    @staticmethod
    def _build_rings(segs):
        """Chain ``(usag, raw_seq)`` edge runs into closed rings using exact node
        identity, following FSPT order and tolerating a flipped edge. A ring is
        emitted when it returns to its start; a connectivity break starts a new
        ring. Each ring keeps the USAG of the edge that opened it."""
        rings, cur, cur_usag = [], [], None
        for usag, seq in segs:
            if not seq:
                continue
            if not cur:
                cur, cur_usag = list(seq), usag
            elif _raw_same(cur[-1], seq[0]):
                cur.extend(seq[1:])
            elif _raw_same(cur[-1], seq[-1]):
                cur.extend(list(reversed(seq))[1:])
            else:
                rings.append((cur_usag, cur))
                cur, cur_usag = list(seq), usag
            if len(cur) >= 4 and _raw_same(cur[0], cur[-1]):
                rings.append((cur_usag, cur))
                cur, cur_usag = [], None
        if len(cur) >= 3:
            rings.append((cur_usag, cur))
        return rings

    def _rings_to_polygons(self, geo_rings):
        """Group closed (usag, ring) pairs into a Polygon or MultiPolygon:
        exteriors are shells, USAG=interior rings are holes assigned to the
        containing shell."""
        exteriors = [r for u, r in geo_rings if u != USAG_INTERIOR]
        interiors = [r for u, r in geo_rings if u == USAG_INTERIOR]
        if not exteriors:                       # degenerate: no marked exterior
            exteriors, interiors = [r for _, r in geo_rings], []
        if len(exteriors) == 1:
            return Polygon([exteriors[0]] + interiors, xy_dp=self.xy_dp)
        polys = [[ext] for ext in exteriors]
        for hole in interiors:
            idx = self._containing(exteriors, hole)
            polys[idx if idx is not None else 0].append(hole)
        return MultiPolygon(polys, xy_dp=self.xy_dp)

    @staticmethod
    def _point_in_ring(pt, ring):
        x, y = pt[0], pt[1]
        inside = False
        n = len(ring)
        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

    def _containing(self, exteriors, hole):
        if not hole:
            return None
        for i, r in enumerate(exteriors):
            if self._point_in_ring(hole[0], r):
                return i
        return None

    # -- unified dispatch -----------------------------------------------------
    def geometry(self, feature):
        """Assemble any feature by its PRIM: 1→point, 2→line, 3→area. None for
        a feature with no primitive or one that can't be assembled."""
        if feature is None:
            return None
        if feature.prim == 1:
            return self.point(feature)
        if feature.prim == 2:
            return self.line(feature)
        if feature.prim == 3:
            return self.area(feature)
        return None

    def iter_geometries(self):
        """Yield ``(feature, geometry)`` for every feature that assembles, across
        all primitives — the full-geometry counterpart to ``iter_points``."""
        for f in self.cell.features:
            g = self.geometry(f)
            if g is not None:
                yield f, g


# ---------------------------------------------------------------------------
# Dependency-free WKT comparator (verification currency)
# ---------------------------------------------------------------------------
_NUM = re.compile(r"[-+]?\d*\.?\d+(?:[eE][-+]?\d+)?")
_TYPE = re.compile(r"\s*([A-Za-z]+)\s*(ZM|Z|M)?")


def parse_wkt(wkt):
    """Lightweight parse → ``(base_type, has_z, coords)``.

    ``base_type`` is upper-case (POINT/MULTIPOINT/LINESTRING/...); ``coords`` is
    a flat, ordered list of float tuples. Dimensionality is taken from an
    explicit ``Z`` tag (ISO WKT — what ``geom.ExportToIsoWkt()`` and this
    module's ``.wkt()`` both emit). When no tag is present, dimension is inferred
    from the coordinate count only if unambiguous (odd count divisible by 3).
    Ring/part structure is flattened — sufficient for the point and line phases;
    the area phase will add structure-aware ring comparison.
    """
    if not wkt:
        return "", False, []
    s = wkt.strip()
    m = _TYPE.match(s)
    base = m.group(1).upper() if m else ""
    has_z = bool(m and m.group(2) and "Z" in m.group(2).upper())
    tail = s[m.end():] if m else s
    nums = [float(x) for x in _NUM.findall(tail)]
    n = len(nums)
    stride = 3 if has_z else 2
    if not has_z and n and n % 2 != 0 and n % 3 == 0:
        stride, has_z = 3, True          # unambiguous untagged 3-D
    coords = [tuple(nums[i:i + stride])
              for i in range(0, n - stride + 1, stride)]
    return base, has_z, coords


def _paren_groups(s):
    """Inner text of each top-level parenthesised group in ``s``."""
    groups, depth, start = [], 0, None
    for i, ch in enumerate(s):
        if ch == "(":
            if depth == 0:
                start = i + 1
            depth += 1
        elif ch == ")":
            depth -= 1
            if depth == 0 and start is not None:
                groups.append(s[start:i])
    return groups


def _ring_coords(text):
    out = []
    for pair in text.split(","):
        nums = _NUM.findall(pair)
        if len(nums) >= 2:
            out.append((float(nums[0]), float(nums[1])))
    return out


def parse_rings(wkt):
    """Parse a POLYGON/MULTIPOLYGON into ``[polygon]`` where each polygon is
    ``[ring]`` and each ring is ``[(x, y)]`` — preserving ring structure (which
    the flat ``parse_wkt`` discards), so areas can be compared structurally."""
    if not wkt:
        return []
    s = wkt.strip()
    m = _TYPE.match(s)
    base = m.group(1).upper() if m else ""
    body = s[m.end():] if m else s
    if base == "POLYGON":
        inner = _paren_groups(body)
        if not inner:
            return []
        return [[_ring_coords(r) for r in _paren_groups(inner[0])]]
    if base == "MULTIPOLYGON":
        outer = _paren_groups(body)
        if not outer:
            return []
        polys = _paren_groups(outer[0])
        return [[_ring_coords(r) for r in _paren_groups(p)] for p in polys]
    return []


def _canonical_ring(ring, dp=7):
    """A ring's identity independent of start vertex and winding: drop the
    closing duplicate, round to ``dp``, then take the lexicographically smallest
    rotation across both directions. Two geometrically identical rings share
    one canonical form."""
    pts = [(round(p[0], dp), round(p[1], dp)) for p in ring]
    if len(pts) > 1 and pts[0] == pts[-1]:
        pts = pts[:-1]
    if not pts:
        return ()

    def rot_min(seq):
        k = min(range(len(seq)), key=lambda i: seq[i])
        return tuple(seq[k:] + seq[:k])

    return min(rot_min(pts), rot_min(list(reversed(pts))))


def rings_match(a_wkt, b_wkt, dp=7):
    """True if two POLYGON/MULTIPOLYGON WKTs describe the same shape — same set
    of polygons, each the same set of rings — independent of ring start vertex,
    winding, and ring order."""
    pa, pb = parse_rings(a_wkt), parse_rings(b_wkt)
    if len(pa) != len(pb):
        return False

    def norm(polys):
        return sorted(tuple(sorted(_canonical_ring(r, dp) for r in poly))
                      for poly in polys)

    return norm(pa) == norm(pb)


def _coords_close(a, b, tol, ordered):
    if len(a) != len(b):
        return False
    if not ordered:
        a, b = sorted(a), sorted(b)
    for ca, cb in zip(a, b):
        if len(ca) != len(cb):
            return False
        if any(abs(u - v) > tol for u, v in zip(ca, cb)):
            return False
    return True


def wkt_match(a_wkt, b_wkt, tol=1e-7, ordered=True):
    """Compare two WKT strings structurally + numerically within ``tol``.

    Returns a dict: ``coord_match`` (the load-bearing result — same coordinates
    within tolerance), ``type_match`` (same base type and Z-ness), ``match`` (an
    alias of ``coord_match``), and a human ``detail`` string. Coordinate equality
    is what proves assembly correctness; a type-only difference (e.g. our
    ``POINT Z`` vs an oracle ``MULTIPOINT Z`` for a single sounding) is reported
    separately rather than failing the coordinate test.
    """
    ta, za, ca = parse_wkt(a_wkt)
    tb, zb, cb = parse_wkt(b_wkt)
    type_match = (ta == tb) and (za == zb)
    # Areas: compare ring structure (rotation/winding/order invariant).
    if ta in ("POLYGON", "MULTIPOLYGON") or tb in ("POLYGON", "MULTIPOLYGON"):
        coord_match = rings_match(a_wkt, b_wkt, dp=7)
    else:
        coord_match = _coords_close(ca, cb, tol, ordered)
        # A line equals its reverse — accept either direction.
        if (not coord_match and ordered
                and ta in ("LINESTRING", "MULTILINESTRING")):
            coord_match = _coords_close(ca, list(reversed(cb)), tol, ordered)
    if coord_match and type_match:
        detail = "ok"
    elif coord_match:
        detail = f"coords match; type differs ({ta}{'Z' if za else ''} vs {tb}{'Z' if zb else ''})"
    elif len(ca) != len(cb):
        detail = f"coordinate count differs ({len(ca)} vs {len(cb)})"
    else:
        detail = "coordinate values differ beyond tolerance"
    return {"match": coord_match, "coord_match": coord_match,
            "type_match": type_match, "detail": detail}
