# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
enc_loader.py
Loads selected S-57 object layers from multiple ENC files and merges
them into a single GeoPackage on disk (one table per layer name).
Layers are placed FLAT into a single "ENC Layers" parent group, ordered by
the priority-number prefix on each layer name (no per-group sub-groups).
"""

import os
import re
import tempfile
import uuid
from collections import OrderedDict
from qgis.core import (
    QgsVectorLayer,
    QgsVectorFileWriter,
    QgsProject,
    QgsMessageLog,
    Qgis,
    QgsProcessingFeedback,
    QgsWkbTypes,
    QgsField,
    QgsFields,
    QgsFeature,
    QgsSpatialIndex,
    QgsGeometry,
    QgsPointXY,
    QgsPoint,
    QgsMultiPoint,
    QgsRectangle,
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsExpressionContextUtils,
)
from qgis.PyQt.QtCore import QVariant
from .enc_config import get_config
from .enc_scanner import read_dsid_full, usage_band


# Geometry-kind suffixes and their draw-priority offset (points drawn on top).
_GEOM_KINDS = {
    "point": 2,
    "line":  1,
    "area":  0,
}


def _geom_kind(layer):
    """Return 'point' | 'line' | 'area' | None for a QgsVectorLayer."""
    gt = QgsWkbTypes.geometryType(layer.wkbType())
    if gt == QgsWkbTypes.PointGeometry:
        return "point"
    if gt == QgsWkbTypes.LineGeometry:
        return "line"
    if gt == QgsWkbTypes.PolygonGeometry:
        return "area"
    return None


def _feature_geom_kind(feat):
    """Return 'point' | 'line' | 'area' | None for a single feature's actual
    geometry. Used to catch the rare case where an S-57 layer is declared as
    one geometry type but a stray feature is another (e.g. a point inside an
    LNDARE area layer) — such a feature must not be written into a table typed
    for the layer's declared kind, or the GeoPackage write drops it noisily."""
    geom = feat.geometry()
    if geom is None or geom.isEmpty():
        return None
    gt = QgsWkbTypes.geometryType(geom.wkbType())
    if gt == QgsWkbTypes.PointGeometry:
        return "point"
    if gt == QgsWkbTypes.LineGeometry:
        return "line"
    if gt == QgsWkbTypes.PolygonGeometry:
        return "area"
    return None


def _held_back_native_reason(fpath, lname, cat_dir, rcids=None):
    """Native, GDAL-free explanation of *why* features of class ``lname`` in the
    cell ``fpath`` are problematic — runs the geomcheck diagnostic and returns a
    short breakdown for that class (e.g. ``"1 degenerate-area"`` or
    ``"2 unassemblable"``), or ``""`` if it finds nothing or cannot run. This is
    the native engine's independent read, so it confirms that a feature GDAL
    held back from the area table is a genuine geometry defect rather than a
    loader artefact. Never raises.

    ``rcids`` scopes the check to the features actually held back: without it the
    geomcheck reports every issue in the class (so a cell with one held-back
    feature but many other merely-OGC-invalid features of the same class would
    annotate the one-feature held-back line with the whole class's defects, which
    overstates it). Passing the held-back RCIDs makes the breakdown match the
    held-back count."""
    try:
        try:
            from .validation import Catalogue
            from .s57.geomcheck import check_paths
        except (ImportError, ValueError):
            from validation import Catalogue
            from s57.geomcheck import check_paths
        out = check_paths([fpath], catalogue=Catalogue(cat_dir), rcids=rcids)
        if not out:
            return ""
        _summary, issues = out[0]
        from collections import Counter
        kinds = Counter(i.kind for i in issues if i.klass == lname)
        return ", ".join(f"{n} {k}" for k, n in sorted(kinds.items()))
    except Exception:
        return ""


class _gdal_quiet:
    """Context manager that silences GDAL's low-level error handler for the
    duration of a block, then restores it. Used to suppress harmless,
    expected per-feature CRITICALs (e.g. dropping a geometrically-incompatible
    feature from a single-type table) without hiding real failures, which are
    still reported via QgsVectorFileWriter result codes. Falls back to a no-op
    if GDAL's Python error API isn't available.
    """

    def __enter__(self):
        self._ok = False
        try:
            from osgeo import gdal
            self._gdal = gdal
            gdal.PushErrorHandler("CPLQuietErrorHandler")
            self._ok = True
        except Exception:
            self._gdal = None
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        if self._ok and self._gdal is not None:
            try:
                self._gdal.PopErrorHandler()
            except Exception:
                pass
        return False


def _log(msg, level=Qgis.Info):
    QgsMessageLog.logMessage(msg, "ENCore GIS", level)


class _StageFeedback(QgsProcessingFeedback):
    """Forwards a single processing algorithm's 0-100% into an overall
    sub-progress range [base, base+span], via subprogress_fn(percent, label)."""

    def __init__(self, subprogress_fn, base, span, label):
        super().__init__()
        self._fn = subprogress_fn
        self._base = base
        self._span = span
        self._label = label

    def setProgress(self, pct):
        if self._fn:
            overall = self._base + (pct / 100.0) * self._span
            self._fn(int(overall), self._label)


