#!/usr/bin/env python3
"""Estonia's thermal growing season from public station data, 1965-2025 and 1926-2025.

Our re-run for https://factchecked.science/estonia-growing-season/ (Fact-Checked Science, October 26, 2026).
It is not the paper's analysis: by its abstract, Saue (2026, International Journal of Climatology,
doi:10.1002/joc.70638) used daily mean air temperatures from 15 Estonian meteorological stations, probably the
Estonian Environment Agency's, which the author's other 2026 data sets say are available on request (our
inference; we could not read the paper's text). This script uses NOAA's GHCN-Daily, which carries daily maximum
and minimum temperatures for Estonian stations, mostly from the European Climate Assessment & Dataset (ECA&D;
Tartu-Toravere's until 1964 are flagged as government-supplied data); true daily means (TAVG) are nearly
complete there only at Tallinn and Tartu-Toravere, so the daily mean is approximated as TG = (TX + TN) / 2 and
checked against TAVG at those two stations: day by day, and with TAVG shifted uniformly by +0.1 C.

Definition (ETCCDI, https://etccdi.pacificclimate.org/list_27_indices.shtml): the season starts on the first
day of the first run of at least 6 days with TG > 5 C (searched 1 January to 30 June) and ends on the first
day of the first run of at least 6 days with TG < 5 C after 1 July (1 January of the next year if there is
none); its length is end minus start, in days. A missing day breaks a run. A year counts only if it has at
most 15 missing days and no month with more than 3 (the climdex rule). Values with a quality flag are dropped.
Trends: Theil-Sen slope and a two-sided Mann-Kendall test (normal approximation, tie-corrected). The network
series is the mean, year by year, of the stations' anomalies from their own 1965-2025 means, for years with at
least 5 of the 11 stations, plus the stations' average mean.

Usage:
    python3 gsl_rerun.py                     # downloads 12 station files (about 6 MB) into ./ghcnd/
    python3 gsl_rerun.py --data DIR          # use files already in DIR (named <station>.csv.gz)
    python3 gsl_rerun.py --out numbers.json  # where to write the results (default gsl-rerun-numbers.json)

Standard library only (Python 3.8 or later). NOAA updates the station files daily; the script prints each
file's SHA-256 and says whether it matches the file we used on October 26, 2026. Results for 1965-2025 should
not change unless NOAA revises past values or quality flags.

Data credit: Menne et al. (2012), doi:10.1175/JTECH-D-11-00103.1, and NOAA NCEI's GHCN-Daily
(doi:10.7289/V5D21VHZ); the maximum and minimum series come mostly from ECA&D (Klein Tank et al. 2002,
doi:10.1002/joc.773). NOAA's GHCN-Daily documentation says some non-US data are exchanged under WMO Resolution 40
(Cg-XII) and adds: "For non-U.S. locations data, the data or any derived product shall not be provided to other
users or be used for the re-export of commercial services."
"""
import argparse
import calendar
import collections
import csv
import datetime
import gzip
import hashlib
import json
import math
import os
import statistics
import sys
import urllib.request

URL = "https://www.ncei.noaa.gov/pub/data/ghcn/daily/by_station/{}.csv.gz"

# The 11 stations with near-complete ECA&D maximum and minimum temperatures for 1965-2025
STATIONS = {
    "EN000026045": "Kunda", "EN000026242": "Tartu-Tõravere", "EN000026249": "Võru",
    "ENE00100914": "Vilsandi", "ENE00175011": "Kuressaare", "ENE00175023": "Kuusiku",
    "ENE00175027": "Narva", "ENE00175035": "Pärnu", "ENE00175051": "Tallinn",
    "ENE00175079": "Türi", "ENE00175087": "Valga",
}
# Seven of them with long records, for 1926-2025
CENTURY = ["EN000026242", "ENE00175035", "ENE00175051", "ENE00175027", "ENE00100914", "EN000026249", "ENE00175023"]
# Daily means (TAVG) for the check of the approximation: Tallinn's TAVG is filed under the WMO number
TAVG_CHECK = [("Tartu-Tõravere", "EN000026242", "EN000026242"), ("Tallinn", "EN000026038", "ENE00175051")]

