#!/usr/bin/env python3
"""A rough flood-frequency estimate for the San Juan River at Pagosa Springs (USGS 09342500).
Fact-Checked Science, "S0 News, November 1, 2026", section 10 ("Our analysis"). Standard library only.

Run it
  python3 -I pagosa_flood_freq.py [folder]
  The folder (default: the script's own) must hold, or will receive:
    usgs_peak_09342500.rdb  https://nwis.waterdata.usgs.gov/nwis/peak?site_no=09342500&agency_cd=USGS&format=rdb
  and must hold the extract published beside this page:
    pagosa-stage-flow-2026-10.csv  USGS's preliminary instantaneous stage and flow at or above 10 ft, October 22-25, 2026,
                                   read November 1, 2026. Preliminary values can be revised. USGS marks 01:00-05:15 MDT on
                                   October 23 as an equipment malfunction, so the maximum, at 05:30, may understate the crest.

Rules, fixed before running
  - The annual peak series is the USGS file, which ends with water year 2025 (a water year starts on October 1). We run it as it
    stands, then with October 2025's 8,570 cfs added as water year 2026 (the peak the Western Water Assessment gives; USGS's
    file doesn't hold it yet), then with October 2026's preliminary peak added as well, as water year 2027.
  - Log-Pearson III by the method of moments on log10 Q with station skew, frequency factors from the Wilson-Hilferty
    approximation. This is a screening number, without the historical-period and low-outlier treatment of Bulletin 17C. The
    Western Water Assessment's USGS PeakFQ run put 8,570 cfs (the peak of October 2025) at 25-40 years and 25,000 cfs (1911) at about 200;
    this fit agrees at the first and gives about ten times as long at the second, so it says little about flows well above 9,000 cfs.
  - A straight line through ln Q against stage, for stages at or above 10 ft, reads what NWS's forecast crest of 17.0 ft would
    be. It extrapolates beyond the highest measured pair (15.02 ft), so it is only a rough reading.
"""
import csv
import math
import os
import sys
import urllib.request
from statistics import NormalDist

PEAK_URL = "https://nwis.waterdata.usgs.gov/nwis/peak?site_no=09342500&agency_cd=USGS&format=rdb"


def fit(q):
    x = [math.log10(v) for v in q]
    n = len(x)
    m = sum(x) / n
    s = math.sqrt(sum((a - m) ** 2 for a in x) / (n - 1))
    g = n / ((n - 1) * (n - 2)) * sum(((a - m) / s) ** 3 for a in x)
    return m, s, g, n


def kt(T, g):
    z = NormalDist().inv_cdf(1 - 1.0 / T)
    if abs(g) < 1e-6:
        return z
    k = g / 6
    return (2 / g) * ((1 + k * z - k * k) ** 3 - 1)


def return_period(Q, m, s, g):
    lo, hi = 1.01, 100000.0
    for _ in range(200):
        mid = math.sqrt(lo * hi)
        if 10 ** (m + kt(mid, g) * s) < Q:
            lo = mid
        else:
            hi = mid
    return mid


def main():
    sys.stdout.reconfigure(encoding="utf-8")
    folder = sys.argv[1] if len(sys.argv) > 1 else os.path.dirname(os.path.abspath(__file__))
    path = os.path.join(folder, "usgs_peak_09342500.rdb")
    if not os.path.exists(path):
        print("downloading", PEAK_URL)
        req = urllib.request.Request(PEAK_URL, headers={"User-Agent": "factchecked.science pagosa_flood_freq.py"})
        with urllib.request.urlopen(req, timeout=120) as r, open(path, "wb") as f:
            f.write(r.read())
    pk = []
    with open(path, encoding="utf-8", errors="replace") as f:
        for line in f:
            if line.startswith("USGS"):
                r = line.split("\t")
                if r[4].strip():
                    pk.append((r[2], float(r[4])))
    with open(os.path.join(folder, "pagosa-stage-flow-2026-10.csv"), encoding="utf-8") as f:
        iv = [(r["time_mdt"], float(r["flow_cfs"]), float(r["stage_ft"])) for r in csv.DictReader(f)]
    top = max(iv, key=lambda x: x[1])
    print(f"Annual peaks in the USGS file: {len(pk)}, {pk[0][0]} to {pk[-1][0]}")
    print("Top peaks:", sorted(pk, key=lambda x: -x[1])[:6])
    print("Preliminary maximum, October 2026:", top)
    base = [v for _, v in pk]
    for name, q in (("Peaks to water year 2025", base), ("Plus October 2025's 8,570 cfs (water year 2026)", base + [8570.0]),
                    ("Plus October 2026's preliminary peak as well (water year 2027)", base + [8570.0, top[1]])):
        m, s, g, n = fit(q)
        print(f"\n{name}: n={n}, mean log10={m:.3f}, sd={s:.3f}, skew={g:.3f}")
        for T in (10, 25, 50, 100, 200, 500):
            print(f"  Q{T} = {10 ** (m + kt(T, g) * s):,.0f} cfs")
        for Q in (8570, 9230, 13660, 25000):
            print(f"  {Q:,} cfs: about a {return_period(Q, m, s, g):,.0f}-year return period")
    pairs = sorted({(g, q) for _, q, g in iv})
    xs = [g for g, _ in pairs]
    ys = [math.log(q) for _, q in pairs]
    xb, yb = sum(xs) / len(xs), sum(ys) / len(ys)
    b = sum((a - xb) * (c - yb) for a, c in zip(xs, ys)) / sum((a - xb) ** 2 for a in xs)
    a0 = yb - b * xb
    print(f"\nStage-discharge pairs: {len(pairs)}, {pairs[0]} to {pairs[-1]}")
    print(f"ln Q = {a0:.3f} + {b:.4f} * stage; 15.02 ft reads {math.exp(a0 + b * 15.02):,.0f} cfs; 17.0 ft reads {math.exp(a0 + b * 17.0):,.0f} cfs (extrapolated)")
    print("\nBetween the Western Water Assessment's figures (8,570 cfs at 25-40 years, 25,000 cfs at about 200), log on log; rough:")
    for Q in (9230, 13660):
        print(f"  {Q:,} cfs: {' to '.join(f'{T0 * (Q / 8570) ** (math.log(200 / T0) / math.log(25000 / 8570)):.0f}' for T0 in (25, 40))} years")


if __name__ == "__main__":
    main()