def load_layers(layer_index, selection, feedback_fn=None, progress_fn=None,
                subprogress_fn=None, target_crs="EPSG:3857",
                want_coverage=False, coverage_usages=None, best_files=None,
                beautiful=False, cancel_fn=None):
    """Load the user's selected object layers into QGIS, styled and reprojected.

    Builds one session GeoPackage from the chosen per-usage layers (merged, and
    in Beautiful mode run through each layer's processing recipe), applies the
    bundled QML styles, and reprojects everything to target_crs. Honours the
    live Chart Control variables and the Fast/Beautiful render mode. Progress and
    status flow back through the optional feedback/progress callbacks (used by
    LoadWorker). Returns the session GeoPackage path and the loaded-layer
    manifest — the values LoadWorker relays via its `finished` signal."""
    # `selection` is a list of (usage_band, code) pairs — the exact per-usage
    # layers the user chose. `coverage_usages` is the set of usage bands the
    # user chose ENC Coverage for. Both are per-usage (option (b)), so we load
    # ONLY what was selected, never every usage a code happens to appear in.
    # `beautiful` is the pre-load render mode: True runs each layer's processing
    # recipe (fix → dissolve → singleparts; seam-free, slower); False (Fast,
    # the default) loads the raw merged geometry + symbology, cell seams and all.
    coverage_usages = set(coverage_usages or set())

    def log(msg, level=Qgis.Info):
        _log(msg, level)
        if feedback_fn:
            feedback_fn(msg)

    def progress(cur, total):
        if progress_fn:
            progress_fn(cur, total)

    cfg = get_config()

    # NOTE: this loader no longer uses GDAL's S-57 driver. Object classes are
    # decoded by the in-house ISO 8211 engine (iso8211/ + s57/) against the
    # bundled catalogue (validation.Catalogue), which natively covers the
    # Inland-ENC classes (e.g. OBJL 17029 `boylat`) WITH their attributes — so
    # the old OGR_S57_OPTIONS / S57_CSV config injection is gone. GDAL is still
    # used only for GeoPackage I/O via QGIS's "ogr" provider (output side).

    # All loaded layers are reprojected to this CRS during the merge write, so
    # everything downstream (recipes, zoomLevel spacing, portrayal, the project
    # CRS QGIS adopts from the first layer) is in one projected, metric CRS.
    dest_crs = QgsCoordinateReferenceSystem(target_crs)
    if not dest_crs.isValid():
        log(f"Target CRS '{target_crs}' invalid — falling back to EPSG:3857",
            Qgis.Warning)
        dest_crs = QgsCoordinateReferenceSystem("EPSG:3857")
    log(f"Target CRS for all layers: {dest_crs.authid()}")

    # Each load writes to its OWN uniquely-named GeoPackage in the temp dir.
    # Earlier versions reused one fixed filename, but on Windows the previous
    # load's layers keep open handles to that file, so a second load could not
    # delete it ([WinError 32] ... used by another process) and appended to the
    # stale file instead. A unique name per load sidesteps the lock entirely:
    # the old file stays owned by its layers, the new load is always clean.
    tmp_dir = tempfile.gettempdir()
    gpkg_path = os.path.join(
        tmp_dir, f"enc_explorer_layers_{uuid.uuid4().hex}.gpkg"
    )

    # Best-effort tidy: remove leftover ENCore GIS GeoPackages from past
    # loads that are no longer locked (e.g. from a previous QGIS session). Any
    # still held open by live layers will fail to delete — that's fine, we just
    # skip them silently rather than logging a scary warning every load.
    removed = 0
    for fname in os.listdir(tmp_dir):
        if fname.startswith("enc_explorer_layers") and (
            fname.endswith(".gpkg") or ".gpkg-" in fname or
            fname.endswith(".gpkg-wal") or fname.endswith(".gpkg-shm") or
            fname.endswith(".gpkg-journal")
        ):
            stale = os.path.join(tmp_dir, fname)
            if stale == gpkg_path:
                continue
            try:
                os.remove(stale)
                removed += 1
            except Exception:
                pass   # still locked by a live layer; leave it be
    if removed:
        log(f"Cleaned up {removed} leftover temp file(s) from earlier loads.")

    log(f"Writing to: {gpkg_path}")

    total_layers = len(selection)

    # Layers are BUILT here (worker thread) but REGISTERED on the main thread by
    # register_layers(), to avoid the QGIS rule that map layers must be added to
    # the project from the main thread. We therefore collect a lightweight
    # MANIFEST describing each finished GeoPackage table — sort_key, table name,
    # display name, and how to style it — rather than the QgsVectorLayer objects
    # themselves. The main thread re-opens each table from the GeoPackage,
    # applies the style, adds it, and stamps the project CRS when the project is
    # empty. (Adding from the worker thread bypassed QGIS's "empty project
    # adopts the first layer's CRS" step, which is what produced the
    # intermittent displaced-vs-basemap bug.)
    manifest = []   # list of {sort_key, table, name, style_code, geom}

    def _manifest_table(layer):
        """The GeoPackage table name behind a file-backed layer, from its URI."""
        src = layer.source()
        if "layername=" in src:
            return src.split("layername=", 1)[1].split("|")[0]
        return None

    # ── v55: resolve buoy↔light relationships ONCE, up front. The result feeds
    #    both the AtoN merge (folding a light onto its buoy row) and the LIGHTS
    #    build (tagging buoyLight), so the two agree no matter which order they
    #    appear in the selection.
    _links = _compute_buoy_light_links(layer_index, selection, log)
    _claimed_lnams = _links["claimed_lnams"]
    _fold_by_usage = _links["fold_by_usage"]

    for sel_idx, (ub, lname) in enumerate(selection):
        # Cooperative cancel: bail out at a layer boundary if the caller asked to
        # stop (e.g. QGIS is closing). Building QgsVectorLayers / writing the GPKG
        # off the main thread must NOT keep running while QGIS tears down — that
        # races the QgsProject/QgsStyle singletons and crashes (heap corruption).
        if cancel_fn is not None and cancel_fn():
            log("[load] cancel requested — stopping at layer boundary.",
                Qgis.Warning)
            break
        # ── Composite Aids-to-Navigation layer ──────────────────────────────
        # "AtoN" is a synthetic selection that merges the floating buoy classes
        # and their topmarks into ONE point layer (see _build_aton_layer). It
        # gathers files containing ANY AtoN-family class at this usage, rather
        # than a single code's files.
        if lname == "AtoN":
            # Match family classes against the scan index case-insensitively:
            # Inland cells index their buoys under lowercase keys (`boylat`)
            # while the family list is the standard uppercase (`BOYLAT`). A cell
            # uses one spelling or the other, so this just widens the net.
            fam_upper = _ATON_BODY_CLASSES_UPPER | {_ATON_TOPMARK_CLASS.upper()}
            aton_files = []
            seen_paths = set()
            for key, files in layer_index.items():
                if key.upper() not in fam_upper:
                    continue
                for f in files:
                    if f in seen_paths:
                        continue
                    if usage_band(os.path.splitext(os.path.basename(f))[0]) == ub:
                        aton_files.append(f)
                        seen_paths.add(f)
            log(f"[{sel_idx + 1}/{total_layers}] U{ub} Aids to Navigation "
                f"(AtoN) — {len(aton_files)} file(s)")
            progress(sel_idx, total_layers)
            if not aton_files:
                log(f"  ✗ No AtoN source files at usage {ub}", Qgis.Warning)
                progress(sel_idx + 1, total_layers)
                continue
            try:
                aton_layer = _build_aton_layer(
                    aton_files, ub, dest_crs, gpkg_path, log,
                    base_priority=cfg.priority_for("BOYLAT") or 80,
                    subprogress_fn=subprogress_fn,
                    light_fold=_fold_by_usage.get(ub),
                )
                if aton_layer is not None:
                    sort_key = ub * 1000 + (cfg.priority_for("BOYLAT") or 80)
                    manifest.append({
                        "sort_key": sort_key,
                        "table": _manifest_table(aton_layer),
                        "name": aton_layer.name(),
                        "style_code": "AtoN",
                        "geom": "point",
                    })
            except Exception as e:
                log(f"  ✗ Could not build AtoN layer (usage {ub}): {e}",
                    Qgis.Warning)
            progress(sel_idx + 1, total_layers)
            continue

        source_files = layer_index.get(lname, [])
        display_name = cfg.display_name(lname)

        # Keep ONLY this usage band's files for this code — the user picked
        # (usage, layer) explicitly, so we never spill into other usages.
        usage_files = [
            f for f in source_files
            if usage_band(os.path.splitext(os.path.basename(f))[0]) == ub
        ]

        log(f"[{sel_idx + 1}/{total_layers}] U{ub} {display_name} ({lname}) — "
            f"{len(usage_files)} file(s)")
        progress(sel_idx, total_layers)

        if not usage_files:
            log(f"  ✗ No source files for {lname} at usage {ub}", Qgis.Warning)
            progress(sel_idx + 1, total_layers)
            continue

        base_priority = cfg.priority_for(lname)

        # (single usage per selected pair — the block below runs once)
        if True:

            # ── Merge this usage's source files into the GeoPackage, SPLIT BY
            #    GEOMETRY. A single S-57 class can hold points, lines AND
            #    polygons; a GeoPackage table holds one geometry type, so we
            #    write each kind to its own table, usage-qualified:
            #    u<usage>_<LNAME>_<point/line/area>.
            kinds = {}   # kind -> {"table","count","first"}

            # Complete, catalogue-driven attribute schema for this class (the
            # full S-57/Inland datamodel — Attribute_A∪B∪C — plus any extras
            # present in the data). Every cell is built against this shared
            # schema so each loaded layer carries its full datamodel (e.g. LIGHTS
            # always has SECTR1/SECTR2) and the GeoPackage merge — first cell
            # creates the table, the rest append with AppendToLayerNoNewFields —
            # never drops an attribute the first cell happened to lack.
            attr_union = _class_attr_union(usage_files, lname)

            for fpath in usage_files:
                fname = os.path.basename(fpath)
                src_layer = _open_ogr_layer(fpath, lname, attr_union)
                if src_layer is None or not src_layer.isValid():
                    continue
                feat_count = src_layer.featureCount()
                if feat_count == 0:
                    continue
                kind = _geom_kind(src_layer)
                if kind is None:
                    log(f"    ? {fname}: unknown geometry for {lname}, skipping",
                        Qgis.Warning)
                    continue

                state = kinds.setdefault(
                    kind, {"table": f"u{ub}_{lname}_{kind}",
                           "count": 0, "first": True}
                )

                # Guard against the rare S-57 quirk where a layer is declared as
                # one geometry kind but holds a stray feature of another (e.g. a
                # point inside an LNDARE area layer). If we wrote that feature
                # into the kind-typed table it would survive the merge but then
                # be rejected — noisily — by a strictly-typed processing step
                # like fixgeometries. So we select only the features whose
                # ACTUAL geometry matches this table's kind and write just those.
                #
                # Only area/line tables are at risk (a lower-dimension stray
                # landing in a higher-dimension table is what breaks); point
                # tables can't hit this, so we skip the per-feature scan for them
                # entirely and keep the common fast path untouched.
                only_selected = False
                write_count = feat_count
                if kind in ("area", "line"):
                    matching_ids, held_kinds, held_rcids = [], [], []
                    rcid_idx = src_layer.fields().lookupField("RCID")
                    for f in src_layer.getFeatures():
                        fk = _feature_geom_kind(f)
                        if fk == kind:
                            matching_ids.append(f.id())
                        else:
                            held_kinds.append(fk or "null")
                            if rcid_idx >= 0:
                                rv = f.attribute(rcid_idx)
                                if rv is not None:
                                    try:
                                        held_rcids.append(int(rv))
                                    except (ValueError, TypeError):
                                        pass
                    if held_kinds:
                        only_selected = True
                        dropped = len(held_kinds)
                        from collections import Counter
                        actual = ", ".join(
                            f"{n} {k}" for k, n in
                            sorted(Counter(held_kinds).items()))
                        _cat_dir = os.path.join(os.path.dirname(__file__),
                                                "s57catalogue")
                        # Scope the native reason to just the held-back features
                        # (their RCIDs) so the breakdown matches the held-back
                        # count rather than the whole class's geometry issues.
                        why = _held_back_native_reason(
                            fpath, lname, _cat_dir,
                            rcids=held_rcids or None)
                        why = f"; native geomcheck: {why}" if why else ""
                        src_layer.selectByIds(matching_ids)
                        log(f"    • {fname}: {dropped} feature(s) held back from "
                            f"{kind} table — actual geometry {actual}, not "
                            f"{kind}{why}", Qgis.Info)
                        write_count = len(matching_ids)

                transform_context = QgsProject.instance().transformContext()
                options = QgsVectorFileWriter.SaveVectorOptions()
                options.driverName = "GPKG"
                options.layerName = state["table"]
                options.fileEncoding = "UTF-8"
                options.onlySelectedFeatures = only_selected
                # Reproject this source layer to the target CRS in the same pass
                # as the merge write. S-57 arrives in EPSG:4326; everything
                # downstream then works in the target's projected metres.
                src_crs = src_layer.crs()
                if state["first"]:
                    # ── CRS DIAGNOSTIC ─ what GDAL actually assigned to this
                    #    S-57 source. If this is not EPSG:4326 (or is empty),
                    #    that's the lead for any global-misplacement issue.
                    _se = src_layer.extent()
                    log(f"    [crs] source {lname}/{kind}: "
                        f"{src_crs.authid() or 'EMPTY'} "
                        f"(valid={src_crs.isValid()}, "
                        f"geographic={src_crs.isGeographic()}) "
                        f"extent=[{_se.xMinimum():.2f},{_se.yMinimum():.2f} .. "
                        f"{_se.xMaximum():.2f},{_se.yMaximum():.2f}] "
                        f"→ target {dest_crs.authid()}", Qgis.Info)
                if src_crs.isValid() and src_crs != dest_crs:
                    options.ct = QgsCoordinateTransform(
                        src_crs, dest_crs, transform_context
                    )
                if state["first"] and not os.path.exists(gpkg_path):
                    options.actionOnExistingFile = \
                        QgsVectorFileWriter.CreateOrOverwriteFile
                elif state["first"]:
                    options.actionOnExistingFile = \
                        QgsVectorFileWriter.CreateOrOverwriteLayer
                else:
                    options.actionOnExistingFile = \
                        QgsVectorFileWriter.AppendToLayerNoNewFields

                # Real write failures are still surfaced via the result code
                # below; the quiet scope only suppresses harmless GDAL chatter.
                with _gdal_quiet():
                    result = QgsVectorFileWriter.writeAsVectorFormatV3(
                        src_layer, gpkg_path, transform_context, options
                    )
                if only_selected:
                    src_layer.removeSelection()
                if result[0] != QgsVectorFileWriter.NoError:
                    log(f"    ✗ Write failed for {state['table']}: {result[1]}",
                        Qgis.Warning)
                    continue
                state["count"] += write_count
                state["first"] = False

            written = {k: v for k, v in kinds.items() if v["count"] > 0}
            if not written:
                continue

            # ── v55: stamp the buoyLight flag on the LIGHTS layer (before it is
            #    loaded/styled), so portrayal can hide buoy-lights already drawn
            #    as lit buoys. Always added (False for all) when there are no
            #    matches, so a QML filter on buoyLight never hits a missing field.
            if lname.upper() == _ATON_LIGHT_CLASS:
                for _k, _state in written.items():
                    _tag_buoylights(gpkg_path, _state["table"],
                                    _claimed_lnams, log)

            # Process + portray + add each geometry piece for this usage.
            for kind, state in written.items():
                actual_table = _find_gpkg_table(gpkg_path, state["table"])
                if actual_table is None:
                    log(f"  ✗ Table {state['table']} not found in GPKG",
                        Qgis.Warning)
                    continue

                # Loaded layer name = "<usage> <priority> <human name> (<geom>)".
                #  - usage: the navigational-purpose band integer (1–8). Leads
                #    the name so the flat layer list self-groups by usage, then
                #    by priority within each usage — matching the descending
                #    stacking (higher usage + higher priority on top).
                #  - priority: base XML priority (0-90), zero-padded to 2 digits.
                #  - geom tag: ALWAYS shown for cross-load consistency.
                # Style lookup is by (code, geom) only — usage does NOT affect
                # symbology, so every usage of a class loads the same QML.
                prio = base_priority + _GEOM_KINDS.get(kind, 0)
                # Sort key for stacking: (usage desc, priority desc). Encode as a
                # single sortable integer so the existing flat-order sort works.
                sort_key = ub * 1000 + prio
                lyr_name = f"U{ub} | {base_priority:02d} {display_name} ({kind})"

                # 1) Processing recipe. Beautiful runs the whole recipe. Fast
                #    keeps only portrayal-ESSENTIAL steps and drops the cosmetic
                #    geometry cleanup (fix/dissolve/singleparts): so SOUNDG's
                #    vertex split + zoom-level field still run (soundings need
                #    them to portray), while area dissolves are left raw and
                #    cell seams stay visible.
                recipe = cfg.recipe_for(lname, kind)
                if recipe and not beautiful:
                    essential = [st for st in recipe
                                 if _step_name(st) not in _BEAUTIFY_STEPS]
                    if len(essential) != len(recipe):
                        log(f"  Fast mode: {lname} ({kind}) — essential steps "
                            f"only ({_recipe_summary(essential) or 'none'}); "
                            f"cosmetic cleanup skipped")
                    recipe = essential
                if recipe:
                    log(f"  Recipe for {lname} ({kind}) [usage {ub}]: "
                        f"{_recipe_summary(recipe)}")
                    result_layer = _run_recipe(
                        gpkg_path, actual_table, lname, kind, lyr_name,
                        recipe, log, subprogress_fn=subprogress_fn
                    )
                    if subprogress_fn:
                        subprogress_fn(0, "")
                    if result_layer is None:
                        log(f"  ✗ Recipe produced nothing for {lname} ({kind})",
                            Qgis.Warning)
                        continue
                else:
                    uri = f"{gpkg_path}|layername={actual_table}"
                    result_layer = QgsVectorLayer(uri, lyr_name, "ogr")
                    if not result_layer.isValid():
                        log(f"  ✗ QgsVectorLayer invalid: {uri}", Qgis.Warning)
                        continue

                # 2) Record a manifest entry; the main thread re-opens this
                #    GeoPackage table, applies the (code, geom) style, and adds
                #    it. Style lookup is by (code, geom) ONLY — usage is
                #    irrelevant to symbology.
                manifest.append({
                    "sort_key": sort_key,
                    "table": actual_table,
                    "name": lyr_name,
                    "style_code": lname,
                    "geom": kind,
                })
                _out_crs = result_layer.crs()
                _oe = result_layer.extent()
                log(f"  ✓ '{lyr_name}' "
                    f"(sort {sort_key}, {result_layer.featureCount()} features)")
                log(f"    [crs] output {lname}/{kind}: "
                    f"{_out_crs.authid() or 'EMPTY'} "
                    f"(valid={_out_crs.isValid()}) "
                    f"extent=[{_oe.xMinimum():.2f},{_oe.yMinimum():.2f} .. "
                    f"{_oe.xMaximum():.2f},{_oe.yMaximum():.2f}]", Qgis.Info)

        progress(sel_idx + 1, total_layers)

    # ── Hybrid: ENC Coverage layer(s), per-usage — but ONLY for the usage
    #    bands the user actually selected coverage for (option (b)). Each gets
    #    its own coverage layer (harbour footprints separate from overview).
    if want_coverage and best_files:
        cov_by_usage = {}
        for dsnm, fpath in best_files.items():
            ub = usage_band(dsnm) or usage_band(
                os.path.splitext(os.path.basename(fpath))[0])
            if ub not in coverage_usages:
                continue   # user didn't pick ENC Coverage under this usage
            cov_by_usage.setdefault(ub, {})[dsnm] = fpath

        for ub in sorted(cov_by_usage.keys()):
            if cancel_fn is not None and cancel_fn():
                log("[load] cancel requested — stopping coverage build.",
                    Qgis.Warning)
                break
            try:
                cov_layer = _build_coverage_layer(
                    cov_by_usage[ub], dest_crs, gpkg_path, log,
                    usage=ub, subprogress_fn=subprogress_fn
                )
                if cov_layer is not None:
                    # Priority 90 within its usage → top of that usage's block.
                    sort_key = ub * 1000 + 90
                    manifest.append({
                        "sort_key": sort_key,
                        "table": _manifest_table(cov_layer),
                        "name": cov_layer.name(),
                        "style_code": "ENC_COVERAGE",
                        "geom": "area",
                    })
            except Exception as e:
                log(f"  ✗ Could not build ENC Coverage (usage {ub}): {e}",
                    Qgis.Warning)

    # Do NOT add layers here — this runs on the worker thread. Hand the manifest
    # back so the dock can register everything on the main thread (see
    # register_layers), which is the only safe place to add map layers and to
    # let an empty project adopt the load CRS.
    log(f"Built {len(manifest)} layer(s); handing to main thread to register.")
    log("Done.")
    return gpkg_path, manifest