# The files we used, fetched October 26, 2026 (UTC 19:22)
OUR_RUN = {
    "EN000026045": "f3e10ce84c18657dc5b917a45662210e2e67425d0bc1155f280bb86751bb00b2",
    "EN000026242": "beb9be0bc8c273617d13e669b1999119d125565942c93ee438fc99abfb69e4fd",
    "EN000026249": "8e0055c6fe1c1d204697b60b59d2f1bf78a4604552f023b0e1e1459458cae6c0",
    "ENE00100914": "55acbfa31ca4b51b4e9fd27d44d5e75869653ac5353ea6328f59b528d9af95eb",
    "ENE00175011": "0cb9a162f72a7f5ae39127c0b9f72df2c7fa095d62186f4f0494d0674e100a90",
    "ENE00175023": "8fd0142feed3ab1466374db491c2e69e0a61c95c23646f5ec120145ee56ad929",
    "ENE00175027": "d0df2ef86cac62ecd8e909529af2516ec8463925b8c1d1fcd48e4eb7f513b1cc",
    "ENE00175035": "7da96f2273e975e09ccec3c59613ff92e44af1bce74f440db10a760af3971a4b",
    "ENE00175051": "810abecfa9f29a907bc4fd303344150326405e8f1bdaacd94e76ab61b34f3369",
    "ENE00175079": "5c2cc7147f20d4dd84071d50b0c026f9e51eb496ab2680d07afd482b156d1cac",
    "ENE00175087": "c9b27ff53a78345f4ccc2f5cb4173abd1a367b4a7cffd5a46ca678110874501a",
    "EN000026038": "0b709e115744f33a57d1fc5b8060d33e4e0e663805f0e2b4bd2299ff7ca3d75c",
}
ONE = datetime.timedelta(days=1)


# ---------- data ----------
def fetch(sid, data_dir, download=True):
    path = os.path.join(data_dir, sid + ".csv.gz")
    if not os.path.exists(path):
        if not download:
            sys.exit(f"missing {path}")
        os.makedirs(data_dir, exist_ok=True)
        req = urllib.request.Request(URL.format(sid), headers={"User-Agent": "gsl_rerun.py (factchecked.science)"})
        with urllib.request.urlopen(req, timeout=180) as r:
            body = r.read()
        with open(path, "wb") as f:
            f.write(body)
    with open(path, "rb") as f:
        sha = hashlib.sha256(f.read()).hexdigest()
    same = "same file as our run" if OUR_RUN.get(sid) == sha else "differs from our run of 26 Oct 2026 (NOAA updates these files)"
    print(f"  {sid}  sha256 {sha[:16]}…  {same}")
    return path


def read_station(path):
    """Return {element: {date: value in degrees C}} for TAVG, TMAX and TMIN, dropping quality-flagged values."""
    d = collections.defaultdict(dict)
    with gzip.open(path, "rt", encoding="ascii", errors="replace") as f:
        for row in csv.reader(f):
            if len(row) < 7:
                continue
            _, date, el, val, _mflag, qflag, _sflag = row[:7]
            if el not in ("TAVG", "TMAX", "TMIN") or qflag.strip():
                continue
            d[el][datetime.date(int(date[:4]), int(date[4:6]), int(date[6:8]))] = int(val) / 10.0
    return d


def tg_proxy(d):
    tx, tn = d.get("TMAX", {}), d.get("TMIN", {})
    return {day: (x + tn[day]) / 2.0 for day, x in tx.items() if day in tn and x >= tn[day]}


# ---------- the index ----------
def year_complete(series, y):
    missing = 0
    for m in range(1, 13):
        n = sum(1 for dd in range(1, calendar.monthrange(y, m)[1] + 1) if datetime.date(y, m, dd) not in series)
        if n > 3:
            return False
        missing += n
    return missing <= 15


