Reading HEC-DSS Time Series into pandas

import logging
import re

import numpy as np
import pandas as pd

log = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")

# The DSS missing-value sentinel. Reading it as a number rather than as a gap
# is how a hydrograph acquires a -3.4e38 spike that ruins every statistic.
DSS_MISSING = -3.4028234663852886e38


def read_dss_series(
    dss_path: str,
    pathname: str,
    stamp: str = "end",
    tz: str | None = None,
) -> pd.Series:
    """
    Read one regular-interval DSS series into a pandas Series.

    Parameters
    ----------
    dss_path : Path to the .dss file.
    pathname : Six-part pathname. Part D should be EMPTY to span all blocks.
    stamp    : "end" if the series is stamped at interval ends (the DSS
               convention for regular series) — the index is shifted back to
               interval starts so it aligns with observed data read from CSV.
    tz       : Optional timezone to localise to.
    """
    from hecdss import HecDss

    parts = pathname.strip("/").split("/")
    if len(parts) != 6:
        raise ValueError(f"expected six pathname parts, got {len(parts)}")
    if parts[3]:
        log.warning("Pathname carries a D part (%r) — this matches ONE storage "
                    "block only. Blanking it to read the whole series", parts[3])
        parts[3] = ""
        pathname = "/" + "/".join(parts) + "/"

    with HecDss(dss_path) as dss:
        record = dss.get(pathname)
        if record is None:
            catalogue = [p for p in dss.get_catalog()
                         if parts[1] in p and parts[2] in p]
            log.error("No record at %s. %d similar pathname(s) exist: %s",
                      pathname, len(catalogue), catalogue[:3])
            raise KeyError(f"no DSS record at {pathname}")

        values = np.asarray(record.values, dtype="float64")
        times = pd.to_datetime(record.times)
        units = getattr(record, "units", "")

    # --- Sentinel to NaN before anything else. Every later statistic is
    # meaningless if a -3.4e38 survives into the array. ---
    n_missing = int((values <= DSS_MISSING / 2).sum())
    values = np.where(values <= DSS_MISSING / 2, np.nan, values)
    if n_missing:
        log.warning("%d missing-value sentinel(s) converted to NaN", n_missing)

    series = pd.Series(values, index=times, name=parts[2])

    # --- Shift end-stamps back to interval starts. Without this the series
    # sits one step later than an observed record read from CSV, and the
    # comparison shows a systematic timing error that is purely conventional. ---
    if stamp == "end" and len(series) > 1:
        step = series.index[1] - series.index[0]
        series.index = series.index - step
        log.info("Shifted the index back by one %s interval to interval starts",
                 step)

    if tz:
        series.index = series.index.tz_localize(tz)

    log.info("%s: %d values, %s to %s, units %r",
             parts[1], len(series), series.index.min(), series.index.max(), units)
    finite = series.dropna()
    if finite.empty:
        raise ValueError(f"{pathname} contains no finite values — the compute "
                         "almost certainly failed despite a clean exit")
    log.info("  min %.3f, max %.3f, mean %.3f", finite.min(), finite.max(),
             finite.mean())
    return series


def read_run_results(
    dss_path: str,
    elements: list[str],
    run_label: str,
    basin: str,
    parameter: str = "FLOW",
    interval: str = "15MIN",
) -> pd.DataFrame:
    """
    Read one parameter for many elements of a single run into one frame.

    Asserts that every requested element produced records, which is the only
    reliable success test for a HEC-HMS compute.
    """
    frames, missing = {}, []
    for element in elements:
        pathname = f"/{basin}/{element}/{parameter}//{interval}/{run_label}/"
        try:
            frames[element] = read_dss_series(dss_path, pathname)
        except (KeyError, ValueError) as exc:
            log.error("%s: %s", element, exc)
            missing.append(element)

    if missing:
        raise RuntimeError(
            f"{len(missing)} of {len(elements)} element(s) produced no results: "
            f"{missing[:5]} — the run did not complete, whatever its exit code")

    df = pd.DataFrame(frames)
    log.info("Assembled %d element(s) × %d timesteps", df.shape[1], df.shape[0])

    # A run whose elements have different lengths means the compute stopped
    # partway; the frame will be full of NaN in the tail rather than short.
    tail_nan = df.tail(1).isna().sum().sum()
    if tail_nan:
        log.error("%d element(s) have no value at the final timestep — the "
                  "compute terminated early", int(tail_nan))
    return df


def catalogue(dss_path: str, contains: str = "") -> list[str]:
    """List the pathnames in a DSS file — the first thing to run when a read fails."""
    from hecdss import HecDss
    with HecDss(dss_path) as dss:
        paths = [p for p in dss.get_catalog() if contains in p]
    log.info("%d pathname(s) in %s matching %r", len(paths), dss_path, contains)
    for p in paths[:10]:
        log.info("  %s", p)
    return paths


# --- Example usage ---
# catalogue("run.dss", contains="FLOW")
# df = read_run_results("run.dss", ["SUBBASIN-1", "SUBBASIN-2", "OUTLET"],
#                       run_label="RUN:100YR", basin="BASIN")

WARNING Pathname carries a D part ('01JAN2024') — blanking it
ERROR   No record at /BASIN/SUBBASIN-3/FLOW//15MIN/RUN:100YR/. 2 similar
        pathname(s) exist: ['/BASIN/SUBBASIN 3/FLOW//15MIN/RUN:100YR/',
        '/BASIN/SUBBASIN-3/FLOW//1HOUR/RUN:100YR/']