#!/usr/bin/env python3
"""Station check: do geomagnetic storms leave a mark on weather-station records?

The code behind section 3, "A check with public data", of "Do geomagnetic storms change the
weather?" on Fact-Checked Science: https://factchecked.science/storms-and-weather/

It asks whether the strong and moderate winter geomagnetic storms of 1958-2024 coincide with
temperature or precipitation anomalies at 10 North American weather stations more often than
randomly timed storms do, and repeats the precipitation test for strong storms in every season.
It is a plausibility test of claims in J. Raeder (2026), "Regional and Seasonal Effects of
Geomagnetic Storms on Terrestrial Weather", Geophysical Research Letters 53(10), e2025GL121097
(https://doi.org/10.1029/2025GL121097), not a replication of that paper.

Run it in an empty folder:

    python3 station_check.py
    python3 station_check.py --cache <download folder> --out <output folder>

By default it downloads into ./station_check_data and writes into ./station_check_results, both
in the folder you run it from. It needs Python 3 with NumPy and pandas. The first run downloads
about 200 MB: NASA's OMNI2 hourly file (184 MB), 72 monthly pages of Kyoto Dst (under 1 MB) and
10 GHCN-Daily station files (about 13 MB). Later runs reuse the downloads and take about a minute.

Data (all public):
    Hourly Dst   WDC for Geomagnetism, Kyoto, final values, 1957-1962 (doi:10.17593/14515-74000)
                 NASA SPDF OMNI2 hourly data set, 1963-2024 (doi:10.48322/1shr-ht18); OMNI takes
                 its Dst from Kyoto (final values to 2020, provisional after)
    Stations     NOAA GHCN-Daily, by-station files (doi:10.7289/V5D21VHZ)
Both are updated over time, so a later run can differ slightly from ours, which used GHCN-Daily
3.35-upd-2026101017 and the OMNI2 file as downloaded on October 11, 2026.

Outputs (in --out):
    station-check-results.json  every number in section 3 of the page, with data versions, settings
                                and sensitivity checks
    storms.csv                  every storm used (both levels, all seasons, 1958-2024)
    station_tests.csv           the 40 winter station tests (10 stations x 2 variables x 2 levels)

Method (every random generator has a fixed seed and the order of the draws is fixed, so a run on
the same data reproduces our numbers exactly):

 1. Hourly Dst. WDC for Geomagnetism, Kyoto, final Dst monthly pages for 1957-1962
    (https://wdc.kugi.kyoto-u.ac.jp/dst_final/YYYYMM/index.html) and NASA SPDF OMNI2
    hourly file for 1963-2024 (word 41 = Dst; OMNI takes Dst from Kyoto).
 2. Storm days. A UT day is a strong-storm day if its lowest hourly Dst is -200 nT or
    below, and a moderate-storm day if its lowest Dst is below -100 and above -200 nT
    (-200 < min <= -100). Days of one level in one season (winter = Dec-Feb, spring =
    Mar-May, summer = Jun-Aug, fall = Sep-Nov) are merged into one storm when they are at
    most 2 days apart, so separate storms have at least two days between them that are not
    days of that level. "Hours" are the hours with Dst in the level's range.
 3. Station data. NOAA GHCN-Daily by-station files; TMAX, TMIN, PRCP values with a blank
    quality flag; daily mean temperature T = (TMAX + TMIN) / 2; years 1958-2024.
    (GHCN-Daily v3.35, released 2026-10-01, renamed Environment Canada stations from
    CA00xxxxxxx to CAN0xxxxxxx; the script tries the new ID first, then the old one.)
 4. Anomalies. For each station and variable: remove a straight-line trend fitted to all
    days, then subtract a 31-day running mean of the day-of-year mean and divide by a
    31-day running mean of the day-of-year standard deviation (so values are in SD units;
    the page shows them x100 as "% of SD"). Before the division, the departure is in the
    variable's own units (K for temperature); the per-day tables report both.
 5. Winter station tests (seed 1). For each station, variable and level, the real value is the
    mean anomaly over all winter storm days with data. 3,000 random sets are drawn, each with
    the same number of storms and the same storm lengths (in days), each pseudo-storm starting
    on a random winter day with data. The plotted range is the 2.5-97.5 percentile range of
    the random sets. p is two-tailed and percentile-based: twice the share of random sets at
    or beyond the real mean on its side, min(1, 2 x min(share >= real, share <= real)).
    p_from_zero, the share of random sets whose mean is at least as far from zero as the real
    one, is kept for comparison: random precipitation means centre below zero, so p_from_zero
    cannot register unusual dryness. For strong-storm precipitation the script also scores a
    completely dry set of storm days (0 mm on every storm day) against the same random sets.
 6. Coherence and precipitation sign (seed 11, continuing): 5,000 random five-storm sets
    (lengths as the real strong winter storms) over all ten stations at once. Then, per station,
    20,000 draws of 5 random winter days for the precipitation skew. Then the lag composite:
    station-mean temperature anomaly from 10 days before to 10 days after the first day of each
    of the 77 moderate winter storms, against 2,000 sets of 77 random winter onset days.
 7. Strong storms and precipitation in every season (seeds 3 and 13): the precipitation test of
    step 5 for strong storms in each season, and the number of stations drier than average
    against 5,000 random sets of storms, as in step 6.
 8. Winter as +-30 days around 21 December, the season window of the paper author's 2024
    conference poster (a separate generator, seed 1): the 40 station tests of step 5 again, with
    the storms and the random start days taken from that window.
 9. Other sensitivity checks: other storm definitions and other year windows.

The script uses only NumPy and pandas besides the standard library.
"""
import argparse
import datetime as dt
import json
import os
import ssl
import urllib.error
import urllib.request
import warnings