def register_layers(gpkg_path, manifest, target_crs=None):
    """MAIN-THREAD registration of the layers built by load_layers().

    Adding map layers to the project must happen on the GUI thread. We re-open
    each finished GeoPackage table here, apply its (code, geom) style, and add
    it to the project root in priority order (ascending sort_key, so the
    highest-priority layer ends up on top). If the project is currently empty,
    we also stamp its CRS to the load CRS — replicating, on the main thread,
    QGIS's normal "empty project adopts the first layer's CRS" behaviour that an
    off-thread add would otherwise skip (the cause of the intermittent
    displaced-against-basemap bug).

    Returns the list of layers added.
    """
    proj = QgsProject.instance()
    cfg = get_config()
    was_empty = len(proj.mapLayers()) == 0

    ordered = sorted(manifest or [], key=lambda e: e.get("sort_key", 0))
    added = []
    for e in ordered:
        table = e.get("table")
        if not table:
            _log(f"  ✗ Manifest entry without a table: {e}", Qgis.Warning)
            continue
        uri = f"{gpkg_path}|layername={table}"
        layer = QgsVectorLayer(uri, e.get("name", table), "ogr")
        if not layer.isValid():
            _log(f"  ✗ Could not open {uri} to register", Qgis.Warning)
            continue
        style_code, geom = e.get("style_code"), e.get("geom")
        if style_code and geom:
            _apply_portrayal(layer, cfg, style_code, geom, _log)
        proj.addMapLayer(layer, True)   # standard, fully-wired add to root top
        added.append(layer)

    # Stamp the project CRS only when the project was empty, so we never
    # override a CRS the user (or a basemap loaded first) already established.
    if was_empty and added:
        crs = None
        if target_crs:
            cand = QgsCoordinateReferenceSystem(target_crs)
            if cand.isValid():
                crs = cand
        if crs is None:
            crs = added[0].crs()
        if crs.isValid():
            proj.setCrs(crs)
            _log(f"Empty project: set project CRS to {crs.authid()}.")

    _log(f"Registered {len(added)} layer(s) in the project root (main thread).")
    return added


def _step_name(step):
    """The algorithm name of a recipe step (string step or single-key dict)."""
    if isinstance(step, str):
        return step
    if isinstance(step, dict):
        for name in step:
            return name
    return ""


# Cosmetic geometry-cleaning steps: they remove cell seams and tidy polygons
# but a layer still portrays without them. The Render=Fast mode drops these and
# keeps every other (portrayal-essential) step — so SOUNDG's vertex split and
# zoom-level field still run in Fast (soundings need them to draw), while area
# dissolves are skipped. Render=Beautiful runs the whole recipe.
_BEAUTIFY_STEPS = {"fix", "dissolve", "singleparts"}


def _recipe_summary(recipe):
    """Short human string for a recipe, e.g. 'fix → dissolve(DRVAL1,DRVAL2) → singleparts'."""
    parts = []
    for step in recipe:
        if isinstance(step, str):
            parts.append(step)
        elif isinstance(step, dict):
            for name, params in step.items():
                if name == "dissolve" and isinstance(params, dict):
                    flds = ",".join(params.get("fields", []))
                    parts.append(f"dissolve({flds})")
                else:
                    parts.append(name)
    return " → ".join(parts)


def _run_recipe(gpkg_path, table_name, lname, kind, display_name, recipe, log,
                subprogress_fn=None):
    """Run a processing recipe on an already-merged GeoPackage table.

    Steps run in memory (TEMPORARY_OUTPUT); only the final result is written
    back to the GeoPackage (overwriting the table). Supported steps:
        "fix"                                  → native:fixgeometries
        {"dissolve": {"fields": [...]}}        → native:dissolve
        "singleparts"                          → native:multiparttosingleparts
    Returns a file-backed QgsVectorLayer, or None on failure.
    """
    try:
        import processing
    except Exception as e:
        log(f"    ✗ Processing framework unavailable: {e}", Qgis.Warning)
        return None

    try:
        uri = f"{gpkg_path}|layername={table_name}"
        current = QgsVectorLayer(uri, table_name, "ogr")
        if not current.isValid():
            log(f"    ✗ Could not open {uri} for recipe", Qgis.Warning)
            return None

        n_steps = max(1, len(recipe))
        for idx, step in enumerate(recipe):
            base = int(idx / n_steps * 90)
            span = int(90 / n_steps)

            if isinstance(step, str):
                name, params = step, {}
            elif isinstance(step, dict):
                name, params = next(iter(step.items()))
            else:
                continue

            if name == "fix":
                log("    → fix geometries")
                if subprogress_fn:
                    subprogress_fn(base, "Fixing geometries…")
                fb = _StageFeedback(subprogress_fn, base, span, "Fixing geometries…")
                # Quiet scope: if the LIST_AS_STRING option didn't fully
                # neutralise a map-typed attribute on this driver build, the
                # output write can still emit a per-feature CRITICAL about an
                # empty "map" value. The geometry result is unaffected, so we
                # suppress that noise; a genuine failure still yields no/empty
                # output, which is handled by the validity check after the loop.
                with _gdal_quiet():
                    current = processing.run(
                        "native:fixgeometries",
                        {"INPUT": current, "OUTPUT": "TEMPORARY_OUTPUT"},
                        feedback=fb,
                    )["OUTPUT"]

            elif name == "dissolve":
                field_names = [f.name() for f in current.fields()]
                want = params.get("fields", []) if isinstance(params, dict) else []
                use = [f for f in want if f in field_names]
                if use:
                    log(f"    → dissolve on {', '.join(use)}")
                else:
                    log("    → dissolve (requested fields absent — dissolving all)",
                        Qgis.Warning)
                if subprogress_fn:
                    subprogress_fn(base, "Dissolving…")
                fb = _StageFeedback(subprogress_fn, base, span, "Dissolving…")
                current = processing.run(
                    "native:dissolve",
                    {"INPUT": current, "FIELD": use, "OUTPUT": "TEMPORARY_OUTPUT"},
                    feedback=fb,
                )["OUTPUT"]

            elif name == "singleparts":
                log("    → multipart to singleparts")
                if subprogress_fn:
                    subprogress_fn(base, "Multipart → singlepart…")
                fb = _StageFeedback(subprogress_fn, base, span, "Multipart → singlepart…")
                current = processing.run(
                    "native:multiparttosingleparts",
                    {"INPUT": current, "OUTPUT": "TEMPORARY_OUTPUT"},
                    feedback=fb,
                )["OUTPUT"]

            elif name == "extractvertices":
                # Explode (multi)geometries into one point per vertex, each
                # carrying its own Z value plus vertex_index/part fields.
                # Essential for S-57 SOUNDG, which stores sounding clusters as
                # multipoint features with per-point depths in Z.
                log("    → extract vertices")
                if subprogress_fn:
                    subprogress_fn(base, "Extracting vertices…")
                fb = _StageFeedback(subprogress_fn, base, span, "Extracting vertices…")
                current = processing.run(
                    "native:extractvertices",
                    {"INPUT": current, "OUTPUT": "TEMPORARY_OUTPUT"},
                    feedback=fb,
                )["OUTPUT"]

            elif name == "zoomlevel":
                # Assign a per-point 'zoomLevel' for density-adaptive sounding
                # portrayal via greedy shoalest-first thinning. params is either
                # {radii: [...]} (fixed metres) or {spacing_multipliers, ...}
                # (adaptive: radii derived from median point spacing, and the
                # reveal-curve project variables are written). See
                # _assign_zoom_levels().
                zparams = params if isinstance(params, dict) else {}
                if subprogress_fn:
                    subprogress_fn(base, "Computing zoom levels…")
                current = _assign_zoom_levels(
                    current, zparams, log,
                    subprogress_fn=subprogress_fn, base=base, span=span,
                )

            else:
                log(f"    ? Unknown recipe step '{name}', skipping", Qgis.Warning)

            if current is None or not current.isValid():
                log(f"    ✗ Step '{name}' produced an invalid layer", Qgis.Warning)
                return None

        # Final and ONLY disk write: persist the processed result back into
        # the GeoPackage (overwriting the merged table for this geometry).
        log("    → writing recipe result to GeoPackage")
        if subprogress_fn:
            subprogress_fn(95, "Writing result…")

        transform_context = QgsProject.instance().transformContext()
        options = QgsVectorFileWriter.SaveVectorOptions()
        options.driverName = "GPKG"
        options.layerName = table_name
        options.fileEncoding = "UTF-8"
        options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer

        # Drop inherited "fid" so GDAL regenerates unique primary keys (dissolve
        # + singleparts can produce duplicate fids → UNIQUE constraint error).
        fields = current.fields()
        keep_attrs = [
            i for i in range(fields.count())
            if fields.at(i).name().lower() != "fid"
        ]
        options.attributes = keep_attrs
        options.layerOptions = ["FID=fid"]

        write_result = QgsVectorFileWriter.writeAsVectorFormatV3(
            current, gpkg_path, transform_context, options
        )
        if write_result[0] != QgsVectorFileWriter.NoError:
            log(f"    ✗ Recipe write failed: {write_result[1]}", Qgis.Warning)
            return None

        if subprogress_fn:
            subprogress_fn(100, "Loading…")

        out_uri = f"{gpkg_path}|layername={table_name}"
        out_layer = QgsVectorLayer(out_uri, display_name, "ogr")
        if not out_layer.isValid():
            log("    ✗ Reloaded recipe result invalid", Qgis.Warning)
            return None
        log(f"    ✓ processed: {out_layer.featureCount()} feature(s)")
        return out_layer

    except Exception as e:
        log(f"    ✗ Recipe failed for {lname} ({kind}): {e}", Qgis.Warning)
        return None


