# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
enc_s102.py
-----------
IHO S-102 (Bathymetric Surface, Edition 3.0.0) support for the Chart Viewer.

S-102 is an S-100 HDF5 *gridded coverage* (not ISO 8211), so — unlike the ENC
vector path — it is read with GDAL, which ships a dedicated, open-source S-102
driver (band 1 = depth, band 2 = uncertainty; georectified grid, S-100 Part 8).

This module:
  • finds S-102 .h5 datasets in the chart source folder,
  • reports their WGS84 coverage extent and tests consistency with the loaded
    S-57/S-101 charts (extent overlap),
  • loads the depth grid as a QgsRasterLayer portrayed with the SAME depth-band
    colours as DEPARE (shallow/safety/deep + the Day/Dusk/Night palette),
  • (helpers for) a safety-depth contour and thinned soundings derived from the
    grid — see enc_s102_derive.

Depth convention: GDAL returns S-102 band 1 as depth (positive down; negative =
drying height above datum), matching DRVAL1 — so the DEPARE band thresholds map
straight across.
"""

import os
import tempfile
import zipfile
import shutil

try:
    from qgis.core import (
        Qgis, QgsMessageLog, QgsProject, QgsRectangle, QgsRasterLayer,
        QgsVectorLayer, QgsField, QgsFeature, QgsGeometry, QgsPoint,
        QgsCoordinateReferenceSystem, QgsCoordinateTransform,
        QgsSingleBandPseudoColorRenderer, QgsRasterShader, QgsColorRampShader)
    from qgis.PyQt.QtGui import QColor
    from qgis.PyQt.QtCore import QVariant
except (ImportError, ValueError):                       # headless import guard
    QColor = None


def _log(msg, level=None):
    try:
        QgsMessageLog.logMessage(msg, "ENCore GIS",
                                 level if level is not None else Qgis.Info)
    except Exception:                                   # noqa: BLE001
        pass


def _gdal():
    from osgeo import gdal
    gdal.UseExceptions()
    return gdal


# ── Detection ───────────────────────────────────────────────────────────────
def is_s102(path):
    """True if GDAL identifies ``path`` as an S-102 dataset."""
    try:
        gdal = _gdal()
        drv = gdal.IdentifyDriver(path)
        if drv is not None and drv.ShortName == "S102":
            return True
    except Exception:                                   # noqa: BLE001
        pass
    # Fallback: S-102 cell-name convention (starts '102') + HDF5 extension.
    b = os.path.basename(path).lower()
    return b.startswith("102") and b.endswith((".h5", ".hdf5"))


def find_s102(folder):
    """All S-102 datasets under ``folder`` (recursively). Empty if none/!exists."""
    out = []
    if not folder or not os.path.isdir(folder):
        return out
    for dirpath, _dirs, files in os.walk(folder):
        for fn in files:
            if fn.lower().endswith((".h5", ".hdf5")):
                p = os.path.join(dirpath, fn)
                if is_s102(p):
                    out.append(p)
    return sorted(out)


_ZIP_CACHE = {}                                         # zip path -> [extracted .h5]


def _extract_s102_from_zip(zip_path):
    """Extract any HDF5 members of a .zip exchange set to a temp dir and return
    the S-102 ones (HDF5/GDAL can't read S-102 straight from /vsizip). Cached."""
    if zip_path in _ZIP_CACHE:
        return _ZIP_CACHE[zip_path]
    out = []
    try:
        with zipfile.ZipFile(zip_path) as z:
            members = [m for m in z.namelist()
                       if m.lower().endswith((".h5", ".hdf5"))]
            if members:
                dest = tempfile.mkdtemp(prefix="enc_s102_zip_")
                for m in members:
                    tgt = os.path.join(dest, os.path.basename(m))
                    with z.open(m) as src, open(tgt, "wb") as dst:
                        shutil.copyfileobj(src, dst)
                    if is_s102(tgt):
                        out.append(tgt)
    except Exception as e:                              # noqa: BLE001
        _log(f"S-102: reading zip failed ({os.path.basename(zip_path)}): {e}",
             Qgis.Warning)
    out = sorted(out)
    _ZIP_CACHE[zip_path] = out
    return out


def datasets_from_source(source):
    """S-102 datasets for a chart SOURCE — a folder (walked) or a .zip exchange
    set (HDF5 members extracted to temp). The .h5 are not in the cell staging
    dir, so we read them from where the charts were loaded."""
    if not source:
        return []
    if os.path.isdir(source):
        return find_s102(source)
    if zipfile.is_zipfile(source):
        return _extract_s102_from_zip(source)
    return []


# ── Coverage extent + consistency ────────────────────────────────────────────
def _open(path):
    """Open an S-102 dataset, descending into the first subdataset if the file
    exposes its coverage that way (GDAL ≥3.12 multi-group files)."""
    gdal = _gdal()
    ds = gdal.Open(path)
    if ds is not None and ds.RasterCount == 0:
        subs = ds.GetSubDatasets()
        if subs:
            ds = gdal.Open(subs[0][0])
    return ds


_DEPTH_CACHE = {}                                       # source path -> temp .tif


def _normalized_depth_tif(path):
    """Write (and cache) a single-band GeoTIFF of POSITIVE-DOWN depth from the
    S-102 grid, so surface/contour/soundings share one consistent convention.

    S-102 producers differ in vertical sign (NOAA: deep = large +; the Dutch
    CSMART set: deep = large −), and GDAL's DEPTH_OR_ELEVATION option is not
    consistent between them. So we detect the convention from the data — if the
    deepest extreme is negative (|min| > |max|), the grid is elevation-style and
    we negate it — and write a clean depth raster (deep = large positive,
    drying = negative, matching DRVAL1)."""
    if path in _DEPTH_CACHE and os.path.exists(_DEPTH_CACHE[path]):
        return _DEPTH_CACHE[path]
    try:
        gdal = _gdal()
        ds = _open(path)
        if ds is None:
            return None
        band = ds.GetRasterBand(1)
        arr = band.ReadAsArray().astype("float32")
        nd = band.GetNoDataValue()
        nd = 1.0e6 if nd is None else float(nd)
        mask = (arr == nd)
        valid = arr[~mask]
        if valid.size and abs(float(valid.min())) > abs(float(valid.max())):
            arr = -arr                                  # elevation → depth
        arr[mask] = nd                                  # restore nodata sentinel
        fd, tif = tempfile.mkstemp(suffix="_s102depth.tif")
        os.close(fd)
        out = gdal.GetDriverByName("GTiff").Create(
            tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32)
        out.SetGeoTransform(ds.GetGeoTransform())
        out.SetProjection(ds.GetProjection())
        ob = out.GetRasterBand(1)
        ob.SetNoDataValue(nd)
        ob.WriteArray(arr)
        ob.SetDescription("depth")
        out.FlushCache()
        out = None
        _DEPTH_CACHE[path] = tif
        return tif
    except Exception as e:                              # noqa: BLE001
        _log(f"S-102 depth normalisation failed for "
             f"{os.path.basename(path)}: {e}", Qgis.Warning)
        return None


_MOSAIC_CACHE = {}                                      # tuple(paths) -> .vrt


def _depth_source(paths):
    """A single GDAL-openable depth source for one OR several S-102 datasets:
    each is normalised to a positive-down depth GeoTIFF, and multiple tiles are
    mosaicked into a VRT — so the surface / contour / soundings cover the WHOLE
    loaded coverage, not just one tile. (VRT mosaic assumes a common CRS, true
    for adjacent tiles in one delivery; falls back to the first tile otherwise.)"""
    if isinstance(paths, str):
        paths = [paths]
    tifs = [t for t in (_normalized_depth_tif(p) for p in (paths or [])) if t]
    if not tifs:
        return None
    if len(tifs) == 1:
        return tifs[0]
    key = tuple(sorted(paths))
    if key in _MOSAIC_CACHE and os.path.exists(_MOSAIC_CACHE[key]):
        return _MOSAIC_CACHE[key]
    try:
        gdal = _gdal()
        fd, vrt = tempfile.mkstemp(suffix="_s102_mosaic.vrt")
        os.close(fd)
        ds = gdal.BuildVRT(vrt, tifs)
        if ds is None:
            raise RuntimeError("BuildVRT returned None")
        ds = None
        _MOSAIC_CACHE[key] = vrt
        return vrt
    except Exception as e:                              # noqa: BLE001
        _log(f"S-102: mosaic failed ({e}); using one tile only", Qgis.Warning)
        return tifs[0]


def extent_wgs84(path):
    """S-102 coverage bounding box as a WGS84 QgsRectangle, or None."""
    try:
        ds = _open(path)
        if ds is None:
            return None
        gt = ds.GetGeoTransform()
        nx, ny = ds.RasterXSize, ds.RasterYSize
        xs = (gt[0], gt[0] + nx * gt[1])
        ys = (gt[3], gt[3] + ny * gt[5])
        src = QgsCoordinateReferenceSystem.fromWkt(ds.GetProjection())
        wgs = QgsCoordinateReferenceSystem("EPSG:4326")
        rect = QgsRectangle()
        rect.setMinimal()
        if src.isValid() and src != wgs:
            xf = QgsCoordinateTransform(src, wgs, QgsProject.instance())
            for x in xs:
                for y in ys:
                    p = xf.transform(x, y)
                    rect.include(p)
        else:
            rect.set(min(xs), min(ys), max(xs), max(ys))
        return rect
    except Exception as e:                              # noqa: BLE001
        _log(f"S-102: extent read failed for {os.path.basename(path)}: {e}",
             Qgis.Warning)
        return None


def metadata(path):
    """Dataset/producer metadata GDAL exposes for an S-102 file.

    Future S-102 editions carry richer discovery metadata (issue date, vertical
    datum, producing agency, data assessment…). This surfaces whatever the
    driver reports today so the viewer/info panel can consume it as the standard
    grows — currently a flat {key: value} dict (may be empty on Ed 3.0.0 files)."""
    try:
        ds = _open(path)
        if ds is None:
            return {}
        md = dict(ds.GetMetadata() or {})
        try:
            band = ds.GetRasterBand(1)
            md.setdefault("band1", band.GetDescription())
            md.setdefault("nodata", band.GetNoDataValue())
        except Exception:                               # noqa: BLE001
            pass
        return md
    except Exception:                                   # noqa: BLE001
        return {}


def coverage_consistent(s102_paths, charts_extent_wgs84, min_overlap=0.05):
    """Which S-102 datasets are consistent with the loaded charts: their WGS84
    coverage must overlap the charts' extent (intersection ≥ min_overlap of the
    S-102 footprint). Returns the list of consistent paths."""
    if charts_extent_wgs84 is None or charts_extent_wgs84.isEmpty():
        return []
    ok = []
    for p in s102_paths:
        ext = extent_wgs84(p)
        if ext is None or ext.isEmpty():
            continue
        inter = ext.intersect(charts_extent_wgs84)
        if inter.isEmpty():
            continue
        frac = (inter.area() / ext.area()) if ext.area() > 0 else 0.0
        if frac >= min_overlap:
            ok.append(p)
            _log(f"S-102 {os.path.basename(p)}: consistent "
                 f"({frac * 100:.0f}% of its footprint overlaps the charts)")
        else:
            _log(f"S-102 {os.path.basename(p)}: only {frac * 100:.0f}% overlap "
                 "— treated as inconsistent", Qgis.Warning)
    return ok


# ── Surface portrayal (DEPARE colours) ───────────────────────────────────────
def load_surface(paths, target_crs, name="S-102 Bathymetric surface"):
    """The S-102 depth surface as a QgsRasterLayer (band 1 = positive-down depth),
    mosaicked across all given datasets. Reprojection to the viewer CRS is
    on-the-fly. Returns the layer (caller styles it) or None."""
    uri = _depth_source(paths)
    if uri is None:
        return None
    rlayer = QgsRasterLayer(uri, name, "gdal")
    if not rlayer.isValid():
        _log(f"S-102: QgsRasterLayer invalid for {uri}", Qgis.Warning)
        return None
    return rlayer


def load_surface_clips(paths, boxes_wgs84, name="S-102 surface (tide)"):
    """Footprint-clipped copies of the depth mosaic — one QgsRasterLayer per
    WGS84 (w, s, e, n) box — so the tide-shifted depth shading can be confined
    to the S-104 footprint: the viewer stacks these ABOVE the base surface
    (which keeps the charted shading) and shades only them at the shifted
    thresholds. Boxes are clamped to the raster; empty overlaps are skipped."""
    try:
        from osgeo import gdal
    except Exception:                                   # noqa: BLE001
        return []
    src = _depth_source(paths)
    ds = gdal.Open(src) if src else None
    if ds is None or not boxes_wgs84:
        return []
    gt = ds.GetGeoTransform()
    xs = (gt[0], gt[0] + ds.RasterXSize * gt[1])
    ys = (gt[3], gt[3] + ds.RasterYSize * gt[5])
    ext = QgsRectangle(min(xs), min(ys), max(xs), max(ys))
    src_crs = QgsCoordinateReferenceSystem.fromWkt(ds.GetProjection())
    wgs = QgsCoordinateReferenceSystem("EPSG:4326")
    xf = (QgsCoordinateTransform(wgs, src_crs, QgsProject.instance())
          if src_crs.isValid() and src_crs != wgs else None)
    out = []
    for i, (w, s, e, n) in enumerate(boxes_wgs84):
        r = QgsRectangle(w, s, e, n)
        try:
            if xf is not None:
                r = xf.transformBoundingBox(r)
        except Exception:                               # noqa: BLE001
            continue
        r = r.intersect(ext)                            # projWin must stay inside
        if r.isEmpty():
            continue
        fd, vrt = tempfile.mkstemp(suffix=f"_s102_tide{i}.vrt")
        os.close(fd)
        try:
            clip = gdal.Translate(vrt, ds, projWin=[
                r.xMinimum(), r.yMaximum(), r.xMaximum(), r.yMinimum()])
            if clip is None:
                continue
            clip = None                                 # flush the VRT
            rl = QgsRasterLayer(vrt, name, "gdal")
            if rl.isValid():
                out.append(rl)
        except Exception as exc:                        # noqa: BLE001
            _log(f"S-102 tide clip failed: {exc}", Qgis.Warning)
    return out


def clip_lines_to_boxes(layer, boxes_wgs84, inside):
    """Clip a derived contour layer IN PLACE to the union of the WGS84
    (w, s, e, n) boxes: ``inside=True`` keeps only what falls within them,
    ``inside=False`` only what falls outside — the two halves of a footprint-
    gated contour (tide-shifted cut inside, charted cut outside). Returns the
    surviving feature count."""
    if layer is None or not boxes_wgs84:
        return 0 if layer is None else layer.featureCount()
    wgs = QgsCoordinateReferenceSystem("EPSG:4326")
    xf = (None if layer.crs() == wgs else
          QgsCoordinateTransform(wgs, layer.crs(), QgsProject.instance()))
    union = None
    for (w, s, e, n) in boxes_wgs84:
        r = QgsRectangle(w, s, e, n)
        if xf is not None:
            r = xf.transformBoundingBox(r)
        g = QgsGeometry.fromRect(r)
        union = g if union is None else union.combine(g)
    from qgis.core import QgsWkbTypes
    dp = layer.dataProvider()
    drop, change = [], {}
    for ft in layer.getFeatures():
        g = ft.geometry()
        ng = g.intersection(union) if inside else g.difference(union)
        if ng is not None and not ng.isEmpty():
            # An intersection can return a collection (touch points etc.);
            # keep only the line parts — the layer is MultiLineString.
            ng = ng.convertToType(QgsWkbTypes.LineGeometry, True)
        if ng is None or ng.isEmpty():
            drop.append(ft.id())
        else:
            change[ft.id()] = ng
    if change:
        dp.changeGeometryValues(change)
    if drop:
        dp.deleteFeatures(drop)
    layer.updateExtents()
    return layer.featureCount()


def merge_contours(dst, src):
    """Append ``src``'s features into ``dst`` (same schema — both come from
    _contour_layer). Either may be None; returns the combined layer, which
    keeps ``dst``'s styling."""
    if src is None:
        return dst
    if dst is None:
        return src
    dst.dataProvider().addFeatures(list(src.getFeatures()))
    dst.updateExtents()
    return dst


def depth_band_shader_items(shallow, safety, deep, bands):
    """Discrete colour-ramp items (upper-edge value → QColor) matching the
    DEPARE bands: <0 drying, 0–shallow, shallow–safety, safety–deep, >deep."""
    def col(label):
        rgb = bands.get(label, (200, 200, 200))
        return QColor(*rgb)
    Item = QgsColorRampShader.ColorRampItem
    big = 1.0e6                                          # also masks nodata edge
    return [
        Item(0.0,           col("above LAT"),     "drying"),
        Item(float(shallow), col("above Shallow"), f"0–{shallow} m"),
        Item(float(safety),  col("above Safety"),  f"{shallow}–{safety} m"),
        Item(float(deep),    col("below Safety"),  f"{safety}–{deep} m"),
        Item(big,            col("Deep"),          f">{deep} m"),
    ]


def apply_depth_shader(rlayer, shallow, safety, deep, bands):
    """Render the S-102 surface (band 1) with the DEPARE discrete depth bands.
    Rebuilds the renderer, so it can be called live on Chart Control / palette
    changes."""
    if rlayer is None or not rlayer.isValid():
        return
    ramp = QgsColorRampShader()
    ramp.setColorRampType(QgsColorRampShader.Discrete)
    ramp.setColorRampItemList(
        depth_band_shader_items(shallow, safety, deep, bands))
    shader = QgsRasterShader()
    shader.setRasterShaderFunction(ramp)
    renderer = QgsSingleBandPseudoColorRenderer(
        rlayer.dataProvider(), 1, shader)
    rlayer.setRenderer(renderer)
    rlayer.triggerRepaint()


# ── Derived products: safety contour + soundings ─────────────────────────────
def _src_xform(ds, target_crs):
    """(src_crs, transform-to-target) for an open S-102 dataset."""
    src = QgsCoordinateReferenceSystem.fromWkt(ds.GetProjection())
    dst = QgsCoordinateReferenceSystem(target_crs)
    if not src.isValid() or not dst.isValid() or src == dst:
        return src, None
    return src, QgsCoordinateTransform(src, dst, QgsProject.instance())


# ── Shoal-safe contour smoothing (IHO S-4: declutter only, SHOAL-BIASED) ─────
# A smoothed depth contour must NEVER indicate deeper water than reality. We
# smooth the depth grid with a nodata-aware box average, then keep the per-cell
# SHOALEST of (original, smoothed). Depth is positive-down, so shoalest = minimum
# ⇒ smoothed ≤ original everywhere ⇒ the "deep-enough" region can only shrink ⇒
# provably shoal-safe. Net effect: the deep side of a contour is smoothed while
# shoal intrusions are preserved/expanded (shoal bias) — the CARIS/HPD look.
_SMOOTH_RADIUS = 1      # box half-width in grid cells
_SMOOTH_ITERS = 4       # averaging passes (Laplacian-style declutter)


def _box_sum(a, r):
    """Sum of ``a`` over each (2r+1)×(2r+1) window (zero-padded edges), O(n) via a
    summed-area table — independent of the window radius."""
    import numpy as np
    cs = np.cumsum(np.cumsum(a, axis=0), axis=1)
    cs = np.pad(cs, ((1, 0), (1, 0)))                   # leading zero row + col
    h, w = a.shape
    r0 = np.arange(h)
    c0 = np.arange(w)
    top = np.clip(r0 - r, 0, h)
    bot = np.clip(r0 + r + 1, 0, h)
    lft = np.clip(c0 - r, 0, w)
    rgt = np.clip(c0 + r + 1, 0, w)
    return (cs[np.ix_(bot, rgt)] - cs[np.ix_(top, rgt)]
            - cs[np.ix_(bot, lft)] + cs[np.ix_(top, lft)])


def _shoal_smooth_band(band, radius=_SMOOTH_RADIUS, iters=_SMOOTH_ITERS):
    """An in-memory single-band GDAL dataset holding the SHOAL-SAFE smoothed depth
    grid of ``band`` (same geotransform/projection/nodata), or None on error."""
    import numpy as np
    gdal = _gdal()
    try:
        arr = band.ReadAsArray().astype("float64")
    except Exception:                                   # noqa: BLE001
        return None
    if arr is None:
        return None
    nodata = band.GetNoDataValue()
    valid = np.isfinite(arr)
    if nodata is not None:
        valid &= (arr != nodata)
    if not valid.any():
        return None
    orig = arr.copy()
    wsum = _box_sum(valid.astype("float64"), radius)    # valid-cell count / window
    cur = np.where(valid, arr, 0.0)
    for _ in range(max(1, iters)):                      # nodata-aware box average
        num = _box_sum(np.where(valid, cur, 0.0), radius)
        cur = np.divide(num, wsum, out=cur.copy(), where=wsum > 0)
    # shoal-biased clamp: never deeper than the real seabed (positive-down ⇒ min)
    out = np.minimum(orig, cur)
    fill = float(nodata) if nodata is not None else 0.0
    out = np.where(valid, out, fill)
    ds = band.GetDataset()
    h, w = arr.shape
    mem = gdal.GetDriverByName("MEM").Create("", w, h, 1, gdal.GDT_Float32)
    mem.SetGeoTransform(ds.GetGeoTransform())
    mem.SetProjection(ds.GetProjection())
    ob = mem.GetRasterBand(1)
    if nodata is not None:
        ob.SetNoDataValue(float(nodata))
    ob.WriteArray(out.astype("float32"))
    return mem


def _contour_layer(paths, depth, target_crs, name, smooth=True, valdco=None):
    """Generate one contour line from the S-102 coverage at ``depth`` → an
    unstyled MultiLineString QgsVectorLayer (VALDCO = depth). Shared by the
    safety contour (DEPCNT-styled) and the extra contour (red). When ``smooth``,
    the depth grid is shoal-safely smoothed first (declutter, never deeper)."""
    try:
        from osgeo import gdal, ogr
    except Exception:                                   # noqa: BLE001
        return None
    src = _depth_source(paths)                          # positive-down depth
    ds = gdal.Open(src) if src else None
    if ds is None:
        return None
    band = ds.GetRasterBand(1)
    smoothed_ds = _shoal_smooth_band(band) if smooth else None
    if smoothed_ds is not None:                         # keep ref alive for contour
        band = smoothed_ds.GetRasterBand(1)
    mem = ogr.GetDriverByName("Memory").CreateDataSource("c")
    clayer = mem.CreateLayer("c", geom_type=ogr.wkbLineString25D)
    clayer.CreateField(ogr.FieldDefn("ID", ogr.OFTInteger))
    clayer.CreateField(ogr.FieldDefn("ELEV", ogr.OFTReal))
    opts = ["ID_FIELD=0", "ELEV_FIELD=1", f"FIXED_LEVELS={float(depth)}"]
    nodata = band.GetNoDataValue()
    if nodata is not None:
        opts.append(f"NODATA={nodata}")
    try:
        gdal.ContourGenerateEx(band, clayer, options=opts)
    except Exception as e:                              # noqa: BLE001
        _log(f"S-102 contour at {depth} m failed: {e}", Qgis.Warning)
        return None
    _src, xf = _src_xform(ds, target_crs)
    dstc = QgsCoordinateReferenceSystem(target_crs)
    layer = QgsVectorLayer(
        f"MultiLineString?crs={dstc.authid() or 'EPSG:4326'}", name, "memory")
    dp = layer.dataProvider()
    dp.addAttributes([QgsField("VALDCO", QVariant.String)])
    layer.updateFields()
    feats = []
    for of in clayer:
        geom = of.GetGeometryRef()
        if geom is None:
            continue
        g = QgsGeometry.fromWkt(geom.ExportToWkt())
        if g is None or g.isEmpty():
            continue
        if xf is not None:
            g.transform(xf)
        ft = QgsFeature(layer.fields())
        ft.setGeometry(g)
        # VALDCO carries the contour's IDENTITY (the emphasis value the DEPCNT QML
        # bolds on), which may differ from the GEOMETRY depth: under the S-104 tide
        # the safety line is CUT at (safety − level) charted metres but is still THE
        # safety contour, so it must keep VALDCO = safety to stay bold.
        ft["VALDCO"] = str(valdco if valdco is not None else depth)
        feats.append(ft)
    dp.addFeatures(feats)
    layer.updateExtents()
    return layer, len(feats)


def safety_contour(paths, safety_depth, target_crs,
                   name="S-102 safety contour", valdco=None):
    """The S-102 safety-depth contour, drawn like the DEPCNT safety contour. The
    line geometry is cut at ``safety_depth`` (the tide-shifted CHARTED depth when
    S-104 is live); pass ``valdco`` = the true safety value so the DEPCNT QML's
    'VALDCO >= @enc_safety_contour' rule still emphasises it as the safety line."""
    res = _contour_layer(paths, safety_depth, target_crs, name, valdco=valdco)
    if res is None:
        return None
    layer, n = res
    _apply_s57_qml(layer, "DEPCNT", "line")
    _log(f"S-102 safety contour: {n} part(s) at {safety_depth} m "
         f"(VALDCO {valdco if valdco is not None else safety_depth})")
    return layer


def extra_contour(paths, depth, target_crs, name="S-102 extra contour"):
    """The S-102 'extra contour' (Chart Control), drawn as a bold RED line."""
    res = _contour_layer(paths, depth, target_crs, name)
    if res is None:
        return None
    layer, n = res
    try:
        from qgis.core import (QgsLineSymbol, QgsSingleSymbolRenderer)
        sym = QgsLineSymbol.createSimple(
            {"line_color": "255,0,0,255", "line_width": "0.4",
             "line_width_unit": "MM"})
        layer.setRenderer(QgsSingleSymbolRenderer(sym))
    except Exception:                                   # noqa: BLE001
        pass
    _log(f"S-102 extra contour: {n} part(s) at {depth} m (red)")
    return layer


def soundings(paths, target_crs, target_count=6000,
              name="S-102 soundings"):
    """Sounding points sampled from the S-102 coverage (depth = Z), thinned with
    the same adaptive shoalest-first algorithm + SOUNDG portrayal as ENC
    soundings. Sub-samples the (mosaicked) grid to ~target_count candidates."""
    try:
        gdal = _gdal()
    except Exception:                                   # noqa: BLE001
        return None
    src = _depth_source(paths)                          # positive-down depth
    ds = gdal.Open(src) if src else None
    if ds is None:
        return None
    band = ds.GetRasterBand(1)
    nx, ny = ds.RasterXSize, ds.RasterYSize
    gt = ds.GetGeoTransform()
    nodata = band.GetNoDataValue()
    step = max(1, int((nx * ny / max(target_count, 1)) ** 0.5))
    arr = band.ReadAsArray()
    _src, xf = _src_xform(ds, target_crs)
    dstc = QgsCoordinateReferenceSystem(target_crs)
    layer = QgsVectorLayer(f"PointZ?crs={dstc.authid() or 'EPSG:4326'}",
                           name, "memory")
    dp = layer.dataProvider()
    feats = []
    for row in range(0, ny, step):
        for col in range(0, nx, step):
            z = float(arr[row, col])
            if nodata is not None and z == nodata:
                continue
            if z != z:                                  # NaN
                continue
            x = gt[0] + (col + 0.5) * gt[1] + (row + 0.5) * gt[2]
            y = gt[3] + (col + 0.5) * gt[4] + (row + 0.5) * gt[5]
            g = QgsGeometry(QgsPoint(x, y, z))
            if xf is not None:
                g.transform(xf)
            ft = QgsFeature()
            ft.setGeometry(g)
            feats.append(ft)
    dp.addFeatures(feats)
    layer.setCrs(dstc)
    layer.updateExtents()
    # Adaptive thinning + SOUNDG portrayal (reuse the ENC machinery).
    try:
        from .enc_loader import _assign_zoom_levels
        from .enc_chartview import _soundg_params
        from .enc_config import get_config
    except (ImportError, ValueError):
        from enc_loader import _assign_zoom_levels
        from enc_chartview import _soundg_params
        from enc_config import get_config
    layer = _assign_zoom_levels(layer, _soundg_params(get_config()), _log)
    layer.setName(name)
    layer.setCrs(dstc)
    _apply_s57_qml(layer, "SOUNDG", "point")
    _log(f"S-102 soundings: {len(feats)} sampled (step {step}px) → thinned")
    return layer


def _apply_s57_qml(layer, code, kind):
    """Apply the bundled ENC QML (DEPCNT/SOUNDG) so derived S-102 products are
    portrayed identically to their ENC counterparts (Chart Control-driven)."""
    try:
        from .enc_loader import _apply_portrayal
        from .enc_config import get_config
    except (ImportError, ValueError):
        from enc_loader import _apply_portrayal
        from enc_config import get_config
    try:
        _apply_portrayal(layer, get_config(), code, kind, _log)
    except Exception as e:                              # noqa: BLE001
        _log(f"S-102 derived portrayal failed for {code}: {e}", Qgis.Warning)