import numpy as np
import pandas as pd

# ---------------------------------------------------------------------------------------
# Configuration
# ---------------------------------------------------------------------------------------
YEARS = (1958, 2024)
SEASONS = {"winter": [12, 1, 2], "spring": [3, 4, 5], "summer": [6, 7, 8], "fall": [9, 10, 11]}
# Level ranges on the daily minimum Dst: lo < min <= hi
LEVELS = {"moderate": (-200.0, -100.0), "strong": (-9999.0, -200.0)}

# Station order matters: it fixes the order of the random draws (as in the original code).
# (current GHCN-Daily ID, ID used on 2026-09-30, display name)
STATIONS = [
    ("USW00014942", "USW00014942", "Omaha"),
    ("CAN05023222", "CA005023222", "Winnipeg"),
    ("USW00014922", "USW00014922", "Minneapolis"),
    ("USW00023062", "USW00023062", "Denver"),
    ("USW00094846", "USW00094846", "Chicago"),
    ("USW00003947", "USW00003947", "Kansas City"),
    ("CAN05060600", "CA005060600", "Churchill, Hudson Bay"),
    ("CAN03031093", "CA003031093", "Calgary"),
    ("USW00003927", "USW00003927", "Dallas–Fort Worth"),
    ("USW00024233", "USW00024233", "Seattle"),
]
N_TEST = 3000      # random sets per station test
N_COHERENCE = 5000 # random sets of storms for the coherence and precipitation-sign figures
N_SKEW = 20000     # random 5-day draws per station for the precipitation skew
N_LAG = 2000       # random onset sets for the lag composite
SEED_TESTS, SEED_FIELD = 1, 11           # steps 5 and 6
SEED_SEASON_TESTS, SEED_SEASON_SIGN = 3, 13  # step 7
SEED_WINDOW = 1                          # step 8 (a separate generator)

OMNI_URL = "https://spdf.gsfc.nasa.gov/pub/data/omni/low_res_omni/omni2_all_years.dat"
KYOTO_URL = "https://wdc.kugi.kyoto-u.ac.jp/dst_final/{y}{m:02d}/index.html"
GHCN_URL = "https://www.ncei.noaa.gov/pub/data/ghcn/daily/by_station/{sid}.csv.gz"
GHCN_VERSION_URL = "https://www.ncei.noaa.gov/pub/data/ghcn/daily/ghcnd-version.txt"


# ---------------------------------------------------------------------------------------
# Downloads (cached)
# ---------------------------------------------------------------------------------------
def fetch(url, path):
    """Download url to path unless it is already there. Returns path, or None on HTTP 404."""
    if os.path.exists(path) and os.path.getsize(path) > 0:
        return path
    os.makedirs(os.path.dirname(path), exist_ok=True)
    print(f"Downloading {url}", flush=True)
    ctx = ssl.create_default_context()
    try:
        with urllib.request.urlopen(url, timeout=600, context=ctx) as r:
            data = r.read()
    except urllib.error.HTTPError as e:
        if e.code == 404:
            return None
        raise
    with open(path + ".part", "wb") as f:
        f.write(data)
    os.replace(path + ".part", path)
    return path


# ---------------------------------------------------------------------------------------
# 1. Hourly Dst
# ---------------------------------------------------------------------------------------
def read_kyoto_month(path, year, month):
    """Parse one Kyoto 'dst_final' monthly page: rows of day + 24 hourly values (4 chars each,
    three groups of 8 separated by a space). Hour index 0 is 00-01 UT. 9999 = missing."""
    txt = open(path, "rb").read().decode("euc-jp", "ignore")
    rows, inpre = [], False
    for line in txt.splitlines():
        if line.startswith("DAY"):
            inpre = True
            continue
        if not inpre:
            continue
        if line.startswith("<"):
            break
        if not line.strip():
            continue
        day, rest = int(line[:3]), line[3:]
        vals = []
        for g in range(3):
            seg = rest[g * 33:g * 33 + 32]
            for k in range(8):
                v = seg[k * 4:k * 4 + 4].strip()
                vals.append(float(v) if v not in ("", "9999") else np.nan)
        for hr, v in enumerate(vals):
            rows.append((dt.datetime(year, month, day, hr), v))
    return rows