def _assign_zoom_levels(layer, params, log, subprogress_fn=None, base=0, span=10):
    """Add an integer 'zoomLevel' field to a point layer for density-adaptive
    sounding portrayal.

    params (dict) selects the radii:
      - {"radii": [...]}                       fixed radii in metres, OR
      - {"spacing_multipliers": [...],         adaptive: radii = multipliers x
         "start_factor": 1000,                 median nearest-neighbour spacing.
         "reveal_floor": 5000}                 Also writes project variables
                                               enc_soundg_start (= spacing x
                                               start_factor) and enc_soundg_ratio
                                               so all points reveal by 1:floor.

    Greedy shoalest-first thinning: points are split by sign of their Z
    (z >= 0 = depths/positives; z < 0 = drying heights/negatives). Within each
    pool, for each radius R (level L = 1, 2, 3, …):

        - sort the surviving (unassigned) points shoalest-first (ascending Z —
          for both pools, the smallest Z is the most dangerous: shallowest
          depth, or biggest drying height as the most-negative value),
        - walk that order keeping a spatial index of points ACCEPTED at this
          level; a candidate is accepted (zoomLevel = L) only if no
          already-accepted point lies within R; otherwise it survives to the
          next, smaller radius.

    Any point still unassigned after the last radius gets level len(radii)+1.
    Lower zoomLevel = more important = shown at coarser zoom. The shoalest
    point in any neighbourhood always survives to level 1, so hazards never
    hide behind decluttering. Returns a new in-memory QgsVectorLayer.
    """
    # Build output fields = input fields + zoomLevel (int).
    in_fields = layer.fields()
    out_fields = QgsFields(in_fields)
    if in_fields.indexFromName("zoomLevel") == -1:
        out_fields.append(QgsField("zoomLevel", QVariant.Int))

    # The greedy distance test runs in METRES. By load time the layer has
    # already been reprojected to the target (projected, metric) CRS, so we use
    # its native coordinates directly. As a defensive fallback, if the layer is
    # somehow still in a geographic CRS (degrees), transform to EPSG:3857 for
    # the distance math while keeping the ORIGINAL geometry for writing back.
    src_crs = layer.crs()
    if src_crs.isValid() and src_crs.isGeographic():
        metric_crs = QgsCoordinateReferenceSystem("EPSG:3857")
        to_metric = QgsCoordinateTransform(
            src_crs, metric_crs, QgsProject.instance()
        )
        log(f"    → zoomlevel: layer is geographic ({src_crs.authid()}); "
            "distances via EPSG:3857")
    else:
        to_metric = None
        log(f"    → zoomlevel: distances in native CRS "
            f"({src_crs.authid() or 'unknown'})")

    # Read every feature once: keep id, original geometry, z, attributes, and
    # a metric xy (EPSG:3857) used only for the distance/index math.
    pts = []   # list of dicts: {geom, z, attrs, xy(metric)}
    for feat in layer.getFeatures():
        geom = feat.geometry()
        if geom is None or geom.isEmpty():
            continue
        try:
            p = geom.constGet()
            z = p.z()
        except Exception:
            z = None
        if z is None or (isinstance(z, float) and (z != z)):  # None or NaN
            z = 0.0
        pt = geom.asPoint()
        mxy = QgsPointXY(pt.x(), pt.y())
        if to_metric is not None:
            try:
                mxy = to_metric.transform(mxy)
            except Exception:
                pass
        pts.append({
            "geom": geom,
            "z": float(z),
            "attrs": feat.attributes(),
            "xy": mxy,
        })

    total = len(pts)
    if total == 0:
        log("    ? zoomlevel: no points to process", Qgis.Warning)
        return layer

    # ── Determine radii: fixed, or adaptive from median point spacing. ──────
    adaptive = bool(params.get("spacing_multipliers"))
    reveal_floor = float(params.get("reveal_floor", 5000))
    start_factor = float(params.get("start_factor", 1000))
    spacing = None
    soundg_start = None

    if adaptive:
        mults = [float(m) for m in params["spacing_multipliers"]]
        spacing = _median_nn_distance(pts, log)
        if spacing is None or spacing <= 0:
            log("    ? zoomlevel: could not measure spacing, using 40 m",
                Qgis.Warning)
            spacing = 40.0
        radii = [m * spacing for m in mults]
        soundg_start = spacing * start_factor
        log(f"    → adaptive: median spacing {spacing:.1f} m → "
            f"radii {', '.join(str(int(r)) for r in radii)} m; "
            f"start 1:{int(soundg_start)}")
    else:
        try:
            radii = [float(r) for r in params.get("radii", [500, 250, 125])]
        except (TypeError, ValueError):
            log("    ? zoomlevel: bad radii, using [500,250,125]", Qgis.Warning)
            radii = [500, 250, 125]
        log(f"    → fixed radii {', '.join(str(int(r)) for r in radii)} m")

    # Split: z >= 0 positives (depths), z < 0 negatives (drying heights).
    positives = [p for p in pts if p["z"] >= 0.0]
    negatives = [p for p in pts if p["z"] < 0.0]
    final_level = len(radii) + 1

    def greedy(pool):
        """Greedy shoalest-first assignment. Mutates each point's 'level'."""
        for p in pool:
            p["level"] = None
        # Each pass works on the still-unassigned survivors.
        for lvl, radius in enumerate(radii, start=1):
            survivors = [p for p in pool if p["level"] is None]
            # Shoalest-first: ascending z (smallest = most dangerous).
            survivors.sort(key=lambda p: p["z"])
            index = QgsSpatialIndex()
            accepted = []
            next_id = 0
            for p in survivors:
                # Is there an accepted point within 'radius'?
                near = index.intersects(
                    _square_around(p["xy"], radius)
                )
                blocked = False
                for aid in near:
                    ap = accepted[aid]
                    if ap["xy"].distance(p["xy"]) <= radius:
                        blocked = True
                        break
                if not blocked:
                    p["level"] = lvl
                    f = QgsFeature()
                    f.setId(next_id)
                    f.setGeometry(QgsGeometry.fromPointXY(p["xy"]))
                    index.addFeature(f)
                    accepted.append(p)
                    next_id += 1
        # Remainder → final level.
        for p in pool:
            if p["level"] is None:
                p["level"] = final_level

    if subprogress_fn:
        subprogress_fn(base + int(span * 0.3), "Ranking soundings…")
    greedy(positives)
    if subprogress_fn:
        subprogress_fn(base + int(span * 0.6), "Ranking drying heights…")
    greedy(negatives)

    # Build the output memory layer.
    crs = layer.crs().authid()
    mem = QgsVectorLayer(f"PointZ?crs={crs}", layer.name(), "memory")
    dp = mem.dataProvider()
    dp.addAttributes([out_fields.at(i) for i in range(out_fields.count())])
    mem.updateFields()

    zoom_idx = out_fields.indexFromName("zoomLevel")
    out_feats = []
    for p in positives + negatives:
        f = QgsFeature(out_fields)
        attrs = list(p["attrs"])
        # Pad / set zoomLevel.
        if len(attrs) < out_fields.count():
            attrs += [None] * (out_fields.count() - len(attrs))
        attrs[zoom_idx] = int(p["level"])
        f.setAttributes(attrs)
        f.setGeometry(p["geom"])
        out_feats.append(f)
    dp.addFeatures(out_feats)
    mem.updateExtents()

    # Quick distribution log.
    dist = {}
    for p in positives + negatives:
        dist[p["level"]] = dist.get(p["level"], 0) + 1
    summary = ", ".join(f"L{k}:{dist[k]}" for k in sorted(dist))
    log(f"    ✓ zoomLevel assigned to {total} point(s) — {summary}")

    # ── Adaptive reveal: derive ratio so the HIGHEST level actually present
    #    lands exactly on the reveal floor (e.g. 1:5000), and write the curve
    #    parameters as project variables for SOUNDG.qml to consume. ──────────
    if adaptive and soundg_start:
        l_max = max(dist) if dist else 1
        if l_max <= 1:
            # Only one level → no curve; reveal everything at the floor.
            ratio = 1.0
            soundg_start = reveal_floor
        else:
            # start * ratio^(l_max-1) = reveal_floor  →  solve for ratio.
            ratio = (reveal_floor / soundg_start) ** (1.0 / (l_max - 1))
        try:
            proj = QgsProject.instance()
            QgsExpressionContextUtils.setProjectVariable(
                proj, "enc_soundg_start", round(float(soundg_start), 2))
            QgsExpressionContextUtils.setProjectVariable(
                proj, "enc_soundg_ratio", round(float(ratio), 5))
            log(f"    ✓ reveal: start 1:{int(soundg_start)}, ratio "
                f"{ratio:.4f}, L{l_max}→1:{int(reveal_floor)} "
                "(vars enc_soundg_start / enc_soundg_ratio)")
        except Exception as e:
            log(f"    ? could not write reveal variables: {e}", Qgis.Warning)

    if subprogress_fn:
        subprogress_fn(base + span, "")
    return mem


def _square_around(pt, r):
    """A square QgsRectangle of half-side r around a point, for index prefilter."""
    return QgsRectangle(pt.x() - r, pt.y() - r, pt.x() + r, pt.y() + r)


