Computing 7Q10 and Low-Flow Statistics in Python

import logging

import numpy as np
import pandas as pd
from scipy import stats

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


def climatic_year(index: pd.DatetimeIndex, start_month: int = 4) -> pd.Series:
    """
    Label each date with its climatic year.

    A climatic year starting in April places the whole northern-hemisphere
    low-flow season inside one year, so a drought spanning the calendar
    autumn is not split across two annual minima.
    """
    return pd.Series(
        index.year - (index.month < start_month).astype(int), index=index
    )


def low_flow_statistics(
    daily_cms: pd.Series,
    window_days: int = 7,
    return_periods=(2, 10, 20),
    climatic_start_month: int = 4,
    min_coverage: float = 0.90,
) -> dict:
    """
    Compute nQm low-flow statistics from a daily discharge record.

    Parameters
    ----------
    daily_cms            : Daily mean discharge, indexed by date.
    window_days          : Averaging window; 7 gives the familiar 7Qm family.
    return_periods       : Recurrence intervals to report.
    climatic_start_month : 4 for the northern hemisphere, 10 for the southern.
    min_coverage         : Fraction of a year that must be present for its
                           minimum to be trusted.
    """
    s = pd.Series(daily_cms).dropna().sort_index()
    if s.empty:
        raise ValueError("empty discharge series")

    # --- The rolling mean must require a full window. min_periods=1 would
    # produce a "seven-day" mean from a single day at the start of a gap,
    # which lands as a spuriously low annual minimum. ---
    rolled = s.rolling(f"{window_days}D", min_periods=window_days).mean()

    cy = climatic_year(s.index, climatic_start_month)
    frame = pd.DataFrame({"rolled": rolled, "cy": cy.to_numpy()})

    annual_min, dropped = {}, []
    for year, grp in frame.groupby("cy"):
        coverage = grp["rolled"].notna().sum() / 365.0
        if coverage < min_coverage:
            dropped.append((int(year), round(float(coverage), 2)))
            continue
        annual_min[int(year)] = float(grp["rolled"].min())

    if dropped:
        log.warning("Dropped %d year(s) below %.0f %% coverage: %s",
                    len(dropped), 100 * min_coverage, dropped[:6])

    minima = pd.Series(annual_min).sort_index()
    n = len(minima)
    if n < 10:
        raise ValueError(f"{n} usable climatic years is too short for a "
                         "low-flow frequency estimate")

    n_zero = int((minima <= 0).sum())
    nonzero = minima[minima > 0]
    p_zero = n_zero / n

    if n_zero:
        log.warning("%d of %d years had a zero %d-day minimum (%.1f %%) — "
                    "using conditional probability rather than deleting them",
                    n_zero, n, window_days, 100 * p_zero)

    logs = np.log10(nonzero.to_numpy(dtype=float))
    mean_log, std_log = float(logs.mean()), float(logs.std(ddof=1))
    skew = float(stats.skew(logs, bias=False))

    results = {}
    for T in return_periods:
        # For LOW flow the target is the NON-exceedance probability: the value
        # undercut once in T years. This is the sign flip that separates a
        # low-flow fit from a flood fit.
        p_nonexceed = 1.0 / T
        if n_zero:
            # Conditional adjustment: the fitted distribution describes only
            # the non-zero years, so rescale the probability into that subset.
            if p_nonexceed <= p_zero:
                results[f"{window_days}Q{T}"] = 0.0
                log.info("%dQ%d = 0.0 m³/s — zero flow is more frequent than "
                         "1 in %d years", window_days, T, T)
                continue
            p_nonexceed = (p_nonexceed - p_zero) / (1.0 - p_zero)

        z = stats.norm.ppf(p_nonexceed)
        if abs(skew) < 1e-6:
            kt = z
        else:
            k = skew / 6.0
            kt = (2.0 / skew) * (((z - k) * k + 1.0) ** 3 - 1.0)
        value = float(10 ** (mean_log + kt * std_log))
        results[f"{window_days}Q{T}"] = value
        log.info("%dQ%d = %.4f m³/s", window_days, T, value)

    # --- Flow-duration percentiles are a different statistic entirely, but
    # they are cheap here and worth reporting alongside for context. ---
    fdc = {f"Q{p}": float(np.percentile(s, 100 - p)) for p in (50, 75, 90, 95, 99)}

    log.info("Record: %d climatic years %d–%d, %d zero-flow years",
             n, int(minima.index.min()), int(minima.index.max()), n_zero)
    log.info("Flow duration: Q50 %.3f, Q95 %.4f m³/s", fdc["Q50"], fdc["Q95"])

    return {
        "statistics": results,
        "flow_duration": fdc,
        "annual_minima": minima,
        "n_years": n,
        "n_zero_years": n_zero,
        "dropped_years": dropped,
    }


# --- Example usage ---
# out = low_flow_statistics(daily_series, window_days=7,
#                           return_periods=(2, 10, 20))
# print(out["statistics"])