def gsl_year(series, y, thr=5.0, run=6):
    """(start day of year, end day of year, length in days), or None if no start before 1 July."""
    jan1, jun30, jul1, dec31 = datetime.date(y, 1, 1), datetime.date(y, 6, 30), datetime.date(y, 7, 1), datetime.date(y, 12, 31)

    def warm(start):
        for k in range(run):
            v = series.get(start + k * ONE)
            if v is None or not v > thr:
                return False
        return True

    def cold(start):
        for k in range(run):
            day = start + k * ONE
            v = series.get(day)
            if day > dec31 or v is None or not v < thr:
                return False
        return True

    sos, day = None, jan1
    while day <= jun30:
        if warm(day):
            sos = day
            break
        day += ONE
    if sos is None:
        return None
    eos, day = None, jul1 + ONE
    while day <= dec31:
        if cold(day):
            eos = day
            break
        day += ONE
    if eos is None:
        eos = datetime.date(y + 1, 1, 1)
    return (sos - jan1).days + 1, (eos - jan1).days + 1, (eos - sos).days


def station_series(series, y0, y1):
    out = {}
    for y in range(y0, y1 + 1):
        if year_complete(series, y):
            g = gsl_year(series, y)
            if g is not None:
                out[y] = g
    return out


# ---------- statistics ----------
def theil_sen(xs, ys):
    return statistics.median((ys[j] - ys[i]) / (xs[j] - xs[i]) for i in range(len(xs)) for j in range(i + 1, len(xs)))


def mann_kendall_p(ys):
    n = len(ys)
    s = sum((ys[j] > ys[i]) - (ys[j] < ys[i]) for i in range(n - 1) for j in range(i + 1, n))
    var = (n * (n - 1) * (2 * n + 5) - sum(t * (t - 1) * (2 * t + 5) for t in collections.Counter(ys).values())) / 18.0
    z = 0.0 if s == 0 else (s - (1 if s > 0 else -1)) / math.sqrt(var)
    return math.erfc(abs(z) / math.sqrt(2))


def trend(res, idx, y0, y1):
    xs = [y for y in sorted(res) if y0 <= y <= y1]
    ys = [res[y][idx] for y in xs]
    return {"n": len(xs), "per_decade": round(10 * theil_sen(xs, ys), 2), "mk_p": round(mann_kendall_p(ys), 4),
            "mean": round(statistics.fmean(ys), 1)}


