#!/usr/bin/env python3 """Every number of ours on factchecked.science/storm-geophysical-effects/, recomputed from public data. python3 -I storm_numbers.py [data-folder] It checks the storm of 20-23 March 2026 against GFZ's Kp and Ap and WDC Kyoto's Dst and SYM-H, works out the dipole (geomagnetic) latitude of the Mikhnevo observatory from IGRF-14, and screens the weather reports of Moscow Domodedovo airport, 51 km away, against three fair-weather criteria for atmospheric electricity. Python 3.10 or later, standard library only. Inputs are read from the data folder (default: ./storm_numbers_data) and downloaded there first if missing; the script prints whether each matches, by SHA-256, the copy we used on October 31, 2026. WDC Kyoto's provisional values can be revised, so a newer download may differ slightly. Outputs, written next to the data folder: storm-numbers.json (every number below) and storm-indices.csv (hourly Dst and SYM-H with each hour's Kp, 19-24 March 2026). Terms of the inputs, which carry over to the outputs: GFZ Kp and Ap, CC BY 4.0 (doi:10.5880/Kp.0001); WDC Kyoto Dst (doi:10.17593/14515-74000) and ASY/SYM (doi:10.14989/267216), provisional values, for scientific use, with no commercial use, acknowledging the observatories that supply the data and the WDC for Geomagnetism, Kyoto, and citing the data DOIs; IGRF-14 (IAGA) and the IRIS station service, public; the Iowa Environmental Mesonet's METAR archive, public domain. """ import csv import hashlib import json import math import re import sys import urllib.request from datetime import datetime, timedelta, timezone from pathlib import Path INPUTS = { "kp.json": ("https://kp.gfz.de/app/json/?start=2026-03-17T00:00:00Z&end=2026-03-26T23:59:59Z&index=Kp", "57535c42a06d20ddcb83fa90aeb3c958053673a8236916de41be17cd0a545f94"), "ap_daily.json": ("https://kp.gfz.de/app/json/?start=2026-03-17T00:00:00Z&end=2026-03-26T23:59:59Z&index=Ap", "b1c967a5b2972682c926ecc7c43fc8b8d31a9b99c49258c2cf85581d5db3bba0"), "dst2603.for.request": ("https://wdc.kugi.kyoto-u.ac.jp/dst_provisional/202603/dst2603.for.request", "a5ad67607697747c1f500f93f3660a700b81201732d0dc4fc85d600d7dca0fa8"), "symh.csv": ("https://wdc.kugi.kyoto-u.ac.jp/hapi/data?id=min_asysym&time.min=2026-03-19T00:00:00Z" "&time.max=2026-03-25T00:00:00Z¶meters=symHValue,versionCode&format=csv", "e9f9db54e777ab5de88940d5aca1ead88b0f216faa50b1e77a3df66f3d5af14e"), "igrf14coeffs.txt": ("https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt", "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c"), "mhv.txt": ("https://service.iris.edu/fdsnws/station/1/query?sta=MHV&level=station&format=text", "a1d7174e19469e411de253ccfca8ac60ba90e51fe63355f3acf02029bed325d6"), "uudd.csv": ("https://mesonet.agron.iastate.edu/cgi-bin/request/asos.py?station=UUDD&data=tmpc&data=dwpc" "&data=drct&data=sknt&data=gust&data=mslp&data=alti&data=wxcodes&data=skyc1&data=skyl1&data=skyc2" "&data=skyl2&data=vsby&data=metar&year1=2026&month1=3&day1=19&year2=2026&month2=3&day2=25" "&tz=Etc%2FUTC&format=onlycomma&latlon=yes&missing=M&trace=T&direct=no&report_type=3&report_type=4", "e51f5c2328899e4a2a213cb804d7626ddbef3e83fcdc4b7b72b9f96a6f385730"), } UTC = timezone.utc # The two stages (our definition, from GFZ's Kp): the two spans, around the two Kp 7 maxima and the two Dst minima, # in which Kp stood at 5- or more. Kp touched 5- and 5o again at 12-18 UT on 23 March, in the recovery, with no new # Dst minimum; the script lists that span too, but it is not a stage. STAGES = [("stage 1", datetime(2026, 3, 20, 18, tzinfo=UTC), datetime(2026, 3, 21, 9, tzinfo=UTC)), ("stage 2", datetime(2026, 3, 22, 3, tzinfo=UTC), datetime(2026, 3, 23, 3, tzinfo=UTC))] def fetch(folder): folder.mkdir(parents=True, exist_ok=True) for name, (url, sha) in INPUTS.items(): path = folder / name if not path.exists(): print(f"downloading {name} ...") req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0 (storm_numbers.py)"}) with urllib.request.urlopen(req, timeout=120) as r: path.write_bytes(r.read()) got = hashlib.sha256(path.read_bytes()).hexdigest() print(f" {name}: {'same file as ours' if got == sha else 'differs from ours (sha256 ' + got[:16] + '...)'}") def kp_symbol(v): thirds = round(v * 3) whole, rest = divmod(thirds, 3) return f"{whole}o" if rest == 0 else (f"{whole}+" if rest == 1 else f"{whole + 1}-") def main(): folder = Path(sys.argv[1] if len(sys.argv) > 1 else "storm_numbers_data").resolve() fetch(folder) out = {"terms": "Kp and Ap: GFZ, CC BY 4.0 (doi:10.5880/Kp.0001). Dst (doi:10.17593/14515-74000) and SYM-H " "(doi:10.14989/267216): WDC for Geomagnetism, Kyoto, provisional values, for scientific use and " "no commercial use, acknowledging the observatories that supply the data and the WDC. " "METARs: Iowa Environmental Mesonet, public domain."} # ---- Kp and Ap (GFZ) ---- kp = json.loads((folder / "kp.json").read_text(encoding="utf-8")) rows = [(datetime.strptime(t, "%Y-%m-%dT%H:%M:%SZ").replace(tzinfo=UTC), v, s) for t, v, s in zip(kp["datetime"], kp["Kp"], kp["status"])] storm = [(t, v) for t, v, _ in rows if datetime(2026, 3, 20, tzinfo=UTC) <= t < datetime(2026, 3, 24, tzinfo=UTC)] kmax = max(v for _, v in storm) print("\nKp, GFZ (definitive), 3-hour intervals, 20-23 March 2026:") for day in range(20, 24): print(f" {day} March: " + " ".join(kp_symbol(v) for t, v in storm if t.day == day)) out["kp_status"] = sorted({s for _, _, s in rows}) out["kp_max"] = {"value": kmax, "symbol": kp_symbol(kmax), "intervals_utc": [f"{t:%Y-%m-%d %H}-{t + timedelta(hours=3):%H}" for t, v in storm if v == kmax]} out["kp_intervals_5minus_or_more"] = [sum(1 for _, v in storm if v >= 14 / 3 - 1e-6), len(storm)] print(f" maximum {kp_symbol(kmax)} in {out['kp_max']['intervals_utc']}; Kp 5- or more in " f"{out['kp_intervals_5minus_or_more'][0]} of {len(storm)} intervals") # Every run of consecutive intervals at 5- or more, and whether it holds a Kp maximum (two are the stages) spans, cur = [], None for t, v in storm: if v >= 14 / 3 - 1e-6: if cur is None: cur = [t, t + timedelta(hours=3), v] else: cur[1], cur[2] = t + timedelta(hours=3), max(cur[2], v) elif cur is not None: spans.append(cur) cur = None if cur is not None: spans.append(cur) out["kp_spans_5minus_or_more"] = [{"from_utc": f"{a:%Y-%m-%d %H:%M}", "to_utc": f"{b:%Y-%m-%d %H:%M}", "highest": kp_symbol(m), "stage": m == kmax} for a, b, m in spans] print(" spans at 5- or more:", json.dumps(out["kp_spans_5minus_or_more"])) ap = json.loads((folder / "ap_daily.json").read_text(encoding="utf-8")) out["ap_daily"] = {d[:10]: v for d, v in zip(ap["datetime"], ap["Ap"])} print(" daily Ap:", out["ap_daily"]) # ---- Dst (WDC Kyoto, provisional, hourly). The file's 24 values a day are the means of the hours 00-01 ... 23-24 # UT; each is labelled here by the hour's start ---- dst = {} for line in (folder / "dst2603.for.request").read_text(encoding="ascii").splitlines(): if line.startswith("DST2603"): day = int(line[8:10]) for h in range(24): dst[datetime(2026, 3, day, h, tzinfo=UTC)] = int(line[20 + 4 * h:24 + 4 * h]) def dmin(t0, t1): sel = {t: v for t, v in dst.items() if t0 <= t < t1} t = min(sel, key=sel.get) return t, sel[t] d1 = dmin(datetime(2026, 3, 20, 12, tzinfo=UTC), datetime(2026, 3, 22, tzinfo=UTC)) d2 = dmin(datetime(2026, 3, 22, tzinfo=UTC), datetime(2026, 3, 24, tzinfo=UTC)) between = {t: v for t, v in dst.items() if d1[0] < t < d2[0]} tb = max(between, key=between.get) day20 = {t: v for t, v in dst.items() if t.day == 20} t20 = max(day20, key=day20.get) n50 = sum(1 for t, v in dst.items() if 20 <= t.day <= 24 and v <= -50) out["dst"] = {"stage1_min": [f"{d1[0]:%Y-%m-%d %H}-{d1[0] + timedelta(hours=1):%H} UT", d1[1]], "stage2_min": [f"{d2[0]:%Y-%m-%d %H}-{d2[0] + timedelta(hours=1):%H} UT", d2[1]], "recovery_between": [f"{tb:%Y-%m-%d %H}-{tb + timedelta(hours=1):%H} UT", between[tb]], "max_20_march": [f"{t20:%Y-%m-%d %H}-{t20 + timedelta(hours=1):%H} UT", day20[t20]], "hours_at_or_below_minus50_20_24_march": n50, "hours_at_or_below_minus100": [f"{t:%Y-%m-%d %H} UT" for t in sorted(dst) if dst[t] <= -100]} print("\nDst, WDC Kyoto (provisional):", json.dumps(out["dst"])) # ---- SYM-H (WDC Kyoto, provisional, 1 minute) ---- sym = [] for r in csv.reader((folder / "symh.csv").read_text(encoding="ascii").splitlines()): v = int(r[1]) if v != 99999: sym.append((datetime.strptime(r[0], "%Y-%m-%dT%H:%M:%SZ").replace(tzinfo=UTC), v, r[2])) def smin(t0, t1): return min(((t, v) for t, v, _ in sym if t0 <= t < t1), key=lambda x: x[1]) s1 = smin(datetime(2026, 3, 20, 12, tzinfo=UTC), datetime(2026, 3, 22, tzinfo=UTC)) s2 = smin(datetime(2026, 3, 22, tzinfo=UTC), datetime(2026, 3, 24, tzinfo=UTC)) smax = max(((t, v) for t, v, _ in sym if t.day == 20), key=lambda x: x[1]) # The sudden rise on 20 March: the largest 6-minute rise before 12 UT rises = [(sym[i][1] - sym[i - 6][1], sym[i - 6][0], sym[i - 6][1], sym[i][1]) for i in range(6, len(sym)) if sym[i][0].day == 20 and sym[i][0].hour < 12] rise = max(rises) out["symh"] = {"version_codes": sorted({c for _, _, c in sym}), "stage1_min": [f"{s1[0]:%Y-%m-%d %H:%M} UT", s1[1]], "stage2_min": [f"{s2[0]:%Y-%m-%d %H:%M} UT", s2[1]], "max_20_march": [f"{smax[0]:%Y-%m-%d %H:%M} UT", smax[1]], "sudden_rise_20_march": {"from": f"{rise[1]:%H:%M} UT", "minutes": 6, "nT": rise[0], "values": [rise[2], rise[3]]}} print("SYM-H, WDC Kyoto (provisional):", json.dumps(out["symh"])) # ---- Mikhnevo's dipole latitude (IGRF-14, 20 March 2026) ---- epoch = 2026 + (31 + 28 + 19) / 365.0 g = {} for line in (folder / "igrf14coeffs.txt").read_text(encoding="ascii").splitlines(): p = line.split() if len(p) > 4 and p[0] in ("g", "h") and p[1] == "1": g[(p[0], int(p[2]))] = float(p[-2]) + (epoch - 2025.0) * float(p[-1]) g10, g11, h11 = g[("g", 0)], g[("g", 1)], g[("h", 1)] b0 = math.sqrt(g10 ** 2 + g11 ** 2 + h11 ** 2) plat = 90 - math.degrees(math.acos(-g10 / b0)) plon = math.degrees(math.atan2(-h11, -g11)) for line in (folder / "mhv.txt").read_text(encoding="utf-8").splitlines(): if line.startswith("GE|MHV|"): mlat, mlon = float(line.split("|")[2]), float(line.split("|")[3]) s = (math.sin(math.radians(mlat)) * math.sin(math.radians(plat)) + math.cos(math.radians(mlat)) * math.cos(math.radians(plat)) * math.cos(math.radians(mlon - plon))) out["mikhnevo"] = {"lat": mlat, "lon": mlon, "dipole_pole": [round(plat, 2), round(plon, 2)], "dipole_latitude": round(math.degrees(math.asin(s)), 1)} print("\nMikhnevo:", json.dumps(out["mikhnevo"])) # ---- Weather at Moscow Domodedovo (UUDD) during the two stages ---- # Fair-weather criteria after Harrison & Nicoll (2018), as paraphrased on the page: FW1 no hydrometeors (rain, # drizzle, snow, fog or mist), aerosol or haze, as shown by a visual range above 2 km; FW2 no cumuliform cloud # and no broken or overcast cloud below 1500 m; FW3 wind of 1-8 m/s. METAR weather codes, visibility, cloud and # wind stand in for the criteria. The wind is the METAR's own value in m/s (IEM's knots, converted back, can # turn a reported 1 m/s into 0.998). wx = list(csv.DictReader((folder / "uudd.csv").read_text(encoding="utf-8").splitlines())) alat, alon = float(wx[0]["lat"]), float(wx[0]["lon"]) dist = 2 * 6371.0 * math.asin(math.sqrt(math.sin(math.radians(alat - mlat) / 2) ** 2 + math.cos(math.radians(alat)) * math.cos(math.radians(mlat)) * math.sin(math.radians(alon - mlon) / 2) ** 2)) hydro = ("RA", "SN", "DZ", "FG", "BR", "GR", "GS", "PL", "SG", "TS", "SH", "FZ", "HZ", "FU", "DU", "SA", "UP") def num(x): try: return float(x) except ValueError: return None def metar_wind(r): m = re.search(r"\b(?:\d{3}|VRB)(\d{2,3})(?:G\d{2,3})?MPS\b", r["metar"]) if m: return float(m.group(1)) k = num(r["sknt"]) return None if k is None else round(k * 0.514444, 1) out["airport"] = {"station": "UUDD", "distance_km": round(dist), "stages": {}} for name, t0, t1 in STAGES: sel = [r for r in wx if f"{t0:%Y-%m-%d %H:%M}" <= r["valid"] < f"{t1:%Y-%m-%d %H:%M}"] res = [] for r in sel: vis = num(r["vsby"]) * 1.609 if num(r["vsby"]) is not None else None codes = "" if r["wxcodes"] in ("M", "") else r["wxcodes"] fw1 = vis is not None and vis > 2 and not any(c in codes for c in hydro) low = any(c in ("BKN", "OVC") and num(h) is not None and num(h) * 0.3048 < 1500 for c, h in ((r["skyc1"], r["skyl1"]), (r["skyc2"], r["skyl2"]))) fw2 = not low and "CB" not in r["metar"] and "TCU" not in r["metar"] wind = metar_wind(r) fw3 = wind is not None and 1 <= wind <= 8 res.append((fw1, fw2, fw3, vis, wind, "BR" in codes)) out["airport"]["stages"][name] = { "window_utc": [f"{t0:%Y-%m-%d %H:%M}", f"{t1:%Y-%m-%d %H:%M}"], "reports": len(sel), "fw1": sum(x[0] for x in res), "fw2": sum(x[1] for x in res), "fw3": sum(x[2] for x in res), "all_three": sum(x[0] and x[1] and x[2] for x in res), "mist_reports": sum(x[5] for x in res), "visibility_km": [round(min(x[3] for x in res), 1), round(max(x[3] for x in res), 1)], "wind_m_s": [round(min(x[4] for x in res), 1), round(max(x[4] for x in res), 1)]} print(f"\nWeather at UUDD, {round(dist)} km from Mikhnevo:", json.dumps(out["airport"]["stages"], indent=1)) # ---- Outputs ---- base = folder.parent (base / "storm-numbers.json").write_text(json.dumps(out, indent=1) + "\n", encoding="utf-8", newline="\n") kp_at = {t: v for t, v, _ in rows} with open(base / "storm-indices.csv", "w", encoding="utf-8", newline="\n") as f: w = csv.writer(f, lineterminator="\n") w.writerow(["hour_start_utc", "dst_nT_provisional", "symh_min_nT_provisional", "kp", "kp_symbol"]) for h in range(6 * 24): t = datetime(2026, 3, 19, tzinfo=UTC) + timedelta(hours=h) k = kp_at[t.replace(hour=t.hour // 3 * 3)] sm = min(v for tt, v, _ in sym if t <= tt < t + timedelta(hours=1)) w.writerow([f"{t:%Y-%m-%dT%H:%MZ}", dst[t], sm, f"{k:.3f}", kp_symbol(k)]) print(f"\nwrote {base / 'storm-numbers.json'} and {base / 'storm-indices.csv'}") if __name__ == "__main__": main()