def load_dst(cache):
    """Hourly Dst 1958-2024 (Kyoto 1958-1962, OMNI 1963-2024) as a Series indexed by UT hour."""
    rows = []
    for y in range(1957, 1963):
        for m in range(1, 13):
            p = fetch(KYOTO_URL.format(y=y, m=m), os.path.join(cache, "kyoto", f"{y}{m:02d}.html"))
            rows += read_kyoto_month(p, y, m)
    ky = pd.Series([v for _, v in rows], index=pd.DatetimeIndex([t for t, _ in rows]))
    omni = fetch(OMNI_URL, os.path.join(cache, "omni", "omni2_all_years.dat"))
    w = pd.read_csv(omni, sep=r"\s+", header=None, usecols=[0, 1, 2, 40], names=["y", "doy", "h", "dst"])
    t = pd.to_datetime(w.y.astype(str) + "-" + w.doy.astype(str), format="%Y-%j") + pd.to_timedelta(w.h, unit="h")
    om = pd.Series(np.where(w.dst < 99999, w.dst, np.nan).astype(float), index=pd.DatetimeIndex(t))
    dst = pd.concat([ky[ky.index.year <= 1962], om[om.index.year >= 1963]]).sort_index()
    dst = dst[(dst.index.year >= YEARS[0]) & (dst.index.year <= YEARS[1])]
    return dst, ky, om


# ---------------------------------------------------------------------------------------
# 2. Storm days and storms
# ---------------------------------------------------------------------------------------
def level_days(dmin, level):
    lo, hi = LEVELS[level]
    return sorted(dmin.index[(dmin > lo) & (dmin <= hi)])


def group_days(days, gap=2):
    """Merge sorted days into storms when consecutive days are at most `gap` days apart."""
    ev = []
    for x in sorted(days):
        if ev and (x - ev[-1][-1]).days <= gap:
            ev[-1].append(x)
        else:
            ev.append([x])
    return ev


def storms_for(dmin, level, months=None, keep=None):
    days = level_days(dmin, level)
    if months is not None:
        days = [x for x in days if x.month in months]
    if keep is not None:
        days = [x for x in days if keep(x)]
    return group_days(days)


def near_dec21(x):
    """True for days within 30 days of 21 December (the author's 2024 poster's winter window)."""
    return any(abs((x - pd.Timestamp(yy, 12, 21)).days) <= 30 for yy in (x.year - 1, x.year))


def storm_table(dst, dmin):
    """One row per storm (both levels, all seasons) for storms.csv."""
    rows = []
    day_of_hour = dst.index.normalize()
    for level, (lo, hi) in LEVELS.items():
        for season, months in SEASONS.items():
            for e in storms_for(dmin, level, months):
                h = dst[day_of_hour.isin(e)]
                inlev = h[(h > lo) & (h <= hi)]
                rows.append({
                    "level": level, "season": season,
                    "first_day": e[0].date().isoformat(), "last_day": e[-1].date().isoformat(),
                    "storm_days": len(e),
                    "start_utc": inlev.index.min().strftime("%Y-%m-%d %H:00"),
                    "end_utc": (inlev.index.max() + pd.Timedelta(hours=1)).strftime("%Y-%m-%d %H:00"),
                    "hours_in_level": int(len(inlev)),
                    "min_dst_nT": int(h.min()),
                    "min_dst_hour_utc": h.idxmin().strftime("%Y-%m-%d %H:00"),
                })
    return pd.DataFrame(rows).sort_values(["first_day", "level"]).reset_index(drop=True)


# ---------------------------------------------------------------------------------------
# 3-4. Station data and anomalies
# ---------------------------------------------------------------------------------------
def load_station(cache, sid_new, sid_old):
    path = fetch(GHCN_URL.format(sid=sid_new), os.path.join(cache, "ghcn", sid_new + ".csv.gz"))
    used = sid_new
    if path is None:
        path = fetch(GHCN_URL.format(sid=sid_old), os.path.join(cache, "ghcn", sid_old + ".csv.gz"))
        used = sid_old
    df = pd.read_csv(path, header=None, names=["id", "date", "el", "val", "m", "q", "s", "t"],
                     dtype={"date": str, "m": str, "q": str, "s": str, "t": str}, low_memory=False)
    df = df[df.q.isna() & df.el.isin(["TMAX", "TMIN", "PRCP"])]
    p = df.pivot_table(index="date", columns="el", values="val", aggfunc="first")
    p.index = pd.to_datetime(p.index, format="%Y%m%d")
    p = p[(p.index.year >= YEARS[0]) & (p.index.year <= YEARS[1])]
    out = pd.DataFrame(index=p.index)
    out["T"] = (p.TMAX + p.TMIN) / 20.0   # tenths of deg C -> deg C, mean of max and min
    out["P"] = p.PRCP / 10.0              # tenths of mm -> mm
    return out, used


def anomaly_fit(x):
    """Fit behind the anomalies: the straight-line trend, and the smoothed day-of-year mean and SD
    of the detrended values (31-day running means, wrapping round the year), one value per day."""
    doy = np.minimum(np.asarray(x.index.dayofyear), 365)
    yr = np.asarray(x.index.year) + doy / 365
    g = np.asarray(x.notna())
    c = np.polyfit(yr[g], x.values[g], 1)
    trend = np.polyval(c, yr)
    r = x - trend
    m, s = r.groupby(doy).mean(), r.groupby(doy).std()
    k = 31
    smooth = lambda v: np.convolve(np.r_[v[-k:], v, v[:k]], np.ones(k) / k, "same")[k:-k]
    mm, ss = pd.Series(smooth(m.values), index=m.index), pd.Series(smooth(s.values), index=s.index)
    return trend, mm.reindex(doy).values, ss.reindex(doy).values


