# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
enc_s104.py
-----------
S-104 Water Level Information — read the water level for a time and feed the
Chart Viewer's tidal-adjustment engine.

S-104 (S-100 HDF5) has NO portrayal catalogue by design: this edition is for
"water level adjustment on ECDIS". So we don't draw it — we read a
representative water-level value for the selected time and ADD it to the wet
layers (depth bands / contours via a project variable, soundings via the
displayed value). Detection mirrors enc_s102 (folder or .zip).

Read path note: GDAL's classic S104 driver reads the metadata but FAILS on the
compound Group_NNN/values ("Wrong data type"), so the values are read with the
MULTIDIM API (RootGroup → WaterLevel → WaterLevel.01 → Group_NNN → values), which
returns a structured array {waterLevelHeight, waterLevelTrend, uncertainty}. The
S-104 grid is coarse and the tide is near-uniform over a chart area, so a single
representative scalar (mean of the valid heights) is the value the engine applies.
"""

import os
import shutil
import tempfile
import zipfile
from datetime import datetime, timedelta

try:
    from qgis.core import Qgis, QgsMessageLog
except (ImportError, ValueError):                       # headless import guard
    Qgis = None

_NODATA = -9999.0


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 (folder or .zip), mirroring enc_s102 / enc_s111 ────────────────
def is_s104(path):
    try:
        drv = _gdal().IdentifyDriver(path)
        if drv is not None and drv.ShortName == "S104":
            return True
    except Exception:                                   # noqa: BLE001
        pass
    b = os.path.basename(path).lower()
    return b.startswith("104") and b.endswith((".h5", ".hdf5"))


def find_s104(folder):
    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_s104(p):
                    out.append(p)
    return sorted(out)


_ZIP_CACHE = {}


def _extract_s104_from_zip(zip_path):
    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_s104_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_s104(tgt):
                        out.append(tgt)
    except Exception as e:                              # noqa: BLE001
        _log(f"S-104: 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-104 datasets for a chart source — a folder (walked, incl. any .zip found)
    or a .zip (HDF5 members extracted). A folder may itself hold S104-*.zip
    exchange sets, so those are unpacked too."""
    if not source:
        return []
    if os.path.isdir(source):
        out = list(find_s104(source))
        for dirpath, _dirs, files in os.walk(source):
            for fn in files:
                if fn.lower().endswith(".zip"):
                    out += _extract_s104_from_zip(os.path.join(dirpath, fn))
        return sorted(set(out))
    if zipfile.is_zipfile(source):
        return _extract_s104_from_zip(source)
    return []


# ── Time axis (same shape as S-111) ──────────────────────────────────────────
def _as_list(paths):
    return list(paths) if isinstance(paths, (list, tuple)) else [paths]


def _parse_dt(s):
    for fmt in ("%Y%m%dT%H%M%SZ", "%Y%m%dT%H%M%S", "%Y-%m-%dT%H:%M:%SZ"):
        try:
            return datetime.strptime(str(s).strip(), fmt)
        except (ValueError, TypeError):
            continue
    return None


def time_axis(paths):
    """UNION time axis over ``paths`` as a sorted list of datetimes (from
    dateTimeOfFirstRecord + timeRecordInterval + numberOfTimes)."""
    gdal = _gdal()
    times = set()
    for p in _as_list(paths):
        try:
            md = gdal.Open(p).GetMetadata()             # classic driver reads md
            first = _parse_dt(md.get("dateTimeOfFirstRecord"))
            interval = int(float(md.get("timeRecordInterval") or 0))
            n = int(float(md.get("numberOfTimes") or 0))
            if first and interval and n:
                for i in range(n):
                    times.add(first + timedelta(seconds=interval * i))
            elif first:
                times.add(first)
        except Exception as e:                          # noqa: BLE001
            _log(f"S-104 time axis failed ({os.path.basename(p)}): {e}",
                 Qgis.Warning)
    return sorted(times)


def _nearest_index(md, when, n):
    if when is None:
        return 0
    first = _parse_dt(md.get("dateTimeOfFirstRecord"))
    interval = int(float(md.get("timeRecordInterval") or 0))
    if first and interval:
        return max(0, min(n - 1, int(round((when - first).total_seconds()
                                           / interval))))
    return 0


def water_level_at(paths, when):
    """A single representative water-level height (m, vs vertical datum) for the
    time nearest ``when``, averaged over the (near-uniform, coarse) S-104 grid of
    the given dataset(s). Read via the multidim API (the classic S104 driver can't
    read the compound values). Returns a float, or 0.0 if nothing readable."""
    import numpy as np
    gdal = _gdal()
    vals = []
    for p in _as_list(paths):
        try:
            md = gdal.Open(p).GetMetadata()
            n = int(float(md.get("numberOfTimes") or 0)) or 1
            idx = _nearest_index(md, when, n)
            ds = gdal.OpenEx(p, gdal.OF_MULTIDIM_RASTER)
            inst = ds.GetRootGroup().OpenGroup("WaterLevel").OpenGroup(
                "WaterLevel.01")
            grp = inst.OpenGroup(f"Group_{idx + 1:03d}")
            arr = grp.OpenMDArray("values").ReadAsArray()
            h = np.asarray(arr["waterLevelHeight"], dtype="float64")
            valid = h[h > _NODATA + 1]
            if valid.size:
                vals.append(float(valid.mean()))
        except Exception as e:                          # noqa: BLE001
            _log(f"S-104 read failed ({os.path.basename(p)}): {e}", Qgis.Warning)
    return float(sum(vals) / len(vals)) if vals else 0.0