def _median_nn_distance(pts, log, sample_cap=20000):
    """Median nearest-neighbour distance (metres) over the points' metric xy.

    Builds a spatial index of all points, then for each point (optionally
    sampled, for very large sets) finds its nearest *other* point and records
    the distance. Returns the median, or None if it can't be computed.
    """
    n = len(pts)
    if n < 2:
        return None
    # Build an index over all points.
    index = QgsSpatialIndex()
    feats = []
    for i, p in enumerate(pts):
        f = QgsFeature()
        f.setId(i)
        f.setGeometry(QgsGeometry.fromPointXY(p["xy"]))
        index.addFeature(f)
        feats.append(f)

    # Sample for speed on huge sets; spacing is uniform enough that a sample
    # gives a stable median.
    step = max(1, n // sample_cap)
    dists = []
    for i in range(0, n, step):
        p = pts[i]
        # nearestNeighbor returns ids incl. self; ask for 2, drop self.
        ids = index.nearestNeighbor(p["xy"], 2)
        best = None
        for j in ids:
            if j == i:
                continue
            d = pts[j]["xy"].distance(p["xy"])
            if best is None or d < best:
                best = d
        if best is not None and best > 0:
            dists.append(best)
    if not dists:
        return None
    dists.sort()
    mid = len(dists) // 2
    if len(dists) % 2:
        return dists[mid]
    return 0.5 * (dists[mid - 1] + dists[mid])


def _qml_symbol_geometry(style_path):
    """Inspect a QML file and return the geometry kind its FEATURE RENDERER is
    authored for: 'point' | 'line' | 'area' | None (unknown).

    QGIS QML symbols declare a type attribute: <symbol type="marker|line|fill">
    (marker→point, line→line, fill→area). We must read only symbols belonging
    to the feature renderer (<renderer-v2>), NOT the whole file — a QML also
    contains symbols for elevation profiles, labeling leader lines, diagrams,
    etc., whose geometry says nothing about the feature portrayal. Reading the
    first symbol anywhere would false-positive on perfectly valid styles.

    Returns None if the renderer or its symbol type can't be determined, in
    which case the caller applies the style without blocking.
    """
    try:
        try:
            from .enc_xml_safe import parse_root
        except (ImportError, ValueError):
            from enc_xml_safe import parse_root
        root = parse_root(style_path)
        renderer = root.find("renderer-v2")
        if renderer is None:
            return None
        # First symbol declared within the renderer subtree.
        for sym in renderer.iter("symbol"):
            t = sym.get("type")
            if t == "marker":
                return "point"
            if t == "line":
                return "line"
            if t == "fill":
                return "area"
        return None
    except Exception:
        return None


def _apply_portrayal(layer, cfg, lname, kind, log):
    """Apply a geometry-specific QML portrayal style to a layer if one exists.

    Lookup is styles/<code>_<geom>.qml only (no bare-code fallback). Before
    applying, a geometry-compatibility guard checks the QML's symbol type
    matches the layer's geometry — so a style authored for one geometry can
    never be applied to another (which would render incorrectly or not at all).
    On mismatch the layer is left unstyled and a warning names the fix.
    """
    # Inland ENC notice marks: CEVNI sign-board portrayal built in code
    # (enc_notmrk) — head-group square zoomed out, SIGNI SVG + ORIENT rotation
    # zoomed in, same-position stack below 1:2500. QML fallback on failure.
    if (lname or "").upper() == "NOTMRK" and kind == "point":
        try:
            try:
                from . import enc_notmrk
            except (ImportError, ValueError):
                import enc_notmrk
            if enc_notmrk.apply(layer):
                log("  ✓ Applied CEVNI notice-mark portrayal")
                return
        except Exception as e:                          # noqa: BLE001
            log(f"  ⚠ CEVNI notice-mark portrayal failed ({e}); "
                "falling back to QML lookup", Qgis.Warning)

    style_path = cfg.style_path_for(lname, kind)
    if not style_path:
        return

    # Geometry-compatibility guard.
    qml_geom = _qml_symbol_geometry(style_path)
    if qml_geom is not None and qml_geom != kind:
        log(f"  ⚠ Style {os.path.basename(style_path)} is authored for "
            f"'{qml_geom}' geometry but the layer is '{kind}'. Left unstyled — "
            f"author a {lname}_{kind}.qml with a {kind}-appropriate symbol.",
            Qgis.Warning)
        return

    try:
        msg, ok = layer.loadNamedStyle(style_path)
        if ok:
            layer.triggerRepaint()
            log(f"  ✓ Applied style {os.path.basename(style_path)}")
        else:
            log(f"  ✗ Failed to load style {os.path.basename(style_path)}: {msg}",
                Qgis.Warning)
    except Exception as e:
        log(f"  ✗ Error applying style {os.path.basename(style_path)}: {e}",
            Qgis.Warning)


_S57_CAT = None


def _s57_catalogue():
    """Bundled S-57/Inland catalogue (code/attr maps), built once. GDAL-free."""
    global _S57_CAT
    if _S57_CAT is None:
        try:
            from .validation import Catalogue
        except (ImportError, ValueError):
            from validation import Catalogue
        _S57_CAT = Catalogue(os.path.join(os.path.dirname(__file__),
                                          "s57catalogue"))
    return _S57_CAT


_RS_CACHE = OrderedDict()
_RS_CACHE_MAX = 96


def _cached_recordset(fpath):
    """Parse a cell's merged record set ONCE and reuse it (LRU by path+mtime).
    Without this the native loader re-parses each cell for every object class
    (N classes x M cells full parses) — the cause of slow object scans / loads
    that don't seem to start. With it, a cell is parsed once per load."""
    try:
        from .s57 import S57RecordSet
    except (ImportError, ValueError):
        from s57 import S57RecordSet
    try:
        key = (fpath, os.path.getmtime(fpath))
    except OSError:
        key = (fpath, 0)
    rs = _RS_CACHE.get(key)
    if rs is not None:
        _RS_CACHE.move_to_end(key)
        return rs
    rs = S57RecordSet.from_paths([fpath])
    _RS_CACHE[key] = rs
    _RS_CACHE.move_to_end(key)
    while len(_RS_CACHE) > _RS_CACHE_MAX:
        _RS_CACHE.popitem(last=False)
    return rs


def _qgs_geom(g):
    """Build a QgsGeometry DIRECTLY from a native s57 geometry value — no WKT
    round-trip (serialising to WKT and re-parsing dominated load time). Z is
    preserved for points/soundings."""
    k = g.kind
    if k == "POINT":
        if g.z is not None:
            return QgsGeometry(QgsPoint(g.x, g.y, g.z))
        return QgsGeometry.fromPointXY(QgsPointXY(g.x, g.y))
    if k == "MULTIPOINT":
        if g.z_present:
            mp = QgsMultiPoint()
            for p in g.points:
                mp.addGeometry(QgsPoint(p[0], p[1], p[2] if len(p) > 2 else 0.0))
            return QgsGeometry(mp)
        return QgsGeometry.fromMultiPointXY(
            [QgsPointXY(p[0], p[1]) for p in g.points])
    if k == "LINESTRING":
        return QgsGeometry.fromPolylineXY(
            [QgsPointXY(x, y) for x, y in g.points])
    if k == "MULTILINESTRING":
        return QgsGeometry.fromMultiPolylineXY(
            [[QgsPointXY(x, y) for x, y in part] for part in g.parts])
    if k == "POLYGON":
        return QgsGeometry.fromPolygonXY(
            [[QgsPointXY(x, y) for x, y in ring] for ring in g.rings])
    if k == "MULTIPOLYGON":
        return QgsGeometry.fromMultiPolygonXY(
            [[[QgsPointXY(x, y) for x, y in ring] for ring in poly]
             for poly in g.polygons])
    return QgsGeometry.fromWkt(g.wkt())   # safety fallback


def _lnam(foid):
    """S-57 LNAM string in GDAL's format: AGEN(4 hex)+FIDN(8 hex)+FIDS(4 hex)."""
    if not foid:
        return None
    agen, fidn, fids = (foid + (0, 0, 0))[:3]
    return f"{agen or 0:04X}{fidn or 0:08X}{fids or 0:04X}"


def _lnam_from_ffpt(v):
    """An FFPT reference LNAM → the same hex string as a feature's LNAM, so
    buoy↔light links compare equal. The native reader decodes the FFPT LNAM as a
    little-endian 8-byte int (AGEN[2] + FIDN[4] + FIDS[2] in file order), so the
    fields are unpacked from the low end. A value already a string is returned
    upper-cased."""
    if v is None:
        return None
    if isinstance(v, str):
        return v.strip().upper() or None
    try:
        n = int(v)
    except (TypeError, ValueError):
        return None
    agen = n & 0xFFFF
    fidn = (n >> 16) & 0xFFFFFFFF
    fids = (n >> 48) & 0xFFFF
    return f"{agen:04X}{fidn:08X}{fids:04X}"


# Native geometry .kind → (QGIS memory geom-type, _feature_geom_kind value).
_KIND_MEM = {"POINT": ("Point", "point"), "MULTIPOINT": ("MultiPointZ", "point"),
             "LINESTRING": ("MultiLineString", "line"),
             "MULTILINESTRING": ("MultiLineString", "line"),
             "POLYGON": ("MultiPolygon", "area"),
             "MULTIPOLYGON": ("MultiPolygon", "area")}


def _class_attr_union(usage_files, layer_name):
    """Complete, ordered attribute-acronym field set for one object class —
    sourced from the CATALOGUE (the full datamodel), not just from data content.

    = the catalogue's permissible attributes for the class
    (``Attribute_A ∪ B ∪ C`` — the same set the Auditor uses for FC conformance)
    UNION any attribute actually present in the band's cells (a safety net for
    Inland extensions or attributes not listed in the bundled catalogue).

    Driving the schema from the catalogue means every loaded layer carries its
    full, predictable datamodel (e.g. LIGHTS always has SECTR1/SECTR2) regardless
    of which cell is written first or whether a given attribute happens to be
    populated anywhere. It also makes the per-band GeoPackage merge lossless: the
    merge writes every cell into ONE table (first cell creates the schema, later
    cells append with ``AppendToLayerNoNewFields``, which silently drops any field
    the first cell lacked), so a shared complete schema across all cells is what
    keeps a later cell's SECTR1/SECTR2 from being discarded.

    Order: permissible attrs first (sorted), then any present-but-not-permissible
    extras (sorted), for determinism."""
    cat = _s57_catalogue()
    c2a, a2acr = cat.code_to_acronym, cat.attr_code_to_acronym
    # Permissible set per the bundled S-57/Inland object catalogue. Match the
    # exact acronym first (Inland lowercase classes have their own rows); fall
    # back to the upper-cased key for the standard classes.
    permissible = cat.permissible_attrs.get(layer_name)
    if not permissible:
        permissible = cat.permissible_attrs.get(layer_name.upper(), frozenset())
    present = set()
    for fp in usage_files:
        try:
            rs = _cached_recordset(fp)
        except Exception:                                 # noqa: BLE001
            continue
        for cell in rs.cells:
            if getattr(cell, "parse_error", None):
                continue
            for f in cell.features:
                if c2a.get(f.objl) == layer_name:
                    for code in (f.attrs or {}):
                        if code in a2acr:
                            present.add(a2acr[code])
    perm = sorted(permissible)
    extras = sorted(present - set(perm))
    return perm + extras


def _native_class_layer(fpath, layer_name, attr_acrs=None):
    """Native (GDAL-free) in-memory QgsVectorLayer for one S-57 object class — a
    faithful drop-in for GDAL's per-class S-57 layer. Carries native geometry +
    identity fields (RCID/PRIM/OBJL/RVER/AGEN/FIDN/FIDS/LNAM/LNAM_REFS) + the
    class's S-57 attributes by acronym, in EPSG:4326 (as GDAL's S-57 layers are).

    ``attr_acrs`` (optional) is a caller-supplied ordered list of attribute
    acronyms to use as the attribute field set — the shared superset across all
    cells of a band (see _class_attr_union), so the per-band GeoPackage merge is
    lossless. When None, the field set is this cell's own attributes (legacy
    single-cell behaviour). Attributes absent from a given cell are written NULL.

    Degenerate/unbuildable features are added with NULL geometry so the loader's
    held-back partition catches them exactly as it does with GDAL. Returns None
    if the class is absent or the cell can't be read (caller falls back to GDAL)."""
    try:
        from .s57 import GeometryBuilder
    except (ImportError, ValueError):
        from s57 import GeometryBuilder
    cat = _s57_catalogue()
    c2a, a2acr = cat.code_to_acronym, cat.attr_code_to_acronym
    # Let a genuine read failure PROPAGATE (so the caller can fall back to GDAL);
    # a successful read with no matching class returns None (class simply absent —
    # the caller must NOT then invoke GDAL for it). Match the EXACT acronym (not
    # case-insensitively) so a cell holding both standard `BOYLAT` (OBJL 17) and
    # Inland `boylat` (OBJL 17029) keeps them as two separate layers, as GDAL does.
    cells = [c for c in _cached_recordset(fpath).cells
             if not getattr(c, "parse_error", None)]
    if not cells:
        raise RuntimeError(f"native reader: no readable cell in {fpath}")
    cell = cells[0]
    feats = [f for f in cell.features if c2a.get(f.objl) == layer_name]
    if not feats:
        return None
    builder = GeometryBuilder(cell)

    # Build geometry per feature; pick the dominant kind for the layer type.
    geoms, kinds = {}, {}
    for f in feats:
        g = None
        try:
            g = builder.geometry(f)
        except Exception:                                 # noqa: BLE001
            g = None
        geoms[id(f)] = g
        if g is not None:
            kinds[g.kind] = kinds.get(g.kind, 0) + 1
    if not kinds:
        dom = "POLYGON" if any(f.prim == 3 for f in feats) else (
            "LINESTRING" if any(f.prim == 2 for f in feats) else "POINT")
    else:
        dom = max(kinds, key=kinds.get)
    mem_type, _kind = _KIND_MEM.get(dom, ("MultiPolygon", "area"))

    # Field set: identity + attribute acronyms. Prefer the caller's shared
    # superset (lossless band merge); else fall back to this cell's own union.
    if attr_acrs is None:
        attr_acrs = []
        seen = set()
        for f in feats:
            for code in (f.attrs or {}):
                acr = a2acr.get(code)
                if acr and acr not in seen:
                    seen.add(acr)
                    attr_acrs.append(acr)
    layer = QgsVectorLayer(f"{mem_type}?crs=EPSG:4326", layer_name, "memory")
    dp = layer.dataProvider()
    id_fields = [("RCID", QVariant.LongLong), ("PRIM", QVariant.Int),
                 ("OBJL", QVariant.Int), ("RVER", QVariant.Int),
                 ("AGEN", QVariant.Int), ("FIDN", QVariant.LongLong),
                 ("FIDS", QVariant.Int), ("LNAM", QVariant.String),
                 ("LNAM_REFS", QVariant.String)]
    dp.addAttributes([QgsField(n, t) for n, t in id_fields]
                     + [QgsField(a, QVariant.String) for a in attr_acrs])
    layer.updateFields()

    out = []
    for f in feats:
        feat = QgsFeature(layer.fields())
        g = geoms.get(id(f))
        if g is not None and g.kind == dom:
            feat.setGeometry(_qgs_geom(g))
        # else: leave NULL geometry → held-back partition flags it (as GDAL does)
        agen, fidn, fids = ((f.foid or (None, None, None)) + (None, None, None))[:3]
        vals = {"RCID": f.rcid, "PRIM": f.prim, "OBJL": f.objl, "RVER": f.rver,
                "AGEN": agen, "FIDN": fidn, "FIDS": fids, "LNAM": _lnam(f.foid),
                "LNAM_REFS": ",".join(
                    s for s in (_lnam_from_ffpt(r[0]) for r in (f.ffpt or []))
                    if s)}
        for code, v in (f.attrs or {}).items():
            acr = a2acr.get(code)
            if acr:
                vals[acr] = None if v is None else str(v)
        feat.setAttributes([vals.get(fld.name()) for fld in layer.fields()])
        out.append(feat)
    dp.addFeatures(out)
    layer.updateExtents()
    return layer if layer.isValid() else None


def _open_ogr_layer(fpath, layer_name, attr_acrs=None):
    """One S-57 object class as a QgsVectorLayer — fully native (in-house ISO 8211
    engine + bundled catalogue), no GDAL S-57 driver. ``attr_acrs`` forwards the
    shared per-band attribute superset for a lossless merge (see
    _class_attr_union). Returns None if the class is absent or the cell can't be
    read (verified render-correct in QGIS, v1.32.x)."""
    try:
        return _native_class_layer(fpath, layer_name, attr_acrs=attr_acrs)
    except Exception as e:                                # noqa: BLE001
        _log(f"  ✗ native read failed for {os.path.basename(fpath)} "
             f"[{layer_name}]: {e}", Qgis.Warning)
        return None


def _iter_ogr_layers_ci(fpath, names_upper):
    """[(actual_acronym, QgsVectorLayer), …] for every object class in the cell
    matching one of names_upper (UPPERCASE), fully native. A cell holding BOTH
    standard `BOYLAT` (OBJL 17) and Inland `boylat` (OBJL 17029) yields both as
    separate layers (distinct acronyms), so co-located topmarks aren't orphaned."""
    out = []
    try:
        c2a = _s57_catalogue().code_to_acronym
        present, seen = [], set()
        for cell in _cached_recordset(fpath).cells:
            if getattr(cell, "parse_error", None):
                continue
            for f in cell.features:
                acr = c2a.get(f.objl)
                if acr and acr.upper() in names_upper and acr not in seen:
                    seen.add(acr)
                    present.append(acr)
        for acr in present:
            lyr = _native_class_layer(fpath, acr)
            if lyr is not None and lyr.featureCount() > 0:
                out.append((acr, lyr))
    except Exception as e:                                # noqa: BLE001
        _log(f"  ✗ native read failed for {os.path.basename(fpath)}: {e}",
             Qgis.Warning)
    return out


def _list_gpkg_tables(gpkg_path):
    """Layer (table) names in a GeoPackage.

    A GeoPackage IS a SQLite database, so we read the layer names straight from
    its ``gpkg_contents`` registry with stdlib ``sqlite3`` (GDAL-free, like the
    rest of the native path). This replaces a per-call
    ``QgsProviderRegistry.querySublayers()`` scan — which opens and probes every
    sublayer's geometry type/extent and, called once per written table during a
    load, dominated load time (≈3.8 s / ~19 % on a 68-table set). The sqlite read
    returns the identical set of table names in ~0.1 ms. The provider scan is
    kept as a fallback for any file without a readable ``gpkg_contents``.
    """
    if not os.path.exists(gpkg_path):
        return []
    import sqlite3
    try:
        con = sqlite3.connect(gpkg_path)
        try:
            rows = con.execute("SELECT table_name FROM gpkg_contents").fetchall()
        finally:
            con.close()
        return [r[0] for r in rows if r and r[0]]
    except Exception:                                     # noqa: BLE001
        try:
            from qgis.core import QgsProviderRegistry
            return [s.name() for s in
                    QgsProviderRegistry.instance().querySublayers(gpkg_path)]
        except Exception:                                 # noqa: BLE001
            return []


def _find_gpkg_table(gpkg_path, desired_name):
    tables = _list_gpkg_tables(gpkg_path)
    if desired_name in tables:
        return desired_name
    lower = desired_name.lower()
    for t in tables:
        if t.lower() == lower:
            return t
    sanitised = lower.replace(" ", "_").replace("-", "_")
    for t in tables:
        if t.lower().replace(" ", "_").replace("-", "_") == sanitised:
            return t
    return None

# ── ENC Coverage hybrid layer ────────────────────────────────────────────────


# DSID metadata fields carried by the coverage layer, in display order.
_COVERAGE_DSID_FIELDS = [
    "EDTN", "UPDN", "ISDT", "DSNM", "PRSP", "PRED", "AGEN", "CSCL", "COMT",
]


# ── Merged Aids-to-Navigation layer ──────────────────────────────────────────
# The S-57 buoy classes whose floating bodies merge into one AtoN layer. (Beacon
# BCNxxx classes are deliberately held out for a later round.)
_ATON_BUOY_CLASSES = ["BOYLAT", "BOYCAR", "BOYSAW", "BOYSPP"]
# Fixed beacons. Aids to navigation just like buoys; they carry topmarks and
# lights through the same LNAM_REFS relationship and co-location, so once they
# are gathered the existing topmark/light passes fold onto them unchanged. They
# go into the SAME hybrid layer — S57CLASS (BCNLAT/…/BCNISD) distinguishes a
# beacon row from a buoy row so portrayal can symbolise each differently.
_ATON_BEACON_CLASSES = ["BCNLAT", "BCNCAR", "BCNSAW", "BCNSPP", "BCNISD"]
# Every floating-or-fixed body class the merge gathers, as an UPPERCASE set.
_ATON_BODY_CLASSES_UPPER = {c.upper()
                            for c in _ATON_BUOY_CLASSES + _ATON_BEACON_CLASSES}
_ATON_TOPMARK_CLASS = "TOPMAR"

# Body attributes carried onto each row (union across buoy AND beacon classes —
# buoys carry BOYSHP, beacons BCNSHP; CATLAM is on BOY/BCNLAT, CATCAM on
# BOY/BCNCAR, CATSPM on BOY/BCNSPP; absent ones stay null). The co-located
# topmark's own attributes are copied onto the SAME row under TT_-prefixed
# names so they don't collide with the body's COLOUR/COLPAT.
_ATON_BUOY_FIELDS = ["BOYSHP", "BCNSHP", "COLOUR", "COLPAT", "CATLAM", "CATCAM",
                     "CATSPM", "MARSYS", "OBJNAM"]
_ATON_TOPMARK_FIELDS = ["TOPSHP", "COLOUR", "COLPAT"]  # read from TOPMAR
# How TOPMAR's fields are renamed on the merged row (avoids clobbering the
# buoy's own COLOUR / COLPAT). The QML's topmark symbol layer filters on these.
_ATON_TT_RENAME = {"TOPSHP": "TT_TOPSHP", "COLOUR": "TT_COLOUR",
                   "COLPAT": "TT_COLPAT"}
# S-57 identity / relationship keys carried onto every merged row. LNAM is the
# feature's unique long-name; LNAM_REFS holds its master/slave links (e.g. a
# buoy pointing at its topmark and light). Carried through so downstream work —
# notably AtoN-vs-WFS validation — can join on the explicit S-57 relationship
# rather than on coordinates.
_ATON_IDENT_FIELDS = ["LNAM", "LNAM_REFS"]

# ── Lit buoys (v55). A LIGHTS feature that belongs to a buoy — the buoy's
# LNAM_REFS contains the light's LNAM AND the two share a position — has its
# light attributes folded onto the buoy's merged row under LT_-prefixed names,
# so one row describes the lit mark (body + topmark + light) and can render an
# "Fl G 3s"-style label. Every light still stays in the standalone LIGHTS layer;
# there it gets a boolean `buoyLight` flag (True for the ones folded here) so
# portrayal can hide buoy-lights that are already drawn as lit buoys, while
# shore / sector / beacon lights keep drawing.
_ATON_LIGHT_CLASS = "LIGHTS"
_ATON_LIGHT_FIELDS = ["COLOUR", "LITCHR", "SIGGRP", "SIGPER", "VALNMR", "CATLIT"]
_ATON_LT_RENAME = {f: "LT_" + f for f in _ATON_LIGHT_FIELDS}
_BUOYLIGHT_FIELD = "buoyLight"

# Strict per-column datatypes for the merged AtoN layer — the SAME production
# data-model policy applied to the S-101 AtoN hybrid (see s101.qgis_loader
# _ATON_FIELDDEFS), so both AtoN layers are mappable in the Source Validator:
# single-valued S-57 enumerations are INTEGER sized to their code domain, the
# light period/range are REAL, and STRING is kept only where the value is text —
# names, identity keys (LNAM/LNAM_REFS), the light signal group, and the
# MULTI-VALUED colour attributes (COLOUR/TT_COLOUR/LT_COLOUR — GDAL returns e.g.
# '3,1', also the form AtoN_point.qml's rules match on). Unlisted names → String.
_ATON_FIELD_TYPES = {
    "BOYSHP": (QVariant.Int, 1, 0), "BCNSHP": (QVariant.Int, 1, 0),
    "COLPAT": (QVariant.Int, 1, 0), "CATLAM": (QVariant.Int, 1, 0),
    "CATCAM": (QVariant.Int, 1, 0), "CATSPM": (QVariant.Int, 2, 0),
    "MARSYS": (QVariant.Int, 2, 0),
    "TT_TOPSHP": (QVariant.Int, 2, 0), "TT_COLPAT": (QVariant.Int, 1, 0),
    "LT_LITCHR": (QVariant.Int, 2, 0), "LT_CATLIT": (QVariant.Int, 2, 0),
    "LT_SIGPER": (QVariant.Double, 6, 1), "LT_VALNMR": (QVariant.Double, 6, 1),
}


def _aton_field(name):
    """A QgsField for an AtoN column at its strict type (default String)."""
    typ, ln, pr = _ATON_FIELD_TYPES.get(name, (QVariant.String, 0, 0))
    return QgsField(name, typ, "", ln, pr)


def _aton_coerce(name, value):
    """Coerce a raw attribute value to its column's declared type — int/real, or
    string. Returns None (→ NULL, not '' and not 0) when the value is missing or
    can't be parsed for a numeric column."""
    typ = _ATON_FIELD_TYPES.get(name, (QVariant.String, 0, 0))[0]
    if value is None:
        return None
    if typ == QVariant.Int:
        try:
            return int(str(value).strip())
        except (TypeError, ValueError):
            return None
    if typ == QVariant.Double:
        try:
            return float(str(value).strip())
        except (TypeError, ValueError):
            return None
    return value


# Rounding (decimal degrees) used to pair a topmark with its buoy by shared
# position. The data has them at identical x,y; rounding only guards against
# floating-point noise. ~6 dp ≈ 0.1 m, far tighter than buoy spacing.
_ATON_PAIR_DP = 6


def _aton_attr(feat, name):
    """Read an attribute by name from a feature, or None if absent/blank.

    The lookup is CASE-INSENSITIVE. Inland-ENC reuses the standard attribute
    acronyms but encodes some of them in lowercase — e.g. the inland `boylat`
    class carries `BOYSHP`/`COLOUR`/`COLPAT` (upper) but `catlam`/`marsys`
    (lower). Trying the exact name first and falling back to a case-insensitive
    match lets one field list (`CATLAM`) pull the value from either spelling,
    so the normalised uppercase output field is populated for standard AND
    Inland buoys alike and the existing QML keeps working untouched.
    """
    try:
        flds = feat.fields()
        idx = flds.indexFromName(name)
        if idx < 0:
            idx = flds.lookupField(name)   # QgsFields.lookupField is case-insensitive
        if idx < 0:
            return None
        v = feat.attribute(idx)
        if v is None:
            return None
        # Normalise QVariant-null and empty strings to None.
        s = str(v).strip()
        return v if s != "" and s.upper() != "NULL" else None
    except Exception:
        return None


def _aton_pos_key(geom):
    """A hashable rounded (x, y) key for a point geometry, or None."""
    try:
        p = geom.asPoint()
        return (round(p.x(), _ATON_PAIR_DP), round(p.y(), _ATON_PAIR_DP))
    except Exception:
        return None


def _norm_lnam(x):
    """Normalise an LNAM to a comparable token (strip quotes/brackets, upper)."""
    if x is None:
        return ""
    return str(x).strip().strip("'\"[] ").upper()


def _parse_lnam_refs(v):
    """LNAM_REFS -> set of normalised LNAM tokens.

    GDAL may expose the references as a Python list (StringList) or as a single
    comma/space-joined string, so handle both.
    """
    if v is None:
        return set()
    if isinstance(v, (list, tuple)):
        items = [str(x) for x in v]
    else:
        items = re.split(r"[,\s]+", str(v).strip().strip("[]"))
    out = set()
    for it in items:
        n = _norm_lnam(it)
        if n:
            out.add(n)
    return out


# ── Reader-agnostic AtoN "smelting" core ─────────────────────────────────────
# These two functions hold the actual fold ALGORITHM, independent of how the
# features were read. The main project loader feeds them OGR-read features; the
# Chart Viewer feeds them native-engine features. Same engine, two front-ends —
# so a buoy, its topmark and its light are coupled identically in both.
def match_buoy_lights(refs_by_pos, lights):
    """Which lights belong to a buoy: co-located with a buoy AND named in that
    buoy's LNAM_REFS (relationship confirmed by co-location — the conservative
    rule). ``lights`` is [(lnam, pos_key, {LT_*: val})]; the caller must have
    already excluded sectored lights (they keep their own LIGHTS row so the
    sector arcs survive). Returns (claimed_lnams:set, fold:{pos_key: {LT_*}}).
    Lights not matched here are simply absent from the result — they stay in the
    LIGHTS layer, fully portrayed. """
    claimed, fold = set(), {}
    for lnam, pos, lt in lights:
        host_refs = refs_by_pos.get(pos)
        if host_refs and lnam in host_refs:
            claimed.add(lnam)
            fold[pos] = lt
    return claimed, fold


def compose_aton_rows(buoy_rows, topmark_rows, light_fold):
    """Fold co-located topmarks (by position) and pre-matched buoy-lights onto
    buoy/beacon body rows; topmarks with no body partner stay as standalone rows
    so nothing is dropped. Mutates each host row's ``attrs`` in place.

    buoy_rows    : [{"geom", "attrs", "pos"}]  attrs carry S57CLASS + idents
    topmark_rows : [{"geom", "pos", "tt", "standalone_attrs"}]
                   tt = {TT_*: val} to fold onto a host; standalone_attrs the
                   full attr dict to use when the topmark has no host buoy.
    light_fold   : {pos_key: {LT_*: val}}  already matched (see match_buoy_lights)

    Returns (rows:[{"geom","attrs"}], stats:{"paired","lit","standalone"}).
    """
    buoys_by_pos = {}
    for i, r in enumerate(buoy_rows):
        if r.get("pos") is not None:
            buoys_by_pos[r["pos"]] = i            # last wins (matches loader)
    n_paired = 0
    standalone = []
    for tr in topmark_rows:
        pos = tr.get("pos")
        host = buoys_by_pos.get(pos) if pos is not None else None
        if host is not None:
            hattrs = buoy_rows[host]["attrs"]
            if hattrs.get("TT_TOPSHP") is None and \
                    hattrs.get("TT_COLOUR") is None:   # first topmark wins
                hattrs.update(tr["tt"])
                n_paired += 1
        else:
            standalone.append({"geom": tr["geom"],
                               "attrs": dict(tr["standalone_attrs"])})
    n_lit = 0
    for pos, lt in (light_fold or {}).items():
        host = buoys_by_pos.get(pos)
        if host is None:
            continue
        hattrs = buoy_rows[host]["attrs"]
        wrote = False
        for k, v in lt.items():
            if v is not None and hattrs.get(k) is None:
                hattrs[k] = v
                wrote = True
        if wrote:
            n_lit += 1
    rows = [{"geom": r["geom"], "attrs": r["attrs"]} for r in buoy_rows]
    rows += standalone
    return rows, {"paired": n_paired, "lit": n_lit,
                  "standalone": len(standalone)}


def _compute_buoy_light_links(layer_index, selection, log):
    """Pre-pass: find every LIGHTS feature that belongs to a buoy.

    A light is a buoy-light when the buoy's LNAM_REFS contains the light's LNAM
    *and* the two share a position (the conservative "relationship confirmed by
    co-location" rule). Returns::

        {"claimed_lnams": set(<light LNAM>, ...),
         "fold_by_usage": {ub: {buoy_pos_key: {LT_field: value, ...}}}}

    Computed once, before any layer is built, so the AtoN merge (which folds the
    light onto the buoy row) and the LIGHTS build (which tags `buoyLight`) draw
    from the SAME result regardless of the order those selections are processed.
    Only usage bands where AtoN or LIGHTS is selected are scanned.
    """
    buoy_fam = _ATON_BODY_CLASSES_UPPER
    relevant_ubs = {ub for (ub, ln) in selection
                    if ln == "AtoN" or ln.upper() == _ATON_LIGHT_CLASS}
    claimed = set()
    fold_by_usage = {}
    if not relevant_ubs:
        return {"claimed_lnams": claimed, "fold_by_usage": fold_by_usage}

    for ub in relevant_ubs:
        # Cells at this usage that hold buoys and/or lights.
        files = set()
        for key, fs in layer_index.items():
            ku = key.upper()
            if ku in buoy_fam or ku == _ATON_LIGHT_CLASS:
                for f in fs:
                    if usage_band(os.path.splitext(os.path.basename(f))[0]) == ub:
                        files.add(f)

        refs_by_pos = {}     # buoy pos_key -> set(referenced LNAMs) — co-location
        lights = []          # (lnam, pos_key, {LT_ field: value})
        for f in files:
            for nm, src in _iter_ogr_layers_ci(f, buoy_fam):
                for feat in src.getFeatures():
                    g = feat.geometry()
                    if g is None or g.isEmpty():
                        continue
                    pos = _aton_pos_key(g)
                    if pos is None:
                        continue
                    refs = _parse_lnam_refs(_aton_attr(feat, "LNAM_REFS"))
                    if refs:
                        refs_by_pos.setdefault(pos, set()).update(refs)
            for nm, src in _iter_ogr_layers_ci(f, {_ATON_LIGHT_CLASS}):
                for feat in src.getFeatures():
                    g = feat.geometry()
                    if g is None or g.isEmpty():
                        continue
                    if QgsWkbTypes.geometryType(g.wkbType()) != \
                            QgsWkbTypes.PointGeometry:
                        continue
                    pos = _aton_pos_key(g)
                    lnam = _norm_lnam(_aton_attr(feat, "LNAM"))
                    if not lnam or pos is None:
                        continue
                    # SECTORED lights are never folded/claimed: the merged AtoN
                    # row carries only LT_COLOUR/LITCHR/SIGGRP/SIGPER/VALNMR — it
                    # cannot represent sector arcs (SECTR1/SECTR2). Folding one
                    # (e.g. a fixed sector light on a beacon, which IS in the AtoN
                    # body family) would tag it buoyLight=True and the LIGHTS
                    # portrayal would hide it, killing the sector display. So a
                    # light with any sector limit stays in the LIGHTS layer.
                    if (_aton_attr(feat, "SECTR1") is not None
                            or _aton_attr(feat, "SECTR2") is not None):
                        continue
                    lt = {}
                    for lf in _ATON_LIGHT_FIELDS:
                        lt[_ATON_LT_RENAME[lf]] = _aton_attr(feat, lf)
                    lights.append((lnam, pos, lt))

        claimed_u, fold = match_buoy_lights(refs_by_pos, lights)
        claimed |= claimed_u
        fold_by_usage[ub] = fold
        log(f"    [lights] U{ub}: {len(lights)} light(s), "
            f"{len(fold)} on buoys", Qgis.Info)

    log(f"[lights] buoy-lights matched: {len(claimed)} across "
        f"{len(relevant_ubs)} usage band(s)", Qgis.Info)
    return {"claimed_lnams": claimed, "fold_by_usage": fold_by_usage}


def _tag_buoylights(gpkg_path, table, claimed_lnams, log):
    """Add a boolean `buoyLight` field to a written LIGHTS GeoPackage table.

    True for lights whose LNAM is in claimed_lnams (the buoy-lights folded into
    the hybrid), False for all others. No light is moved or dropped — portrayal
    uses the flag to hide buoy-lights already shown as lit buoys while keeping
    shore / sector / beacon lights drawing.
    """
    try:
        actual = _find_gpkg_table(gpkg_path, table)
        if actual is None:
            return
        lyr = QgsVectorLayer(f"{gpkg_path}|layername={actual}", actual, "ogr")
        if not lyr.isValid():
            return
        lnam_idx = lyr.fields().lookupField("LNAM")   # case-insensitive
        if lnam_idx < 0:
            log(f"    [lights] {actual}: no LNAM field — cannot tag "
                f"{_BUOYLIGHT_FIELD}", Qgis.Warning)
            return
        if lyr.fields().lookupField(_BUOYLIGHT_FIELD) < 0:
            if not lyr.dataProvider().addAttributes(
                    [QgsField(_BUOYLIGHT_FIELD, QVariant.Bool)]):
                log(f"    [lights] {actual}: could not add "
                    f"{_BUOYLIGHT_FIELD}", Qgis.Warning)
                return
            lyr.updateFields()
        bl_idx = lyr.fields().lookupField(_BUOYLIGHT_FIELD)
        lyr.startEditing()
        n_true = 0
        for feat in lyr.getFeatures():
            is_bl = _norm_lnam(feat.attribute(lnam_idx)) in claimed_lnams
            if is_bl:
                n_true += 1
            lyr.changeAttributeValue(feat.id(), bl_idx, bool(is_bl))
        lyr.commitChanges()
        log(f"    [lights] {actual}: {_BUOYLIGHT_FIELD}=True on {n_true} of "
            f"{lyr.featureCount()} light(s)", Qgis.Info)
    except Exception as e:
        log(f"    [lights] tag failed for {table}: {e}", Qgis.Warning)


def _build_aton_layer(usage_files, ub, dest_crs, gpkg_path, log,
                      base_priority=80, subprogress_fn=None, light_fold=None):
    """Build ONE merged Aids-to-Navigation point layer for a usage band.

    Floating buoy bodies (BOYLAT/BOYCAR/BOYSAW/BOYSPP) and fixed beacon bodies
    (BCNLAT/BCNCAR/BCNSAW/BCNSPP/BCNISD) become rows carrying their attributes,
    with S57CLASS recording each row's class. Each co-located TOPMAR feature
    (matched by identical point position) has its TOPSHP/COLOUR/COLPAT copied
    onto the body's row under TT_-prefixed names, so one row describes the whole
    mark — enabling a combined maptip and a topmark symbol whose vertical offset
    can key off the body's BOYSHP/BCNSHP. Topmarks with no body partner are kept
    as their own rows so no data is silently dropped.

    `light_fold` (from _compute_buoy_light_links) is {buoy_pos_key: {LT_field:
    value}} for lights that belong to a buoy; each is folded onto the buoy at
    that position under LT_-prefixed names, so a lit buoy carries its light's
    colour / character / period / range for an "Fl G 3s"-style label.

    Returns a styled, file-backed QgsVectorLayer, or None if nothing was built.
    """
    # Schema: class tag + buoy fields + renamed topmark fields + light fields
    # + identity keys + SRCFILE.
    fields = QgsFields()
    fields.append(_aton_field("S57CLASS"))
    for fname in _ATON_BUOY_FIELDS:
        fields.append(_aton_field(fname))
    for fname in _ATON_TOPMARK_FIELDS:
        fields.append(_aton_field(_ATON_TT_RENAME[fname]))
    for fname in _ATON_LIGHT_FIELDS:
        fields.append(_aton_field(_ATON_LT_RENAME[fname]))
    for fname in _ATON_IDENT_FIELDS:
        fields.append(_aton_field(fname))
    fields.append(_aton_field("SRCFILE"))

    mem = QgsVectorLayer("Point?crs=EPSG:4326", "AtoN", "memory")
    mem_dp = mem.dataProvider()
    mem_dp.addAttributes(fields.toList())
    mem.updateFields()

    # ── Pass 1: collect buoy bodies, indexed by position so topmarks can pair.
    #    Gather EVERY case-variant buoy layer in each cell (a single cell can
    #    hold both standard `BOYLAT` and Inland `boylat`); opening by name alone
    #    would return just one and strand the other's topmarks.
    _buoy_fam_upper = _ATON_BODY_CLASSES_UPPER
    buoy_rows = []          # list of dicts: {geom, attrs{...}, pos}
    buoys_by_pos = {}       # pos_key -> index into buoy_rows (last wins)
    _crs_reported = False
    for fpath in usage_files:
        fname = os.path.basename(fpath)
        for nm, src in _iter_ogr_layers_ci(fpath, _buoy_fam_upper):
            cls = nm.upper()   # canonical class tag (BOYLAT/boylat -> BOYLAT)
            if not _crs_reported:
                # ── CRS DIAGNOSTIC ─ the AtoN merge ASSUMES buoy geometries
                #    are EPSG:4326 (it copies them raw into a 4326 memory
                #    layer). If GDAL reports anything else here, that
                #    assumption is the bug — the geometries would be
                #    mislabelled before reprojection to the target CRS.
                _scrs = src.crs()
                log(f"    [crs] AtoN source {nm}: "
                    f"{_scrs.authid() or 'EMPTY'} "
                    f"(valid={_scrs.isValid()}, "
                    f"geographic={_scrs.isGeographic()}) — "
                    f"merge assumes EPSG:4326", Qgis.Info)
                _crs_reported = True
            for f in src.getFeatures():
                g = f.geometry()
                if g is None or g.isEmpty():
                    continue
                if QgsWkbTypes.geometryType(g.wkbType()) != \
                        QgsWkbTypes.PointGeometry:
                    continue
                attrs = {"S57CLASS": cls, "SRCFILE": fname}
                for bf in _ATON_BUOY_FIELDS:
                    attrs[bf] = _aton_attr(f, bf)
                for idf in _ATON_IDENT_FIELDS:
                    attrs[idf] = _aton_attr(f, idf)
                row = {"geom": QgsGeometry(g), "attrs": attrs,
                       "pos": _aton_pos_key(g)}
                buoy_rows.append(row)
                if row["pos"] is not None:
                    buoys_by_pos[row["pos"]] = len(buoy_rows) - 1

    # ── Pass 2: collect topmarks; pair each onto a co-located buoy if present,
    #    otherwise keep it as a standalone topmark row.
    _tm_fam_upper = {_ATON_TOPMARK_CLASS.upper()}
    topmark_rows = []
    for fpath in usage_files:
        fname = os.path.basename(fpath)
        for nm, src in _iter_ogr_layers_ci(fpath, _tm_fam_upper):
            for f in src.getFeatures():
                g = f.geometry()
                if g is None or g.isEmpty():
                    continue
                if QgsWkbTypes.geometryType(g.wkbType()) != \
                        QgsWkbTypes.PointGeometry:
                    continue
                tt = {_ATON_TT_RENAME[tf]: _aton_attr(f, tf)
                      for tf in _ATON_TOPMARK_FIELDS}
                ident = {idf: _aton_attr(f, idf)
                         for idf in _ATON_IDENT_FIELDS}
                topmark_rows.append(
                    {"geom": QgsGeometry(g), "pos": _aton_pos_key(g), "tt": tt,
                     "standalone_attrs": {"S57CLASS": _ATON_TOPMARK_CLASS,
                                          "SRCFILE": fname, **tt, **ident}})

    if not buoy_rows and not topmark_rows:
        return None

    # ── Fold topmarks (by position) + pre-matched buoy-lights onto bodies via
    #    the shared reader-agnostic composer (also used by the Chart Viewer).
    rows, stats = compose_aton_rows(buoy_rows, topmark_rows, light_fold)
    log(f"[Merged] Aids to Navigation — {len(buoy_rows)} buoy(s), "
        f"{stats['paired']} with topmark, {stats['lit']} lit, "
        f"{stats['standalone']} lone topmark(s)")

    # ── Build features.
    feats = []
    for row in rows:
        feat = QgsFeature(mem.fields())
        feat.setGeometry(row["geom"])
        for k, v in row["attrs"].items():
            try:
                # Coerce to the column's strict type; missing → NULL (not '').
                feat.setAttribute(k, _aton_coerce(k, v))
            except Exception:
                pass
        feats.append(feat)
    mem_dp.addFeatures(feats)
    mem.updateExtents()

    # ── Write to the shared GeoPackage, reprojected to dest_crs.
    table = f"u{ub}_AtoN_point"
    # Priority 80 within the usage → above most areas/lines, with topmarks
    # drawn from the same row on top of the buoy body.
    lyr_name = f"U{ub} | {base_priority:02d} Aids to Navigation (point)"
    transform_context = QgsProject.instance().transformContext()
    options = QgsVectorFileWriter.SaveVectorOptions()
    options.driverName = "GPKG"
    options.layerName = table
    options.fileEncoding = "UTF-8"
    src_crs = mem.crs()
    if src_crs.isValid() and src_crs != dest_crs:
        options.ct = QgsCoordinateTransform(src_crs, dest_crs,
                                            transform_context)
    if os.path.exists(gpkg_path):
        options.actionOnExistingFile = \
            QgsVectorFileWriter.CreateOrOverwriteLayer
    else:
        options.actionOnExistingFile = \
            QgsVectorFileWriter.CreateOrOverwriteFile
    with _gdal_quiet():
        result = QgsVectorFileWriter.writeAsVectorFormatV3(
            mem, gpkg_path, transform_context, options
        )
    if result[0] != QgsVectorFileWriter.NoError:
        log(f"  ✗ AtoN write failed: {result[1]}", Qgis.Warning)
        return None

    actual_table = _find_gpkg_table(gpkg_path, table)
    if actual_table is None:
        log("  ✗ AtoN table not found after write.", Qgis.Warning)
        return None

    uri = f"{gpkg_path}|layername={actual_table}"
    out_layer = QgsVectorLayer(uri, lyr_name, "ogr")
    if not out_layer.isValid():
        log("  ✗ AtoN layer invalid after write.", Qgis.Warning)
        return None

    # Apply styles/AtoN_point.qml (buoy body + offset topmark symbol layers).
    cfg = get_config()
    sp = cfg.style_path_for("AtoN", "point")
    if sp:
        try:
            out_layer.loadNamedStyle(sp)
            out_layer.triggerRepaint()
            log(f"  ✓ Applied style {os.path.basename(sp)}")
        except Exception as e:
            log(f"  ⚠ AtoN style failed to load: {e}", Qgis.Warning)
    else:
        log("  ⚠ No styles/AtoN_point.qml found — layer left unstyled.",
            Qgis.Warning)

    log(f"  ✓ Aids to Navigation ready ({out_layer.featureCount()} feature(s)).")
    _ocrs = out_layer.crs()
    _oe = out_layer.extent()
    log(f"    [crs] AtoN output: {_ocrs.authid() or 'EMPTY'} "
        f"(valid={_ocrs.isValid()}) "
        f"extent=[{_oe.xMinimum():.2f},{_oe.yMinimum():.2f} .. "
        f"{_oe.xMaximum():.2f},{_oe.yMaximum():.2f}]", Qgis.Info)
    return out_layer


def _build_coverage_layer(best_files, dest_crs, gpkg_path, log,
                          usage=0, subprogress_fn=None):
    """Build an ENC Coverage hybrid layer for one usage band: one polygon per
    chart.

    For each winning chart file:
      • read M_COVR area features, keep CATCOV=1 (coverage available),
      • dissolve them into a single coverage polygon for that chart,
      • attach DSID metadata (EDTN UPDN ISDT DSNM PRSP PRED AGEN CSCL COMT;
        blank where the file doesn't provide it), plus CATCOV and a SRCFILE
        column naming the source .000.

    The result is reprojected to dest_crs, written to the shared GeoPackage,
    styled if a QML exists, and returned as a file-backed QgsVectorLayer.
    Returns None if no coverage geometry was found.
    """
    if not best_files:
        log("  ENC Coverage: no winning files to build from.", Qgis.Warning)
        return None

    log(f"[Hybrid] ENC Coverage — {len(best_files)} chart(s)")

    # Build the field schema: CATCOV + SRCFILE + the DSID fields.
    fields = QgsFields()
    fields.append(QgsField("SRCFILE", QVariant.String))
    fields.append(QgsField("CATCOV", QVariant.Int))
    for fname in _COVERAGE_DSID_FIELDS:
        # CSCL/EDTN/UPDN are conceptually numeric but S-57 may carry oddities;
        # store as string to never lose content, except CSCL as int-ish string.
        fields.append(QgsField(fname, QVariant.String))
    # Resolved S-62 producer name for AGEN (added right after the raw code).
    fields.append(QgsField("AGENCY", QVariant.String))

    # In-memory polygon layer in the S-57 native CRS (geographic, EPSG:4326).
    mem = QgsVectorLayer("MultiPolygon?crs=EPSG:4326", "ENC Coverage", "memory")
    mem_dp = mem.dataProvider()
    mem_dp.addAttributes(fields.toList())
    mem.updateFields()

    n_charts = len(best_files)
    built = 0
    feats = []
    for idx, (dsnm_key, fpath) in enumerate(sorted(best_files.items())):
        if subprogress_fn:
            subprogress_fn(int(idx / max(1, n_charts) * 90),
                           "Building coverage…")
        src = _open_ogr_layer(fpath, "M_COVR")
        if src is None or not src.isValid():
            continue

        # Collect CATCOV=1 area geometries from this chart.
        geoms = []
        catcov_seen = None
        for f in src.getFeatures():
            g = f.geometry()
            if g is None or g.isEmpty():
                continue
            if QgsWkbTypes.geometryType(g.wkbType()) != QgsWkbTypes.PolygonGeometry:
                continue
            catcov = None
            try:
                ci = f.fields().indexFromName("CATCOV")
                if ci >= 0:
                    catcov = f.attribute(ci)
            except Exception:
                catcov = None
            # default-filter to CATCOV=1 (coverage available)
            if catcov is not None and int(catcov) != 1:
                continue
            catcov_seen = 1 if catcov is None else int(catcov)
            geoms.append(g)

        if not geoms:
            continue

        # Dissolve into one polygon per chart.
        merged = geoms[0]
        for g in geoms[1:]:
            try:
                merged = merged.combine(g)
            except Exception:
                pass
        if merged is None or merged.isEmpty():
            continue
        # Normalise to multipolygon for a uniform schema.
        try:
            merged.convertToMultiType()
        except Exception:
            pass

        # DSID metadata (thorough read).
        meta = read_dsid_full(fpath)

        feat = QgsFeature(mem.fields())
        feat.setGeometry(merged)
        feat.setAttribute("SRCFILE", os.path.basename(fpath))
        feat.setAttribute("CATCOV", catcov_seen if catcov_seen is not None else 1)
        for fname in _COVERAGE_DSID_FIELDS:
            feat.setAttribute(fname, meta.get(fname, "") or "")
        # Resolve the S-62 producer code (AGEN) to a readable name. Graceful:
        # blank when AGEN is missing or not in the S-62 list.
        agency = get_config().agency_name(meta.get("AGEN"))
        feat.setAttribute("AGENCY", agency or "")
        feats.append(feat)
        built += 1

    if not feats:
        log("  ENC Coverage: no CATCOV=1 coverage polygons found.",
            Qgis.Warning)
        return None

    mem_dp.addFeatures(feats)
    mem.updateExtents()
    log(f"  ENC Coverage: {built} chart coverage polygon(s) assembled.")

    # Write the memory layer to the shared GeoPackage, reprojected to dest_crs.
    # Table + layer names are usage-qualified so each usage's coverage is its
    # own layer (and its own GPKG table).
    cov_table = f"u{usage}_ENC_Coverage" if usage else "ENC_Coverage"
    cov_name = (f"U{usage} | 90 ENC Coverage (area)" if usage
                else "90 ENC Coverage (area)")
    transform_context = QgsProject.instance().transformContext()
    options = QgsVectorFileWriter.SaveVectorOptions()
    options.driverName = "GPKG"
    options.layerName = cov_table
    options.fileEncoding = "UTF-8"
    src_crs = mem.crs()
    if src_crs.isValid() and src_crs != dest_crs:
        options.ct = QgsCoordinateTransform(src_crs, dest_crs, transform_context)
    if os.path.exists(gpkg_path):
        options.actionOnExistingFile = \
            QgsVectorFileWriter.CreateOrOverwriteLayer
    else:
        options.actionOnExistingFile = \
            QgsVectorFileWriter.CreateOrOverwriteFile
    with _gdal_quiet():
        result = QgsVectorFileWriter.writeAsVectorFormatV3(
            mem, gpkg_path, transform_context, options
        )
    if result[0] != QgsVectorFileWriter.NoError:
        log(f"  ✗ ENC Coverage write failed: {result[1]}", Qgis.Warning)
        return None

    actual_table = _find_gpkg_table(gpkg_path, cov_table)
    if actual_table is None:
        log("  ✗ ENC Coverage table not found after write.", Qgis.Warning)
        return None

    uri = f"{gpkg_path}|layername={actual_table}"
    out_layer = QgsVectorLayer(uri, cov_name, "ogr")
    if not out_layer.isValid():
        log("  ✗ ENC Coverage layer invalid after write.", Qgis.Warning)
        return None

    # Apply a style if the user shipped styles/ENC_COVERAGE_area.qml or
    # styles/M_COVR_area.qml (geometry-suffixed, per the new convention).
    cfg = get_config()
    for code in ("ENC_COVERAGE", "M_COVR"):
        sp = cfg.style_path_for(code, "area")
        if sp:
            try:
                out_layer.loadNamedStyle(sp)
                out_layer.triggerRepaint()
                log(f"  ✓ Applied coverage style {os.path.basename(sp)}")
            except Exception:
                pass
            break

    log(f"  ✓ ENC Coverage ready ({out_layer.featureCount()} feature(s)).")
    return out_layer