def anomalies(x):
    """Standardised anomaly: detrend (straight line), remove the smoothed day-of-year mean,
    divide by the smoothed day-of-year SD."""
    trend, day_mean, day_sd = anomaly_fit(x)
    return (x - trend - day_mean) / day_sd


def departures(x):
    """The same anomaly in the variable's own units (K for temperature): not divided by the SD."""
    trend, day_mean, day_sd = anomaly_fit(x)
    return x - trend - day_mean


def anomaly_if(x, value):
    """The standardised anomaly each day would have had if its value had been `value` (0 mm: a dry day)."""
    trend, day_mean, day_sd = anomaly_fit(x)
    return pd.Series((value - trend - day_mean) / day_sd, index=x.index)


# ---------------------------------------------------------------------------------------
# 5. Station test against randomly timed storms
# ---------------------------------------------------------------------------------------
def p_values(vals, obs):
    """p: two-tailed, percentile-based. p_from_zero: share of random sets at least as far from zero."""
    return dict(p=float(min(1.0, 2 * min(np.mean(vals >= obs), np.mean(vals <= obs)))),
                p_from_zero=float(np.mean(np.abs(vals) >= abs(obs))))


def station_test(an, storms, can_start, rng, ndraw, dry=None):
    """Mean anomaly over the storm days with data, against `ndraw` random sets of pseudo-storms with the
    same number and lengths of storms. `can_start` gives, for a DatetimeIndex, a boolean mask of the
    days a pseudo-storm may start on (days of the same season). `dry`, if given, is the anomaly each
    day would have had with no precipitation; a completely dry set of storm days is then scored too."""
    an = an.dropna()
    have = set(an.index)
    evd = [[x for x in e if x in have] for e in storms]
    evd = [e for e in evd if e]
    obs = float(np.mean([an[x] for e in evd for x in e]))
    lens = np.array([len(e) for e in evd])
    pool = an.index[can_start(an.index)]
    # Daily array covering the series plus padding, NaN where there is no anomaly.
    t0 = an.index.min()
    full = pd.date_range(t0, an.index.max() + pd.Timedelta(days=int(lens.max()) + 1), freq="D")
    arr = np.full(len(full), np.nan)
    arr[(an.index - t0).days] = an.values
    pool_pos = np.asarray((pool - t0).days)
    # Same draw order as the original: for each random set, one start per storm in turn.
    idx = rng.integers(0, len(pool_pos), size=ndraw * len(lens)).reshape(ndraw, len(lens))
    starts = pool_pos[idx]
    tot = np.zeros(ndraw)
    cnt = np.zeros(ndraw)
    for j, L in enumerate(lens):
        for d in range(L):
            v = arr[starts[:, j] + d]
            ok = ~np.isnan(v)
            tot[ok] += v[ok]
            cnt[ok] += 1
    with np.errstate(invalid="ignore", divide="ignore"):
        vals = tot / cnt
    lo, hi = np.percentile(vals, [2.5, 97.5])
    res = dict(obs=obs, lo=float(lo), hi=float(hi), **p_values(vals, obs),
               percentile=float(np.mean(vals < obs) * 100), storms=len(evd), days=int(lens.sum()),
               null_median=float(np.median(vals)))
    if dry is not None:
        dry_obs = float(np.mean([dry[x] for e in evd for x in e]))
        pv = p_values(vals, dry_obs)
        res.update(dry_obs=dry_obs, dry_p=pv["p"], dry_p_from_zero=pv["p_from_zero"])
    return res


def stations_drier(P, storms, months, rng, n):
    """Number of stations drier than average over the storm days, and the same count for `n` random sets
    of storms with the same lengths starting on days of the same months (one start per storm in turn)."""
    days = [d for e in storms for d in e]
    real = P.reindex(days).mean()
    pool = P.index[P.index.month.isin(months)]
    lens = [len(e) for e in storms]
    idx = rng.integers(0, len(pool), size=n * len(lens)).reshape(n, len(lens))
    pos = {d: i for i, d in enumerate(P.index)}
    V = P.values
    negs = np.zeros(n, int)
    for i in range(n):
        rows = np.array([pos.get(pool[idx[i, j]] + pd.Timedelta(days=d), -1)
                         for j, L in enumerate(lens) for d in range(L)])
        vv = np.where(rows[:, None] >= 0, V[np.maximum(rows, 0)], np.nan)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", RuntimeWarning)  # a station with no data in a set gives NaN
            negs[i] = int(np.sum(np.nanmean(vv, axis=0) < 0))
    return real, negs