def network(per, idx, stations, y0, y1):
    means = {s: statistics.fmean(v[idx] for v in per[s].values()) for s in stations}
    grand = statistics.fmean(means.values())
    ser = {}
    for y in range(y0, y1 + 1):
        vals = [per[s][y][idx] - means[s] for s in stations if y in per[s]]
        if len(vals) >= max(3, len(stations) // 2):
            ser[y] = statistics.fmean(vals)
    return ser, grand


def slope_of(ser, y0, y1):
    xs = [y for y in sorted(ser) if y0 <= y <= y1]
    ys = [ser[y] for y in xs]
    return {"per_decade": round(10 * theil_sen(xs, ys), 2), "mk_p": round(mann_kendall_p(ys), 4), "n": len(xs)}


def day_label(doy):
    return (datetime.date(2001, 1, 1) + datetime.timedelta(days=round(doy) - 1)).strftime("%d %b")


# ---------- main ----------
def main():
    ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
    ap.add_argument("--data", default="ghcnd", help="folder for the station files (default ./ghcnd)")
    ap.add_argument("--out", default="gsl-rerun-numbers.json", help="results file (default gsl-rerun-numbers.json)")
    ap.add_argument("--no-download", action="store_true", help="fail instead of downloading missing files")
    a = ap.parse_args()
    if hasattr(sys.stdout, "reconfigure"):
        sys.stdout.reconfigure(encoding="utf-8")

    print("Station files (GHCN-Daily, NOAA NCEI):")
    ids = list(STATIONS) + ["EN000026038"]
    data = {sid: read_station(fetch(sid, a.data, not a.no_download)) for sid in ids}
    out = {"made_by": "gsl_rerun.py, Fact-Checked Science, https://factchecked.science/estonia-growing-season/",
           "definition": "ETCCDI growing season length, TG = (TX + TN) / 2 from GHCN-Daily (ECA&D maximum and minimum), 6-day runs, 5 C",
           "stations": {}, "network": {}, "proxy_check": {}, "century": {}}

    # 1965-2025
    print("\n1965-2025, 11 stations:")
    per = {}
    for sid, name in STATIONS.items():
        res = station_series(tg_proxy(data[sid]), 1965, 2025)
        per[sid] = res
        g, s, e = trend(res, 2, 1965, 2025), trend(res, 0, 1965, 2025), trend(res, 1, 1965, 2025)
        out["stations"][sid] = {"name": name, "years": len(res), "GSL": g, "SOS": s, "EOS": e,
                                "mean_start": day_label(s["mean"]), "mean_end": day_label(e["mean"]),
                                "series": {str(y): list(v) for y, v in sorted(res.items())}}
        print(f"  {name:15s} {len(res):2d} years  length {g['mean']:5.1f} d, {g['per_decade']:+5.2f} d/decade (p={g['mk_p']:.4f});"
              f"  start {day_label(s['mean'])} {s['per_decade']:+5.2f};  end {day_label(e['mean'])} {e['per_decade']:+5.2f}")
    for key, idx in (("GSL", 2), ("SOS", 0), ("EOS", 1)):
        sl = [out["stations"][s][key]["per_decade"] for s in STATIONS]
        out["network"][key + "_mean_of_station_slopes"] = round(statistics.fmean(sl), 2)
        out["network"][key + "_station_slope_range"] = [min(sl), max(sl)]
        out["network"][key + "_stations_p_below_0.05"] = sum(1 for s in STATIONS if out["stations"][s][key]["mk_p"] < 0.05)
        ser, grand = network(per, idx, list(STATIONS), 1965, 2025)
        out["network"][key + "_series"] = {str(y): round(v + grand, 2) for y, v in sorted(ser.items())}
        out["network"][key] = slope_of(ser, 1965, 2025)
        print(f"  {key}: mean of station slopes {statistics.fmean(sl):+.2f} d/decade; network series "
              f"{out['network'][key]['per_decade']:+.2f} d/decade (p={out['network'][key]['mk_p']:.4f}, {out['network'][key]['n']} years)")
    gs, grand = network(per, 2, list(STATIONS), 1965, 2025)
    out["network"]["GSL_change_1965_2025_days"] = round(out["network"]["GSL"]["per_decade"] * 6, 1)
    out["network"]["GSL_mean_1965_1994"] = round(statistics.fmean(v for y, v in gs.items() if y <= 1994) + grand, 1)
    out["network"]["GSL_mean_1996_2025"] = round(statistics.fmean(v for y, v in gs.items() if y >= 1996) + grand, 1)
    out["network"]["GSL_sd_year_to_year"] = round(statistics.pstdev(gs.values()), 1)
    out["network"]["GSL_mean"] = round(grand, 1)
    out["network"]["periods"] = {f"{y0}-{y1}": slope_of(gs, y0, y1) for (y0, y1) in
                                 [(1965, 2025), (1965, 2024), (1965, 2015), (1971, 2025), (1981, 2025), (1991, 2025),
                                  (1965, 1994), (1995, 2025), (1971, 2000), (1991, 2020)]}
    print(f"  30-year means: 1965-1994 {out['network']['GSL_mean_1965_1994']} d, 1996-2025 {out['network']['GSL_mean_1996_2025']} d;"
          f" year-to-year SD {out['network']['GSL_sd_year_to_year']} d")
    for p, v in out["network"]["periods"].items():
        print(f"    {p}: {v['per_decade']:+.2f} d/decade (p={v['mk_p']:.3f}, {v['n']} years)")

    # the approximation against true daily means
    print("\nTAVG against (TX+TN)/2:")
    for label, sid_tavg, sid_txtn in TAVG_CHECK:
        tav, prox = data[sid_tavg].get("TAVG", {}), tg_proxy(data[sid_txtn])
        r1, r2 = station_series(tav, 1965, 2025), station_series(prox, 1965, 2025)
        com = sorted(set(r1) & set(r2))
        diff = {k: round(statistics.fmean(r2[y][i] - r1[y][i] for y in com), 2) for k, i in (("SOS", 0), ("EOS", 1), ("GSL", 2))}
        rmse = math.sqrt(statistics.fmean((r2[y][2] - r1[y][2]) ** 2 for y in com))
        # day by day, April-October 1965-2025: the average gap and how far apart the two series are on a typical day
        days = [prox[d] - tav[d] for d in tav if d in prox and 1965 <= d.year <= 2025 and 4 <= d.month <= 10]
        bias = statistics.fmean(days)
        sd = statistics.pstdev(days)
        mad = statistics.fmean(abs(x) for x in days)
        share1 = sum(1 for x in days if abs(x) >= 1.0) / len(days)
        # a uniform shift of +0.1 C applied to TAVG, at the data's 0.1 C resolution (rounded, so that 4.9 + 0.1 is
        # exactly 5.0 and not a float a hair above the threshold): how much does that alone move the seasons?
        shifted = {d: round(v + 0.1, 1) for d, v in tav.items()}
        r3 = station_series(shifted, 1965, 2025)
        com3 = sorted(set(r1) & set(r3))
        rms_shift = math.sqrt(statistics.fmean((r3[y][2] - r1[y][2]) ** 2 for y in com3))
        changed = sum(1 for y in com3 if r3[y][2] != r1[y][2])
        t_tav = round(10 * theil_sen(com, [r1[y][2] for y in com]), 2)
        t_prox = round(10 * theil_sen(com, [r2[y][2] for y in com]), 2)
        out["proxy_check"][label] = {"years": len(com), "proxy_minus_tavg_days": diff, "gsl_rmse_days": round(rmse, 1),
                                     "daily_mean_bias_apr_oct_C": round(bias, 2), "daily_sd_apr_oct_C": round(sd, 2),
                                     "daily_mean_abs_diff_apr_oct_C": round(mad, 2), "share_days_1C_or_more": round(share1, 3),
                                     "days_apr_oct": len(days),
                                     "uniform_plus_0.1C_shift": {"years": len(com3), "gsl_rms_change_days": round(rms_shift, 1),
                                                                 "years_changed": changed},
                                     "gsl_trend_tavg": t_tav, "gsl_trend_proxy": t_prox}
        print(f"  {label}: {len(com)} years; start {diff['SOS']:+.1f} d, end {diff['EOS']:+.1f} d, length {diff['GSL']:+.1f} d;"
              f" RMSE {rmse:.1f} d; trend {t_tav:+.2f} (TAVG) vs {t_prox:+.2f} (proxy) d/decade")
        print(f"     daily, Apr-Oct ({len(days)} days): mean {bias:+.2f} C, SD {sd:.2f} C, mean absolute {mad:.2f} C,"
              f" {100 * share1:.0f}% of days 1 C or more apart")
        print(f"     TAVG shifted uniformly by +0.1 C: length changes {rms_shift:.1f} d RMS, in {changed} of {len(com3)} years")

    # 1926-2025
    print("\n1926-2025, 7 long records:")
    cper = {}
    for sid in CENTURY:
        res = station_series(tg_proxy(data[sid]), 1901, 2025)
        cper[sid] = res
        t = trend(res, 2, 1926, 2025)
        out["century"][sid] = {"name": STATIONS[sid], "GSL_1926_2025": t,
                               "series": {str(y): list(v) for y, v in sorted(res.items())}}
        print(f"  {STATIONS[sid]:15s} {t['n']} years: {t['per_decade']:+.2f} d/decade (p={t['mk_p']:.4f})")
    cs, cgrand = network(cper, 2, CENTURY, 1926, 2025)
    cn = slope_of(cs, 1926, 2025)
    cn.update({k: round(statistics.fmean(v for y, v in cs.items() if lo <= y <= hi) + cgrand, 1)
               for k, (lo, hi) in {"mean_1926_1955": (1926, 1955), "mean_1956_1985": (1956, 1985), "mean_1996_2025": (1996, 2025)}.items()})
    cn["series"] = {str(y): round(v + cgrand, 2) for y, v in sorted(cs.items())}
    out["century"]["network"] = cn
    print(f"  network: {cn['per_decade']:+.2f} d/decade (p={cn['mk_p']:.4f}, {cn['n']} years); 30-year means "
          f"{cn['mean_1926_1955']}, {cn['mean_1956_1985']}, {cn['mean_1996_2025']} d")

    with open(a.out, "w", encoding="utf-8", newline="\n") as f:
        json.dump(out, f, ensure_ascii=False, indent=1)
    print(f"\nwrote {a.out}")


if __name__ == "__main__":
    main()