def _attr_f(group, name):
    """Read a multidim group attribute as float, or None."""
    try:
        a = group.GetAttribute(name)
        return float(a.ReadAsDouble()) if a is not None else None
    except Exception:                                   # noqa: BLE001
        return None


def footprints(paths):
    """WGS84 bounding boxes ``[(xmin, ymin, xmax, ymax)]`` of the S-104 grids —
    the footprint the water-level adjustment APPLIES to. Read from the
    instance group's bounding-box attributes (multidim API; the classic driver
    does not surface them). Used to draw the coverage outline and to GATE the
    tide-adjusted portrayal: outside these boxes the chart must stay un-tided."""
    gdal = _gdal()
    out = []
    for p in _as_list(paths):
        try:
            ds = gdal.OpenEx(p, gdal.OF_MULTIDIM_RASTER)
            inst = ds.GetRootGroup().OpenGroup("WaterLevel").OpenGroup(
                "WaterLevel.01")
            w = _attr_f(inst, "westBoundLongitude")
            e = _attr_f(inst, "eastBoundLongitude")
            s = _attr_f(inst, "southBoundLatitude")
            n = _attr_f(inst, "northBoundLatitude")
            if None not in (w, e, s, n) and w < e and s < n:
                out.append((w, s, e, n))
        except Exception as exc:                        # noqa: BLE001
            _log(f"S-104 footprint failed ({os.path.basename(p)}): {exc}",
                 Qgis.Warning)
    return out


def water_level_series(paths, lon, lat):
    """Water-level time series at the grid node nearest (lon, lat) — a list of
    (datetime, height_m vs datum), over the full record. Among ``paths`` it picks
    the dataset whose regular grid actually covers the point (falling back to the
    nearest grid origin), then reads that single node across every Group_NNN.
    Returns [] if nothing is readable. (~0.6 s for 325 steps on the coarse grid.)"""
    import numpy as np
    gdal = _gdal()
    best = None
    for p in _as_list(paths):
        try:
            md = gdal.Open(p).GetMetadata()
            first = _parse_dt(md.get("dateTimeOfFirstRecord"))
            interval = int(float(md.get("timeRecordInterval") or 0))
            n = int(float(md.get("numberOfTimes") or 0))
            ds = gdal.OpenEx(p, gdal.OF_MULTIDIM_RASTER)
            inst = (ds.GetRootGroup().OpenGroup("WaterLevel")
                    .OpenGroup("WaterLevel.01"))
            olat = _attr_f(inst, "gridOriginLatitude")
            olon = _attr_f(inst, "gridOriginLongitude")
            slat = _attr_f(inst, "gridSpacingLatitudinal")
            slon = _attr_f(inst, "gridSpacingLongitudinal")
            nlat = int(_attr_f(inst, "numPointsLatitudinal") or 0)
            nlon = int(_attr_f(inst, "numPointsLongitudinal") or 0)
            if (None in (olat, olon, slat, slon) or not (nlat and nlon)
                    or not (first and interval and n)):
                continue
            # nearest node (value array is [lat_row i, lon_col j]) + coverage test
            i = int(round((lat - olat) / slat)) if slat else 0
            j = int(round((lon - olon) / slon)) if slon else 0
            inside = (0 <= i < nlat and 0 <= j < nlon)
            i = max(0, min(nlat - 1, i))
            j = max(0, min(nlon - 1, j))
            dist = ((olat + i * slat) - lat) ** 2 + ((olon + j * slon) - lon) ** 2
            cand = {"inst": inst, "i": i, "j": j, "first": first,
                    "interval": interval, "n": n, "inside": inside, "dist": dist}
            if (best is None
                    or (cand["inside"] and not best["inside"])
                    or (cand["inside"] == best["inside"] and dist < best["dist"])):
                best = cand
        except Exception as e:                          # noqa: BLE001
            _log(f"S-104 series read failed ({os.path.basename(p)}): {e}",
                 Qgis.Warning)
    if best is None:
        return []
    inst, i, j = best["inst"], best["i"], best["j"]
    series = []
    for k in range(best["n"]):
        try:
            grp = inst.OpenGroup(f"Group_{k + 1:03d}")
            arr = grp.OpenMDArray("values").ReadAsArray()
            h = float(np.asarray(arr["waterLevelHeight"])[i, j])
            if h > _NODATA + 1:
                series.append(
                    (best["first"] + timedelta(seconds=best["interval"] * k), h))
        except Exception:                               # noqa: BLE001
            continue
    return series