# ---------------------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------------------
def main():
    ap = argparse.ArgumentParser(description=__doc__.split("\n")[0])
    ap.add_argument("--cache", default="station_check_data",
                    help="folder for (or with) the raw downloads (default: ./station_check_data)")
    ap.add_argument("--out", default="station_check_results",
                    help="folder for station-check-results.json, storms.csv and station_tests.csv "
                         "(default: ./station_check_results)")
    a = ap.parse_args()
    os.makedirs(a.out, exist_ok=True)
    pct = lambda v: round(float(v) * 100, 1)
    r3 = lambda v: round(float(v), 3)

    # ---- Dst and storms
    dst, ky, om = load_dst(a.cache)
    dmin = dst.groupby(dst.index.normalize()).min()
    counts = {}
    for level, (lo, hi) in LEVELS.items():
        for season, months in SEASONS.items():
            ev = storms_for(dmin, level, months)
            hrs = int(((dst > lo) & (dst <= hi) & dst.index.month.isin(months)).sum())
            counts[f"{level}_{season}"] = {"storms": len(ev), "hours": hrs}
    storms = storm_table(dst, dmin)
    storms.to_csv(os.path.join(a.out, "storms.csv"), index=False)
    win = SEASONS["winter"]
    in_winter = lambda idx: idx.month.isin(win)
    R4 = storms_for(dmin, "strong", win)
    R3 = storms_for(dmin, "moderate", win)

    # ---- Stations and the 40 winter tests (seed 1)
    rng = np.random.default_rng(SEED_TESTS)
    tests, A, X, DRY, meta = [], {}, {}, {}, []
    for sid_new, sid_old, name in STATIONS:
        x, used = load_station(a.cache, sid_new, sid_old)
        X[name] = x
        DRY[name] = anomaly_if(x["P"], 0.0)
        for v in ["T", "P"]:
            an = anomalies(x[v])
            A[(name, v)] = an
            for lab, ev in [("strong", R4), ("moderate", R3)]:
                dry = DRY[name] if (v == "P" and lab == "strong") else None
                r = station_test(an, ev, in_winter, rng, N_TEST, dry=dry)
                tests.append(dict(station=name, ghcn_id=used, variable=v, level=lab, **r))
        meta.append(dict(station=name, ghcn_id=used, ghcn_id_on_2026_09_30=sid_old,
                         first_T=str(x["T"].first_valid_index().date()), last_T=str(x["T"].last_valid_index().date()),
                         first_P=str(x["P"].first_valid_index().date()), last_P=str(x["P"].last_valid_index().date())))
    TT = pd.DataFrame(tests)
    out_tests = TT.drop(columns=["dry_obs", "dry_p", "dry_p_from_zero"])
    for c in ["obs", "lo", "hi", "null_median"]:
        out_tests[c + "_pct_sd"] = (out_tests[c] * 100).round(1)
    out_tests.drop(columns=["obs", "lo", "hi", "null_median"]).to_csv(os.path.join(a.out, "station_tests.csv"), index=False)

    T = pd.DataFrame({k[0]: v for k, v in A.items() if k[1] == "T"})
    P = pd.DataFrame({k[0]: v for k, v in A.items() if k[1] == "P"})
    TK = pd.DataFrame({name: departures(X[name]["T"]) for name in X})
    PMM = pd.DataFrame({name: X[name]["P"] for name in X})

    # ---- Storm by storm (10-station average of the day's anomaly, averaged over storm days)
    per_storm = []
    for e in R4:
        h = dst[dst.index.normalize().isin(e)]
        per_storm.append(dict(first_day=str(e[0].date()), last_day=str(e[-1].date()), min_dst=int(h.min()),
                              hours_le_minus200=int((h <= -200).sum()),
                              T_pct_sd=round(float(T.reindex(e).mean(axis=1).mean()) * 100),
                              P_pct_sd=round(float(P.reindex(e).mean(axis=1).mean()) * 100),
                              stations_with_T=int(T.reindex(e).notna().any().sum())))

    def day_detail(days):
        """Each station's temperature (K and % of SD) and precipitation (mm and % of SD) on the given days."""
        out = []
        for d in days:
            rows = []
            for name in X:
                g = lambda df: df[name].get(d, np.nan)
                rows.append(dict(station=name,
                                 T_K=None if pd.isna(g(TK)) else round(float(g(TK)), 1),
                                 T_pct_sd=None if pd.isna(g(T)) else pct(g(T)),
                                 P_mm=None if pd.isna(g(PMM)) else round(float(g(PMM)), 1),
                                 P_pct_sd=None if pd.isna(g(P)) else pct(g(P))))
            out.append(dict(day=str(d.date()), min_dst=int(dmin[d]), stations=rows))
        return out

    # ---- Coherence and precipitation sign (seed 11; draw order as the original)
    rng2 = np.random.default_rng(SEED_FIELD)
    pool = T.index[T.index.month.isin(win)]
    lens = [len(e) for e in R4]
    idx = rng2.integers(0, len(pool), size=N_COHERENCE * len(lens)).reshape(N_COHERENCE, len(lens))
    Tv, Pv = T.values, P.values
    tpos = {d: i for i, d in enumerate(T.index)}
    same, negP = np.zeros(N_COHERENCE, int), np.zeros(N_COHERENCE, int)
    for i in range(N_COHERENCE):
        rows = []
        for j, L in enumerate(lens):
            st = pool[idx[i, j]]
            rows += [tpos.get(st + pd.Timedelta(days=d), -1) for d in range(L)]
        rows = np.array(rows)
        tt = np.where(rows[:, None] >= 0, Tv[np.maximum(rows, 0)], np.nan)
        pp = np.where(rows[:, None] >= 0, Pv[np.maximum(rows, 0)], np.nan)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", RuntimeWarning)  # a station with no data in a set gives NaN
            tm, pm = np.nanmean(tt, axis=0), np.nanmean(pp, axis=0)
        same[i] = max(np.sum(tm < 0), np.sum(tm > 0))
        negP[i] = np.sum(pm < 0)
    r4days = [d for e in R4 for d in e]
    obsT, obsP = T.reindex(r4days).mean(), P.reindex(r4days).mean()
    obs_same = int(max((obsT < 0).sum(), (obsT > 0).sum()))
    obs_negP = int((obsP < 0).sum())
    skew = {}
    for c in P.columns:
        v = P[c][P.index.month.isin(win)].dropna().values
        draws = v[rng2.choice(len(v), (N_SKEW, 5))].mean(1)
        skew[c] = round(float(np.mean(draws < 0)), 3)

    # ---- Lag composite around the 77 moderate winter storm onsets (seed 11, continuing)
    Tm = T.mean(axis=1)
    tm_full = Tm.reindex(pd.date_range(Tm.index.min() - pd.Timedelta(days=10), Tm.index.max() + pd.Timedelta(days=10)))
    base = tm_full.index[0]
    tmv = tm_full.values
    onsets = [e[0] for e in R3]
    lags = list(range(-10, 11))
    opos = np.array([(o - base).days for o in onsets])
    lag_obs = [float(np.nanmean(tmv[opos + L])) for L in lags]
    ridx = rng2.integers(0, len(pool), size=N_LAG * len(onsets)).reshape(N_LAG, len(onsets))
    rpos = np.asarray((pool - base).days)[ridx]
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", RuntimeWarning)
        nul = np.array([[np.nanmean(tmv[rpos[i] + L]) for L in lags] for i in range(N_LAG)])
    lag_lo, lag_hi = np.percentile(nul, 2.5, axis=0), np.percentile(nul, 97.5, axis=0)

    # ---- Strong storms and precipitation in every season (seeds 3 and 13)
    rng3 = np.random.default_rng(SEED_SEASON_TESTS)
    rng4 = np.random.default_rng(SEED_SEASON_SIGN)
    by_season = {}
    for season, months in SEASONS.items():
        ev = storms_for(dmin, "strong", months)
        can_start = lambda ix, m=months: ix.month.isin(m)
        rows = []
        for name in P.columns:
            r = station_test(P[name], ev, can_start, rng3, N_TEST, dry=DRY[name])
            rows.append(dict(station=name, real=pct(r["obs"]), random_lo=pct(r["lo"]), random_hi=pct(r["hi"]),
                             p=r3(r["p"]), p_from_zero=r3(r["p_from_zero"]), storms=r["storms"],
                             completely_dry=pct(r["dry_obs"]), completely_dry_p=r3(r["dry_p"]),
                             completely_dry_p_from_zero=r3(r["dry_p_from_zero"])))
        real, negs = stations_drier(P, ev, months, rng4, N_COHERENCE)
        n_dry = int((real < 0).sum())
        by_season[season] = {
            "strong_storms": len(ev), "storm_days": len([d for e in ev for d in e]),
            "stations_drier": n_dry, "wetter_stations": [k for k, v in real.items() if v > 0],
            "share_random_sets_at_least_as_many_drier": r3(np.mean(negs >= n_dry)),
            "median_random_stations_drier": int(np.median(negs)),
            "tests_p_below_0.05": [dict(station=t["station"], real=t["real"], p=t["p"]) for t in rows if t["p"] < 0.05],
            "tests_p_from_zero_below_0.05": [dict(station=t["station"], real=t["real"], p_from_zero=t["p_from_zero"])
                                             for t in rows if t["p_from_zero"] < 0.05],
            "stations_where_a_completely_dry_set_gives_p_below_0.05": sum(t["completely_dry_p"] < 0.05 for t in rows),
            "station_tests": rows}

    # ---- Winter as +-30 days around 21 December, the 2024 poster's window (separate generator, seed 1)
    all_days = pd.date_range(f"{YEARS[0] - 1}-01-01", f"{YEARS[1] + 1}-12-31", freq="D")
    window_days = all_days[[near_dec21(d) for d in all_days]]
    in_window = lambda ix: ix.isin(window_days)
    R4w = storms_for(dmin, "strong", keep=near_dec21)
    R3w = storms_for(dmin, "moderate", keep=near_dec21)
    rng5 = np.random.default_rng(SEED_WINDOW)
    wrows = []
    for sid_new, sid_old, name in STATIONS:
        for v in ["T", "P"]:
            for lab, ev in [("strong", R4w), ("moderate", R3w)]:
                r = station_test(A[(name, v)], ev, in_window, rng5, N_TEST)
                wrows.append(dict(station=name, variable=v, level=lab, real=pct(r["obs"]), random_lo=pct(r["lo"]),
                                  random_hi=pct(r["hi"]), p=r3(r["p"]), p_from_zero=r3(r["p_from_zero"]),
                                  storms=r["storms"]))
    wdays = [d for e in R4w for d in e]
    wT, wP = T.reindex(wdays).mean(), P.reindex(wdays).mean()
    hour_in_window = pd.Series(dst.index.normalize(), index=dst.index).map(near_dec21).values
    window = {
        "winter": "+-30 days of 21 December",
        "strong_storms": [dict(first_day=str(e[0].date()), last_day=str(e[-1].date()), min_dst=int(dmin[e].min()),
                               T_pct_sd=round(float(T.reindex(e).mean(axis=1).mean()) * 100),
                               P_pct_sd=round(float(P.reindex(e).mean(axis=1).mean()) * 100)) for e in R4w],
        "strong_storm_hours_le_minus200": int(((dst <= -200) & hour_in_window).sum()),
        "moderate_storms": len(R3w),
        "n_tests": len(wrows),
        "n_p_below_0.05": sum(t["p"] < 0.05 for t in wrows),
        "tests_p_below_0.05": [dict(station=t["station"], variable=t["variable"], level=t["level"], real=t["real"],
                                    p=t["p"]) for t in wrows if t["p"] < 0.05],
        "n_p_from_zero_below_0.05": sum(t["p_from_zero"] < 0.05 for t in wrows),
        "stations_colder": int((wT < 0).sum()), "stations_drier": int((wP < 0).sum()),
        "station_tests": wrows,
        "strong_storm_days": day_detail(wdays),
    }

    # ---- Other sensitivity checks
    dst5723 = pd.concat([ky, om[om.index.year >= 1963]]).sort_index()
    dst5723 = dst5723[(dst5723.index.year >= 1957) & (dst5723.index.year <= 2023)]
    dmin5723 = dst5723.groupby(dst5723.index.normalize()).min()
    desc = lambda ev, dm: [f"{e[0].date()} ({int(dm[e].min())} nT)" for e in ev]
    djf = dst.index.month.isin(win)
    sens = {
        "strong_winter_storms_DJF_1958_2024": desc(R4, dmin),
        "strong_winter_storms_DJF_1957_2023": desc(storms_for(dmin5723, "strong", win), dmin5723),
        "strong_winter_storms_pm30d_of_Dec21_1958_2024": desc(R4w, dmin),
        "hours_le_minus200_DJF_1957_2023": int(((dst5723 <= -200) & dst5723.index.month.isin(win)).sum()),
        "hours_le_minus200_pm30d_of_Dec21_1958_2024": int(((dst <= -200) & hour_in_window).sum()),
        "winter_hours_total_DJF_1958_2024": int(djf.sum()),
        "winter_hours_R1_bin_0_to_minus50_DJF_1958_2024": int(((dst > -50) & (dst <= 0) & djf).sum()),
        "winter_hours_Dst_below_0_DJF_1958_2024": int(((dst < 0) & djf).sum()),
        "winter_hours_R1_bin_0_to_minus50_pm30d_of_Dec21": int(((dst > -50) & (dst <= 0) & hour_in_window).sum()),
        "moderate_winter_storms_within_2_days_of_a_strong_storm_day": [
            str(e[0].date()) for e in R3 if any(abs((d - r).days) <= 2 for d in e for e2 in R4 for r in e2)],
    }
    # storms defined by peak: runs of days <= -100 nT, classed by their lowest Dst
    by_peak = {}
    for season, months in SEASONS.items():
        ev = group_days([x for x in dmin.index[dmin <= -100] if x.month in months])
        by_peak[season] = {"moderate": sum(1 for e in ev if -200 < dmin[e].min() <= -100),
                           "strong": sum(1 for e in ev if dmin[e].min() <= -200)}
    sens["storms_by_peak_runs_of_days_le_minus100"] = by_peak
    # hour-based grouping: hours in the level merged when less than 48 h apart
    by_hours = {}
    for level, (lo, hi) in LEVELS.items():
        for season, months in SEASONS.items():
            h = dst[(dst > lo) & (dst <= hi) & dst.index.month.isin(months)].index
            n = 1 + int(np.sum(np.diff(h.values).astype("timedelta64[h]").astype(int) > 48)) if len(h) else 0
            by_hours[f"{level}_{season}"] = n
    sens["storms_by_hour_grouping_gap_gt_48h"] = by_hours

    # ---- Assemble
    strongT = TT[(TT.level == "strong") & (TT.variable == "T")]
    strongP = TT[(TT.level == "strong") & (TT.variable == "P")]
    station_rows = lambda df: [dict(station=r.station, real=pct(r.obs), random_lo=pct(r.lo), random_hi=pct(r.hi),
                                    p=r3(r.p), p_from_zero=r3(r.p_from_zero), percentile=round(r.percentile, 1),
                                    storms=int(r.storms), days=int(r.days),
                                    real_inside_random_range=bool(r.lo <= r.obs <= r.hi)) for r in df.itertuples()]
    ghcn_version = None
    try:
        p = fetch(GHCN_VERSION_URL, os.path.join(a.cache, "ghcn", "ghcnd-version.txt"))
        ghcn_version = open(p).read().split("\n")[0].strip()
    except Exception:
        pass
    res = {
        "generated_utc": dt.datetime.now(dt.timezone.utc).strftime("%Y-%m-%d %H:%M"),
        "data": {"dst": "WDC Kyoto final Dst 1958-1962; NASA SPDF OMNI2 hourly (word 41 = Dst) 1963-2024",
                 "omni_file_last_year": int(om.index.year.max()), "ghcn_daily_version": ghcn_version,
                 "missing_dst_hours_1958_2024": int(dst.isna().sum()), "stations": meta},
        "settings": {"years": YEARS, "winter_months": win, "levels_daily_min_dst": LEVELS,
                     "merge_gap_days": 2, "n_random_per_station_test": N_TEST, "n_random_coherence": N_COHERENCE,
                     "n_random_skew_per_station": N_SKEW, "n_random_lag": N_LAG,
                     "seeds": {"winter_tests": SEED_TESTS, "coherence_skew_lag": SEED_FIELD,
                               "seasons_tests": SEED_SEASON_TESTS, "seasons_drier_counts": SEED_SEASON_SIGN,
                               "solstice_window_tests": SEED_WINDOW},
                     "p": "two-tailed, percentile-based: min(1, 2 x min(share of random sets >= real, share <= real))",
                     "p_from_zero": "share of random sets whose mean is at least as far from zero as the real one"},
        "storm_counts": counts,
        "strong_winter_storms": per_storm,
        "strong_winter_storm_days": day_detail(r4days),
        "chart_temperature_strong_winter": station_rows(strongT),
        "chart_precipitation_strong_winter": station_rows(strongP),
        "tests_40": {"n": len(TT), "n_p_below_0.05": int((TT.p < 0.05).sum()),
                     "below_0.05": TT.loc[TT.p < 0.05, ["station", "variable", "level", "p"]].to_dict("records"),
                     "n_p_from_zero_below_0.05": int((TT.p_from_zero < 0.05).sum()),
                     "below_0.05_p_from_zero": TT.loc[TT.p_from_zero < 0.05, ["station", "variable", "level", "p_from_zero"]].to_dict("records"),
                     "single_station_random_range_temperature_strong": [pct(strongT.lo.min()), pct(strongT.lo.max()),
                                                                         pct(strongT.hi.min()), pct(strongT.hi.max())],
                     "single_station_random_range_precip_strong": [pct(strongP.lo.min()), pct(strongP.lo.max()),
                                                                    pct(strongP.hi.min()), pct(strongP.hi.max())],
                     "tests": [dict(station=r.station, variable=r.variable, level=r.level, real=pct(r.obs),
                                    random_lo=pct(r.lo), random_hi=pct(r.hi), p=r3(r.p), p_from_zero=r3(r.p_from_zero),
                                    storms=int(r.storms), days=int(r.days)) for r in TT.itertuples()]},
        "coherence": {"real_same_sign_stations": obs_same,
                      "real_signs": {k: ("+" if v > 0 else "-") for k, v in obsT.items()},
                      "share_random_sets_9_or_10_same_sign": round(float(np.mean(same >= 9)), 3)},
        "precipitation_sign": {"real_stations_drier": obs_negP,
                               "real_wetter_stations": [k for k, v in obsP.items() if v > 0],
                               "share_random_sets_at_least_as_many_drier": round(float(np.mean(negP >= obs_negP)), 3),
                               "median_random_stations_drier": int(np.median(negP)),
                               "share_5_random_winter_days_drier_by_station": skew,
                               "range_5_random_winter_days_drier": [min(skew.values()), max(skew.values())]},
        "precipitation_power_strong_winter": {
            "share_random_sets_all_10_stations_drier": r3(np.mean(negP >= 10)),
            "share_random_sets_at_least_9_stations_drier": r3(np.mean(negP >= 9)),
            "completely_dry_storm_days": [dict(station=r.station, real=pct(r.obs), completely_dry=pct(r.dry_obs),
                                               p=r3(r.dry_p), p_from_zero=r3(r.dry_p_from_zero))
                                          for r in strongP.itertuples()],
            "stations_where_a_completely_dry_set_gives_p_below_0.05": int((strongP.dry_p < 0.05).sum()),
            "stations_where_a_completely_dry_set_gives_p_from_zero_below_0.05": int((strongP.dry_p_from_zero < 0.05).sum())},
        "lag_composite_moderate_winter": {"n_onsets": len(onsets), "lags_days": lags,
                                          "real": [round(v * 100, 1) for v in lag_obs],
                                          "random_lo": [round(v * 100, 1) for v in lag_lo],
                                          "random_hi": [round(v * 100, 1) for v in lag_hi],
                                          "lags_outside_random_range": [L for L, v, l, h in zip(lags, lag_obs, lag_lo, lag_hi)
                                                                        if v > h or v < l]},
        "precipitation_strong_storms_by_season": by_season,
        "winter_as_pm30d_of_Dec21": window,
        "sensitivity": sens,
    }
    json.dump(res, open(os.path.join(a.out, "station-check-results.json"), "w", encoding="utf-8"), indent=1,
              ensure_ascii=False, default=str)
    print(json.dumps({k: res[k] for k in ["storm_counts", "strong_winter_storms", "coherence", "precipitation_sign",
                                          "precipitation_power_strong_winter"]}, indent=1, ensure_ascii=False, default=str))
    print("Winter tests with p below 0.05:", res["tests_40"]["below_0.05"])
    for s, v in by_season.items():
        print(f"{s}: {v['strong_storms']} strong storms, {v['stations_drier']} of 10 stations drier "
              f"({v['share_random_sets_at_least_as_many_drier']:.0%} of random sets as dry or drier)")


if __name__ == "__main__":
    main()