"""exercise_1b.py Exercise I-B. Signal, or a change in what was watching. Separates a change in the hazard record from a change in the system that recorded it, using the United States tornado database, where the observing system changed on a documented schedule. Nothing in this exercise tests whether hazards have intensified. It tests whether a count can be read as a rate, which is a question about denominators and detection. python exercise_1b.py """ from __future__ import annotations import hashlib import io import sys import urllib.request from pathlib import Path import matplotlib import numpy as np import pandas as pd matplotlib.use("Agg") import matplotlib.pyplot as plt # noqa: E402 sys.path.insert(0, str(Path(__file__).resolve().parent)) import common as C # noqa: E402 SPC = "https://www.spc.noaa.gov/wcm/data/1950-2023_actual_tornadoes.csv" POPFILE = "us_population_1950_2023.csv" # columns yr and pop, frozen SNAPFILE = "spc_tornadoes_1950_2023.csv" RADAR = (1990, 1997) # WSR-88D national deployment window STRONG = 2 # F/EF2 and above YEARS = (1950, 2023) def freeze(): p = C.SNAP / SNAPFILE if p.exists(): return pd.read_csv(p, low_memory=False) print(" downloading the SPC tornado database ...", flush=True) with urllib.request.urlopen(SPC, timeout=900) as resp: raw = resp.read() df = pd.read_csv(io.BytesIO(raw), low_memory=False) keep = [c for c in ("om", "yr", "mo", "dy", "st", "mag", "inj", "fat", "slat", "slon", "len", "wid") if c in df.columns] df = df[keep] df.to_csv(p, index=False) return df def national_population(): """National population by year, read from a frozen file in snapshots. The exercise needs a denominator that moves with time. Where no frozen series sits in snapshots the exercise still runs, and the rate panel is reported as unavailable, with nothing substituted for it. Nothing here touches the network. """ p = C.SNAP / POPFILE if not p.exists(): return None pop = pd.read_csv(p) if not {"yr", "pop"} <= set(pop.columns) or len(pop) < 10: return None return pop[["yr", "pop"]].drop_duplicates("yr") def main(): C.ensure_dirs() df = freeze() digest = hashlib.sha256((C.SNAP / SNAPFILE).read_bytes()).hexdigest() df = df[(df.yr >= YEARS[0]) & (df.yr <= YEARS[1])] print("=" * 74) print("EXERCISE I-B: SIGNAL, OR A CHANGE IN WHAT WAS WATCHING") print(f" source NOAA Storm Prediction Center severe weather database") print(f" snapshot {SNAPFILE} sha256 {digest[:32]}") print(f" records {len(df):,} tornado reports, {YEARS[0]} to {YEARS[1]}") print("=" * 74) print("\n This exercise tests the observing system and the denominator.") print(" It does not test whether hazards have intensified.") out = {"records": int(len(df)), "sha256": digest, "years": list(YEARS), "radar_window": list(RADAR)} raw = df.groupby("yr").size() strong = df[df.mag >= STRONG].groupby("yr").size().reindex(raw.index, fill_value=0) print("\n[1] RAW ANNUAL COUNTS") dec = raw.groupby((raw.index // 10) * 10).mean() for d0, v in dec.items(): print(f" {d0}s mean {v:>7.0f} reports per year") out["raw_by_decade"] = {int(k): round(float(v), 1) for k, v in dec.items()} print(f"\n[2] RESTRICTED TO F/EF{STRONG} AND STRONGER") decs = strong.groupby((strong.index // 10) * 10).mean() for d0, v in decs.items(): print(f" {d0}s mean {v:>7.0f} reports per year") out["strong_by_decade"] = {int(k): round(float(v), 1) for k, v in decs.items()} pre = slice(YEARS[0], RADAR[0] - 1) post = slice(RADAR[1] + 1, YEARS[1]) r_pre, r_post = raw.loc[pre].mean(), raw.loc[post].mean() s_pre, s_post = strong.loc[pre].mean(), strong.loc[post].mean() print(f"\n[3] ACROSS THE RADAR DEPLOYMENT WINDOW {RADAR[0]} TO {RADAR[1]}") print(f" {'series':<34}{'before':>10}{'after':>10}{'change':>10}") print(f" {'all reports':<34}{r_pre:>10.0f}{r_post:>10.0f}" f"{100*(r_post-r_pre)/r_pre:>9.0f}%") print(f" {f'F/EF{STRONG} and stronger':<34}{s_pre:>10.0f}" f"{s_post:>10.0f}{100*(s_post-s_pre)/s_pre:>9.0f}%") out["radar_change"] = { "all_before": round(float(r_pre), 1), "all_after": round(float(r_post), 1), "all_pct": round(100 * (r_post - r_pre) / r_pre, 1), "strong_before": round(float(s_pre), 1), "strong_after": round(float(s_post), 1), "strong_pct": round(100 * (s_post - s_pre) / s_pre, 1)} print(" The weak reports carry the increase and the strong ones do not.") print(" Strong reports depend less on detection, and their ratings") print(" depend on damage surveys and on the 2007 change from the") print(" Fujita to the Enhanced Fujita scale.") pop = national_population() rate = None if pop is not None: pp = pop.set_index("yr")["pop"].reindex( range(YEARS[0], YEARS[1] + 1)).interpolate().ffill().bfill() rate = raw / (pp.loc[raw.index] / 1e6) print("\n[4] PER MILLION RESIDENTS") rdec = rate.groupby((rate.index // 10) * 10).mean() for d0, v in rdec.items(): print(f" {d0}s mean {v:>6.2f} reports per million residents") out["rate_by_decade"] = {int(k): round(float(v), 2) for k, v in rdec.items()} else: print("\n[4] PER MILLION RESIDENTS") print(f" No population series in snapshots ({POPFILE}). The rate") print(" panel is omitted, and no series is substituted for it.") out["rate_by_decade"] = None print("\n[5] WHAT THIS CANNOT DETERMINE") for line in ("whether tornado frequency changed", "whether tornado intensity changed", "whether the F to EF rating change in 2007 shifted counts", "whether population growth changed exposure or only " "observation"): print(f" - {line}") # ---- figure ----------------------------------------------------------- panels = 3 if rate is not None else 2 fig, axes = plt.subplots(panels, 1, figsize=(9, 3.1 * panels), sharex=True) series = [(raw, "All tornado reports", "#1f3864"), (strong, f"F/EF{STRONG} and stronger", "#b2182b")] if rate is not None: series.append((rate, "All reports per million residents", "#4d4d4d")) for ax, (s, t, col) in zip(np.atleast_1d(axes), series): ax.plot(s.index, s.values, color=col, lw=1.4) ax.axvspan(*RADAR, color="#fdd0a2", alpha=0.5) ax.set_title(t, loc="left", fontweight="bold", fontsize=11) ax.set_ylabel("reports") np.atleast_1d(axes)[0].text( RADAR[1] + 1, np.atleast_1d(axes)[0].get_ylim()[1] * 0.9, " WSR-88D deployment", fontsize=8, color="#8c510a") np.atleast_1d(axes)[-1].set_xlabel("year") fig.suptitle(f"{'Three' if rate is not None else 'Two'} readings of " "one record", fontsize=13, fontweight="bold") np.atleast_1d(axes)[-1].text( 0, -0.32, "Source: NOAA Storm Prediction Center severe weather " f"database, {YEARS[0]} to {YEARS[1]}, {len(df):,} reports. " "No panel tests hazard intensification.", transform=np.atleast_1d(axes)[-1].transAxes, fontsize=7, color="#666") fig.tight_layout() fig.savefig(C.FIG / "exercise_1b_signal.png", dpi=140, bbox_inches="tight") plt.close(fig) tbl = pd.DataFrame({"year": raw.index, "all_reports": raw.values, f"ef{STRONG}_plus": strong.values}) if rate is not None: tbl["per_million"] = rate.values C.write_handoff("exercise_1b", tbl, { "period": f"{YEARS[0]} to {YEARS[1]}", "intensity_threshold": f"magnitude >= {STRONG}", "radar_window": list(RADAR), "denominator": ("national population, interpolated between published " "estimates" if rate is not None else "unavailable, no frozen " "population series in snapshots"), "rating_scale_change": "the F scale became the EF scale in 2007 and " "this analysis does not adjust for it", "unrated_reports": "magnitude below zero or missing is treated as " "weak and stays in the all-reports series", }, [ "Population controls for exposure and for the number of potential " "observers at once. Which of the two does the rate panel actually " "isolate?", "The rating scale changed in 2007. Does the strong-report series " "cross that boundary cleanly?", "Reports cluster where people and roads are. Does a national total " "hide a spatial reporting bias that a state series would show?", "What would you need in order to say something about the hazard " "itself, and is any of it in this file?", ]) C.save_checkvalues("exercise_1b", out) print(f"\nwrote {C.FIG / 'exercise_1b_signal.png'}") if __name__ == "__main__": main()