"""Re-run the ozone-hole numbers on the page "S0 News, October 24, 2026" (factchecked.science). What it prints: 1. the days of 2026 on which the Antarctic ozone hole's area beat or tied the 1979-2025 record for the date, in CAMS's analysis (daily means and six-hourly values) and in NASA Ozone Watch's daily file; 2. how 2026 ranks for each date from 10 to 24 October (CAMS) and how often it was above the 95th percentile of 1979-2023 (CAMS) or the 90th percentile of NASA's climatology; 3. the days on which the 50 hPa minimum temperature was at or below its record low for the date, in both files; 4. the correlation, 1990-2026, between year-to-year changes in the hole's mean area (7 Sep-13 Oct, CAMS) and in the mean 50 hPa minimum temperature (1 Sep-13 Oct, CAMS) or the mean Ap index (1 Aug-13 Oct, GFZ); 5. how much of 2026's change in area the temperature change accounts for (a least-squares line through the 36 changes), whether what it leaves over follows Ap, how 2026's mean Ap ranks, and the stormy days (daily Ap of 30 or more) in its 1 Aug-13 Oct window. Download the five public files first (any folder) and pass the folder as the only argument: CAMS (Copernicus Atmosphere Monitoring Service), https://sites.ecmwf.int/data/cams/ozone_monitoring/data/ cams_ozone_monitoring_sh_ozone_area.csv cams_ozone_monitoring_sh_50hPa_temperature_minimum.csv NASA Ozone Watch, https://ozonewatch.gsfc.nasa.gov/meteorology/figures/ ozone/to3areas_2026_toms+omi+omps.txt merra2/temperature/tmins_50_2026_merra2.txt GFZ Potsdam (CC BY 4.0), https://kp.gfz.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt Run: python3 -I ozone_checks.py FOLDER Needs only Python 3. The page used the files as served on October 26, 2026 (UTC): CAMS to 24 October 18:00 UTC, NASA to 23 October; both agencies revise recent values, so later downloads can differ slightly. """ import csv import datetime as dt import os import re import statistics as st import sys from collections import defaultdict MISSING = -9999.0 def cams(path): six = {} with open(path, newline="") as f: r = csv.reader(f) next(r) for row in r: if len(row) >= 2 and row[0]: try: six[dt.datetime.strptime(row[0], "%Y-%m-%dT%H:%M")] = float(row[1]) except ValueError: pass days = defaultdict(list) for t, v in six.items(): days[t.date()].append(v) return six, {d: sum(v) / 4 for d, v in days.items() if len(v) == 4} def nasa(path): rows = {} with open(path, encoding="utf-8", errors="replace") as f: for line in f: m = re.match(r"^(\d{4})-(\d{2})-(\d{2})\s+(.*)$", line.strip()) if m: rows[dt.date(int(m.group(1)), int(m.group(2)), int(m.group(3)))] = [float(x) for x in m.group(4).split()] return rows def gfz_ap(path): ap = {} with open(path, encoding="utf-8", errors="replace") as f: for line in f: if line.startswith("#") or not line.strip(): continue p = line.split() try: ap[dt.date(int(p[0]), int(p[1]), int(p[2]))] = int(p[23]) except (ValueError, IndexError): pass return ap def pct(sample, q): s = sorted(sample) k = q * (len(s) - 1) f = int(k) return s[f] + (s[min(f + 1, len(s) - 1)] - s[f]) * (k - f) def pearson(x, y): mx, my = st.mean(x), st.mean(y) return sum((a - mx) * (b - my) for a, b in zip(x, y)) / (sum((a - mx) ** 2 for a in x) * sum((b - my) ** 2 for b in y)) ** 0.5 def main(): folder = sys.argv[1] if len(sys.argv) > 1 else "." p = lambda name: os.path.join(folder, name) A6, A = cams(p("cams_ozone_monitoring_sh_ozone_area.csv")) T6, T = cams(p("cams_ozone_monitoring_sh_50hPa_temperature_minimum.csv")) NA = nasa(p("to3areas_2026_toms+omi+omps.txt")) NT = nasa(p("tmins_50_2026_merra2.txt")) ap = gfz_ap(p("Kp_ap_Ap_SN_F107_since_1932.txt")) print(f"CAMS area to {max(A6)}; NASA area to {max(d for d, v in NA.items() if v[0] != MISSING)}; Ap to {max(ap)}") print("\n1. Records for the date, 2026 (area, million km2)") rec = {} for d, v in A.items(): if 1979 <= d.year <= 2025 and not (d.month == 2 and d.day == 29): k = (d.month, d.day) if k not in rec or v > rec[k][0]: rec[k] = (v, d.year) hits = [d for d in sorted(A) if d.year == 2026 and rec.get((d.month, d.day), (0,))[0] > 0 and A[d] >= rec[(d.month, d.day)][0]] print(f" CAMS, daily mean at or above the 1979-2025 record: {len(hits)} days") for d in hits: r = rec[(d.month, d.day)] print(f" {d} {A[d]:.2f} (record {r[0]:.2f}, {r[1]})") rec6 = {} for t, v in A6.items(): if 1979 <= t.year <= 2025: k = (t.month, t.day, t.hour) if k not in rec6 or v > rec6[k][0]: rec6[k] = (v, t.year) six = [t for t in sorted(A6) if t.year == 2026 and rec6.get((t.month, t.day, t.hour), (0,))[0] > 0 and A6[t] >= rec6[(t.month, t.day, t.hour)][0]] print(f" CAMS, six-hourly values at or above the record for their step: {len(six)} values on {len({t.date() for t in six})} days") nh = [(d, v[0], v[7]) for d, v in sorted(NA.items()) if v[0] != MISSING and v[7] > 0 and v[0] >= v[7]] print(f" NASA, daily value at or above its 1979-2025 maximum: {len(nh)} days") for d, v, m in nh: print(f" {d} {v:.2f} (maximum {m:.2f})") print("\n2. Late October") for day in range(10, 25): d = dt.date(2026, 10, day) vals = sorted(((A[dt.date(y, 10, day)], y) for y in range(1979, 2027) if dt.date(y, 10, day) in A), reverse=True) print(f" CAMS {d}: {A[d]:.2f}, rank {[y for v, y in vals].index(2026) + 1} of {len(vals)} (largest {vals[0][0]:.2f}, {vals[0][1]})") win = [d for d in sorted(A) if d.year == 2026 and d >= dt.date(2026, 9, 17)] over = [d for d in win if A[d] > pct([A[dt.date(y, d.month, d.day)] for y in range(1979, 2024) if dt.date(y, d.month, d.day) in A], 0.95)] print(f" CAMS above the 1979-2023 95th percentile for the date, 17 Sep to {win[-1]}: {len(over)} of {len(win)} days") nw = [d for d, v in sorted(NA.items()) if d >= dt.date(2026, 9, 18) and v[0] != MISSING] print(f" NASA above its 90th percentile, 18 Sep to {nw[-1]}: {sum(NA[d][0] > NA[d][6] for d in nw)} of {len(nw)} days") print("\n3. Record-low 50 hPa minimum temperature for the date, 2026") trec = {} for d, v in T.items(): if 1979 <= d.year <= 2025 and not (d.month == 2 and d.day == 29): k = (d.month, d.day) if k not in trec or v < trec[k]: trec[k] = v th = [d for d in sorted(T) if d.year == 2026 and d.month in (8, 9, 10) and T[d] <= trec[(d.month, d.day)]] print(f" CAMS (daily mean), August-October: {len(th)} days: {[str(d) for d in th]}") an = [(d, v) for d, v in sorted(NT.items()) if v[0] != MISSING and (len(v) < 9 or v[8] == 0)] nth = [d for d, v in an if v[0] <= v[1] and d >= dt.date(2026, 8, 1)] print(f" NASA file, analysed days from 1 August (values from GEOS FP, climatology MERRA-2): {len(nth)} days") print("\n4. What moves the hole from one year to the next, 1990-2026") rows = [] for y in range(1990, 2027): a = st.mean(A[d] for d in A if d.year == y and dt.date(y, 9, 7) <= d <= dt.date(y, 10, 13)) t = st.mean(T[d] for d in T if d.year == y and dt.date(y, 9, 1) <= d <= dt.date(y, 10, 13)) k = st.mean(v for d, v in ap.items() if dt.date(y, 8, 1) <= d <= dt.date(y, 10, 13)) rows.append((y, a, t, k)) da = [rows[i][1] - rows[i - 1][1] for i in range(1, len(rows))] dtm = [rows[i][2] - rows[i - 1][2] for i in range(1, len(rows))] dap = [rows[i][3] - rows[i - 1][3] for i in range(1, len(rows))] print(f" {len(da)} changes: r(area, temperature) = {pearson(dtm, da):+.2f}; r(area, Ap) = {pearson(dap, da):+.2f}") print(f" levels: r(area, temperature) = {pearson([x[2] for x in rows], [x[1] for x in rows]):+.2f}; " f"r(area, Ap) = {pearson([x[3] for x in rows], [x[1] for x in rows]):+.2f}") for y, a, t, k in rows[-3:]: print(f" {y}: mean area {a:.2f} million km2, mean temperature {t:.2f} K, mean Ap {k:.1f}") print("\n5. How much of 2026's change the cold accounts for") a0, b = fit(dtm, da) pred = a0 + b * dtm[-1] res = [y - (a0 + b * x) for x, y in zip(dtm, da)] sd = st.stdev(res) print(f" line through the {len(da)} changes: {b:+.2f} million km2 per K (intercept {a0:+.2f})") print(f" 2025 to 2026: temperature {dtm[-1]:+.2f} K, area {da[-1]:+.2f}; the line predicts {pred:+.2f} ({pred / da[-1]:.0%}); " f"left over {res[-1]:+.2f}, {res[-1] / sd:.1f} times the standard deviation of the residuals ({sd:.2f})") print(f" what the line leaves over, against the Ap changes: r = {pearson(dap, res):+.2f}") l0, lb = fit([x[2] for x in rows], [x[1] for x in rows]) lres = [x[1] - (l0 + lb * x[2]) for x in rows] print(f" in levels (area against temperature, 1990-2026), 2026 is {lres[-1]:+.2f} above its line, " f"rank {sorted(lres, reverse=True).index(lres[-1]) + 1} of {len(lres)}; against Ap: r = {pearson([x[3] for x in rows], lres):+.2f}") means = sorted((x[3], x[0]) for x in rows) print(f" 2026 mean Ap {rows[-1][3]:.1f}: number {[y for _, y in means].index(2026) + 1} from the lowest of {len(means)} years") storms = [(str(d), v) for d, v in sorted(ap.items()) if dt.date(2026, 8, 1) <= d <= dt.date(2026, 10, 13) and v >= 30] print(f" days with daily Ap of 30 or more, 1 Aug-13 Oct 2026: {storms}") def fit(x, y): """Least-squares line y = a + b x; returns (a, b).""" mx, my = st.mean(x), st.mean(y) b = sum((p - mx) * (q - my) for p, q in zip(x, y)) / sum((p - mx) ** 2 for p in x) return my - b * mx, b if __name__ == "__main__": main()