For one Oklahoma tornado, do the two report archives agree on when and where it was, and what did radar, lightning, and the model analysis show at that place and time?

import gc
import subprocess
import sys
import time
from datetime import UTC, datetime, timedelta, timezone
from importlib.metadata import version
from pathlib import Path

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pyproj import Geod

from usdata import pull, verify

# One figure style for every usdata notebook, so previews look alike.
plt.rcParams.update(
    {
        "figure.figsize": (8, 4.5),
        "figure.dpi": 120,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "axes.grid": True,
        "grid.alpha": 0.3,
        "font.size": 10,
    }
)
manifest = Path("dataset.yaml")
notebook_started = time.perf_counter()
pd.set_option("display.width", 130)
pd.set_option("display.max_columns", 20)
print("Executed (UTC):", datetime.now(UTC).isoformat(timespec="seconds"))
print(
    "; ".join(
        f"{name} {version(name)}"
        for name in ("usdata", "pandas", "numpy", "xarray", "xradar", "eccodes", "pyproj")
    )
)

geod = Geod(ellps="WGS84")
RADIUS_KM = 25.0
HALF_DEGREE = 0.5

print(manifest.read_text().strip())

started = time.perf_counter()
result = pull(manifest)
pull_seconds = time.perf_counter() - started

pinned = pd.DataFrame(
    {
        "source": key,
        "dataset": items[0].asset.dataset_id,
        "assets": len(items),
        "bytes": sum(item.provenance.size or 0 for item in items),
        "first_asset": items[0].asset.id,
    }
    for key, items in result.by_source.items()
)
total_bytes = int(pinned.bytes.sum())
print(f"\nPull seconds: {pull_seconds:.1f}")
print(f"Assets pinned: {len(result.fetched)}")
print(f"Total pinned bytes: {total_bytes:,} ({total_bytes / 1e6:.1f} MB)")
pinned

columns = [
    "EVENT_ID",
    "EVENT_TYPE",
    "STATE",
    "CZ_NAME",
    "CZ_TIMEZONE",
    "BEGIN_DATE_TIME",
    "END_DATE_TIME",
    "BEGIN_LAT",
    "BEGIN_LON",
    "END_LAT",
    "END_LON",
    "TOR_F_SCALE",
    "TOR_LENGTH",
    "TOR_WIDTH",
    "INJURIES_DIRECT",
    "DEATHS_DIRECT",
    "SOURCE",
]
reports_item = result.one("reports")
events = reports_item.open_csv(usecols=columns)
print("Records in the 2024 archive:", len(events))
for derived in events.attrs["usdata"]["derived"]:
    print(f"{derived['column']} from {derived['source']}: {derived['rule']}")
    print(f"  rows the rule could not parse: {derived['unparsed']}")

(report,) = events.loc[events.EVENT_ID.eq("1184052")].to_dict("records")
if pd.isna(report["BEGIN_UTC"]) or pd.isna(report["END_UTC"]):
    raise ValueError(f"No UTC time was derived from {report['CZ_TIMEZONE']!r}")
REPORT_UTC = report["BEGIN_UTC"].to_pydatetime()
END_UTC = report["END_UTC"].to_pydatetime()
LAT, LON = float(report["BEGIN_LAT"]), float(report["BEGIN_LON"])
END_LAT, END_LON = float(report["END_LAT"]), float(report["END_LON"])
_, _, straight_m = geod.inv(LON, LAT, END_LON, END_LAT)

for key, value in report.items():
    print(f"{key:18s} {value}")
print()
print("Report start, local:", report["BEGIN_DATE_TIME"], report["CZ_TIMEZONE"])
print("Report start, UTC:  ", REPORT_UTC.isoformat())
print("Report end, UTC:    ", END_UTC.isoformat())
print(f"Reported duration:   {(END_UTC - REPORT_UTC).total_seconds() / 60:.0f} minutes")
print(f"Start to end on the ground: {straight_m / 1000:.2f} km")
print(
    f"Reported path length: {report['TOR_LENGTH']} miles = {report['TOR_LENGTH'] * 1.609344:.2f} km"
)
print(f"Reported path width: {report['TOR_WIDTH']} yards = {report['TOR_WIDTH'] * 0.9144:.0f} m")

del events
gc.collect();

spc_item = result.one("spc")
spc = spc_item.open_csv(dtype={"stf": "string", "f1": "string"})
CST = timezone(timedelta(hours=-6))
local_date = REPORT_UTC.astimezone(CST).strftime("%Y-%m-%d")
print("SPC rows in 2024:", len(spc), "| whole-track rows:", int(spc.sg.eq(1).sum()))

candidates = spc.loc[spc.date.eq(local_date) & spc.sg.eq(1)].copy()
_, _, candidate_m = geod.inv(
    np.full(len(candidates), LON),
    np.full(len(candidates), LAT),
    candidates.slon.values,
    candidates.slat.values,
)
candidates["km_from_report"] = np.round(candidate_m / 1000, 2)
candidates = candidates.sort_values("km_from_report")
print(f"Whole-track tornadoes dated {local_date}: {len(candidates)}")
print(
    candidates.head(3)[
        ["om", "st", "time", "tz", "mag", "slat", "slon", "len", "wid", "km_from_report"]
    ].to_string(index=False)
)

spc_match = candidates.iloc[0]
if int(spc_match["tz"]) != 3:
    raise ValueError(f"Review the SPC timezone code before converting: {spc_match['tz']}")
spc_utc = (
    datetime.strptime(f"{spc_match['date']} {spc_match['time']}", "%Y-%m-%d %H:%M:%S")
    .replace(tzinfo=CST)
    .astimezone(UTC)
)
agreement = pd.DataFrame(
    [
        {
            "quantity": "start time (UTC)",
            "storm_events": REPORT_UTC.isoformat(),
            "spc": spc_utc.isoformat(),
        },
        {
            "quantity": "start latitude (degrees north)",
            "storm_events": LAT,
            "spc": float(spc_match["slat"]),
        },
        {
            "quantity": "start longitude (degrees east)",
            "storm_events": LON,
            "spc": float(spc_match["slon"]),
        },
        {
            "quantity": "rating",
            "storm_events": report["TOR_F_SCALE"],
            "spc": f"EF{int(spc_match['mag'])}",
        },
        {
            "quantity": "path length (miles)",
            "storm_events": report["TOR_LENGTH"],
            "spc": float(spc_match["len"]),
        },
        {
            "quantity": "path width (yards)",
            "storm_events": report["TOR_WIDTH"],
            "spc": float(spc_match["wid"]),
        },
    ]
)
print()
print(
    f"Separation: {spc_match['km_from_report']:.2f} km, "
    f"{(spc_utc - REPORT_UTC).total_seconds() / 60:.0f} minutes"
)
del spc, candidates
gc.collect()
agreement

volumes = result.by_source["radar"]
for volume in volumes:
    offset = (volume.asset.time.start - REPORT_UTC).total_seconds()
    print(
        f"{volume.asset.id}  scan start {volume.asset.time.start:%H:%M:%S} UTC"
        f"  {offset:+.0f} s from the report"
    )

chosen = min(volumes, key=lambda item: abs((item.asset.time.start - REPORT_UTC).total_seconds()))
scan_offset = (chosen.asset.time.start - REPORT_UTC).total_seconds()
print("\nChosen:", chosen.asset.id)

radar = chosen.open_nexrad(sweep=0)
sweep = radar["sweep_0"].to_dataset()
radar_lat, radar_lon = float(radar.latitude), float(radar.longitude)
elevation = float(sweep.sweep_fixed_angle)
bearing, _, radar_m = geod.inv(radar_lon, radar_lat, LON, LAT)
effective_radius_km = 4 / 3 * 6371.0
range_km = radar_m / 1000
beam_km = (
    np.sqrt(
        range_km**2
        + effective_radius_km**2
        + 2 * range_km * effective_radius_km * np.sin(np.deg2rad(elevation))
    )
    - effective_radius_km
)
print(f"KTLX at {radar_lat:.4f} N {radar_lon:.4f} E, sweep elevation {elevation:.2f} degrees")
print(f"Report is {range_km:.2f} km from KTLX on a bearing of {bearing % 360:.1f} degrees")
print(f"Beam centre over the report: {beam_km:.2f} km above the radar (4/3 earth refraction)")
print("Sweep rays span", str(sweep.time.min().values), "to", str(sweep.time.max().values))

near = sweep.DBZH.sel(range=slice(None, 90_000))
azimuths = near.azimuth.values
ranges = near.range.values
order = np.argsort(azimuths)
azimuths, reflectivity = azimuths[order], near.values[order]


# Cell boundaries, so pcolormesh takes the mesh as given instead of inferring edges
# from coordinates that turn back on themselves as the antenna sweeps past north.
azimuth_step, range_step = float(np.diff(azimuths).mean()), float(np.diff(ranges).mean())
azimuth_edges = np.concatenate([azimuths - azimuth_step / 2, [azimuths[-1] + azimuth_step / 2]])
range_edges = np.concatenate([ranges - range_step / 2, [ranges[-1] + range_step / 2]])
ground = np.outer(np.ones_like(azimuth_edges), range_edges * np.cos(np.deg2rad(elevation)))
gate_lon, gate_lat, _ = geod.fwd(
    np.full(ground.shape, radar_lon),
    np.full(ground.shape, radar_lat),
    np.outer(azimuth_edges, np.ones_like(range_edges)),
    ground,
)
centre_ground = np.outer(np.ones_like(azimuths), ranges * np.cos(np.deg2rad(elevation)))
centre_lon, centre_lat, _ = geod.fwd(
    np.full(centre_ground.shape, radar_lon),
    np.full(centre_ground.shape, radar_lat),
    np.outer(azimuths, np.ones_like(ranges)),
    centre_ground,
)
_, _, gate_m = geod.inv(
    np.full(centre_ground.shape, LON), np.full(centre_ground.shape, LAT), centre_lon, centre_lat
)
gate_km = gate_m / 1000
in_radius = gate_km <= RADIUS_KM
print(
    f"\nUnits: {near.attrs['units']}; gates kept: {reflectivity.size}, "
    f"within {RADIUS_KM:.0f} km of the report: {int(in_radius.sum())}"
)
max_dbz = float(np.nanmax(np.where(in_radius, reflectivity, np.nan)))
print(f"Maximum reflectivity within {RADIUS_KM:.0f} km: {max_dbz:.1f} dBZ")
del radar, sweep, near
gc.collect();

rotation_rows = []
rotation_field = None
for grid_item in result.by_source["rotation"]:
    grid = grid_item.open()
    (variable,) = grid.data_vars
    latitudes = grid.latitude.values
    longitudes = grid.longitude.values - 360.0
    rows = np.flatnonzero(np.abs(latitudes - LAT) <= HALF_DEGREE)
    cols = np.flatnonzero(np.abs(longitudes - LON) <= HALF_DEGREE)
    box_lat = latitudes[rows.min() : rows.max() + 1]
    box_lon = longitudes[cols.min() : cols.max() + 1]
    box = grid[variable].values[rows.min() : rows.max() + 1, cols.min() : cols.max() + 1]
    mesh_lon, mesh_lat = np.meshgrid(box_lon, box_lat)
    _, _, cell_m = geod.inv(np.full(box.shape, LON), np.full(box.shape, LAT), mesh_lon, mesh_lat)
    cell_km = cell_m / 1000
    circle = cell_km <= RADIUS_KM
    peak = float(np.max(np.where(circle, box, -np.inf)))
    tied = np.argwhere(circle & (box == peak))
    peak_lat = float(box_lat[tied[:, 0]].mean())
    peak_lon = float(box_lon[tied[:, 1]].mean())
    _, _, peak_m = geod.inv(LON, LAT, peak_lon, peak_lat)
    stamp = grid_item.asset.time.start
    rotation_rows.append(
        {
            "window_utc": f"{(stamp - timedelta(minutes=30)):%H:%M} to {stamp:%H:%M}",
            "peak_in_25km": peak,
            "tied_cells": len(tied),
            "peak_lat": round(peak_lat, 3),
            "peak_lon": round(peak_lon, 3),
            "km_from_report": round(peak_m / 1000, 1),
            "at_report_cell": float(box[np.unravel_index(int(cell_km.argmin()), cell_km.shape)]),
        }
    )
    if stamp.minute == 42:
        rotation_field = (box_lat.copy(), box_lon.copy(), box.copy())
        rotation_units = grid[variable].attrs["units"]
        rotation_name = grid[variable].attrs["long_name"]
        registry_attrs = grid.attrs["usdata"]["registry_attrs"]
        rotation_variable = variable
    del grid, box, cell_km, circle, mesh_lon, mesh_lat
    gc.collect()

rotation = pd.DataFrame(rotation_rows)
print("Decoded variable:", rotation_variable)
print("long_name:", rotation_name)
print("units:", rotation_units)
print("Filled from the registry entry:", registry_attrs)
rotation

glm_columns = [
    "flash_lat",
    "flash_lon",
    "flash_time_offset_of_first_event",
    "flash_quality_flag",
]
frames = []
for glm_item in result.by_source["lightning"]:
    detections = glm_item.open()
    frames.append(detections[glm_columns].reset_coords()[glm_columns].to_dataframe())
    detections.close()
    del detections
flashes = pd.concat(frames, ignore_index=True)
_, _, flash_m = geod.inv(
    np.full(len(flashes), LON),
    np.full(len(flashes), LAT),
    flashes.flash_lon.values,
    flashes.flash_lat.values,
)
flashes["km_from_report"] = flash_m / 1000

nearby = flashes.loc[flashes.km_from_report <= RADIUS_KM]
box_flashes = flashes.loc[
    flashes.flash_lat.between(LAT - HALF_DEGREE, LAT + HALF_DEGREE)
    & flashes.flash_lon.between(LON - HALF_DEGREE, LON + HALF_DEGREE)
]
print(f"Flashes in the whole field of view, 04:29 to 04:49 UTC: {len(flashes):,}")
print(f"Within the half-degree box: {len(box_flashes):,}")
print(f"Within {RADIUS_KM:.0f} km of the report: {len(nearby):,}")
print("Quality flags among those:", nearby.flash_quality_flag.value_counts().to_dict())

per_minute = nearby.flash_time_offset_of_first_event.dt.floor("min").value_counts().sort_index()
mark = pd.Timestamp(REPORT_UTC).tz_localize(None)
before, after = per_minute.loc[per_minute.index < mark], per_minute.loc[per_minute.index >= mark]
print(f"\nMean rate in the ten minutes before the report: {before.mean():.1f} flashes per minute")
print(f"Mean rate from the report onward:                 {after.mean():.1f} flashes per minute")
print(f"Busiest minute: {per_minute.idxmax():%H:%M} UTC with {per_minute.max()} flashes")

fig, ax = plt.subplots(layout="constrained")
ax.bar(per_minute.index, per_minute.values, width=pd.Timedelta(minutes=0.8), color="#39568c")
ax.axvline(mark, color="crimson", linewidth=1.4)
ax.annotate(
    "report begins",
    xy=(mark, per_minute.max()),
    xytext=(4, -2),
    textcoords="offset points",
    color="crimson",
    fontsize=8,
)
ax.set(
    xlabel="UTC minute (2024-05-07)",
    ylabel="Flashes per minute",
    title=f"GOES-16 GLM flashes within {RADIUS_KM:.0f} km of report 1184052",
)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
plt.show()
plt.close(fig)
del frames, flashes
gc.collect();

environment = result.one("environment")
summary = environment.inspect()
inventory = pd.DataFrame(
    {
        "shortName": message.short_name,
        "typeOfLevel": message.type_of_level,
        "level": message.level,
        "units": message.units,
        "name": message.name,
    }
    for message in summary.grib2.messages
)
print(
    f"{summary.asset_id}: {summary.format.value}, {summary.size:,} bytes, {len(inventory)} messages"
)
print("\nEvery message the partial fetch brought back:")
print(inventory.to_string())

started = time.perf_counter()
fields = environment.open()
messages = fields.attrs["usdata"]["messages"]
cape_name = next(name for name, message in messages.items() if message["shortName"] == "cape")
helicity_name = next(name for name, message in messages.items() if message["shortName"] == "hlcy")
analysis_utc = environment.asset.time.start
lead = REPORT_UTC - analysis_utc
print(f"\nDecode seconds: {time.perf_counter() - started:.1f}")
print("Variable names, from attrs['usdata']['messages']:")
for name, message in messages.items():
    print(f"  {name}: message {message['file_index']} of the fetched file, step {message['step']}")
print("Analysis valid:", analysis_utc.isoformat())
print(f"Report follows the analysis by {lead} ({lead.total_seconds() / 60:.0f} minutes)")

latitudes = fields.latitude.values
stored = fields.longitude.values
longitudes = np.where(stored > 180, stored - 360, stored)
rows, cols = np.where((np.abs(latitudes - LAT) <= 0.6) & (np.abs(longitudes - LON) <= 0.75))
window = (slice(rows.min(), rows.max() + 1), slice(cols.min(), cols.max() + 1))
model_lat, model_lon = latitudes[window], longitudes[window]
_, _, model_m = geod.inv(
    np.full(model_lat.shape, LON), np.full(model_lat.shape, LAT), model_lon, model_lat
)
model_km = model_m / 1000
nearest = np.unravel_index(int(model_km.argmin()), model_km.shape)
circle = model_km <= RADIUS_KM
print(
    f"Nearest grid cell: {model_lat[nearest]:.4f} N {model_lon[nearest]:.4f} E, "
    f"{model_km[nearest]:.2f} km from the report (grid spacing is 3 km); "
    f"{int(circle.sum())} cells within {RADIUS_KM:.0f} km"
)

environment_rows = []
for label, variable in (
    ("surface CAPE", cape_name),
    ("0-3 km storm-relative helicity", helicity_name),
):
    values = fields[variable].values[window]
    environment_rows.append(
        {
            "field": label,
            "variable": variable,
            "units": fields[variable].attrs["units"],
            "at_report_cell": round(float(values[nearest]), 1),
            "median_in_25km": round(float(np.median(values[circle])), 1),
            "max_in_25km": round(float(values[circle].max()), 1),
        }
    )
del fields
gc.collect()
pd.DataFrame(environment_rows)

box_lat, box_lon, box_values = rotation_field
fig, axes = plt.subplots(1, 2, layout="constrained")

ax = axes[0]
mesh = ax.pcolormesh(
    gate_lon, gate_lat, reflectivity, cmap="turbo", shading="flat", vmin=-10, vmax=65
)
ax.scatter(
    box_flashes.flash_lon,
    box_flashes.flash_lat,
    s=4,
    c="white",
    alpha=0.55,
    linewidths=0,
    label="GLM flash centroid",
)
ax.plot([LON, END_LON], [LAT, END_LAT], color="black", linewidth=1.6, label="Reported path")
ax.scatter(
    LON, LAT, marker="*", s=200, c="magenta", edgecolors="black", zorder=5, label="Report start"
)
ax.scatter(radar_lon, radar_lat, marker="^", s=70, c="black", zorder=5, label="KTLX")
ax.legend(loc="lower left", fontsize=6.5, facecolor="0.8")
ax.grid(False)
ax.set(
    xlim=(LON - HALF_DEGREE, LON + HALF_DEGREE),
    ylim=(LAT - HALF_DEGREE, LAT + HALF_DEGREE),
    xlabel="Longitude (degrees east)",
    ylabel="Latitude (degrees north)",
    title=f"KTLX {elevation:.2f}° sweep, {chosen.asset.time.start:%H:%M:%S} UTC",
    aspect=1 / np.cos(np.deg2rad(LAT)),
)
fig.colorbar(mesh, ax=ax, label="Reflectivity (dBZ)", shrink=0.8)

ax = axes[1]
mesh = ax.pcolormesh(box_lon, box_lat, box_values, cmap="magma", shading="nearest", vmin=0, vmax=20)
ax.scatter(LON, LAT, marker="*", s=200, c="cyan", edgecolors="black", zorder=5)
ax.plot([LON, END_LON], [LAT, END_LAT], color="cyan", linewidth=1.2)
ax.grid(False)
ax.set(
    xlim=(LON - HALF_DEGREE, LON + HALF_DEGREE),
    ylim=(LAT - HALF_DEGREE, LAT + HALF_DEGREE),
    xlabel="Longitude (degrees east)",
    title="MRMS mid-level shear, 04:12 to 04:42 UTC",
    aspect=1 / np.cos(np.deg2rad(LAT)),
)
fig.colorbar(mesh, ax=ax, label=f"Azimuthal shear ({rotation_units})", shrink=0.8)
for ax in axes:
    ax.title.set_fontsize(9)
fig.suptitle("Tornado report 1184052, Oklahoma County, 7 May 2024")
plt.show()
plt.close(fig)
del gate_lon, gate_lat, centre_lon, centre_lat, reflectivity, gate_km, in_radius
gc.collect();

answer = pd.DataFrame(
    [
        {
            "source": "reports",
            "dataset": "noaa:storm-events",
            "what it gave": "event 1184052: EF1, 04:39:00 to 04:42:00 UTC, "
            f"{LAT} N {LON} E, {report['TOR_LENGTH']} mi by {report['TOR_WIDTH']} yd, "
            f"{report['INJURIES_DIRECT']} injuries",
        },
        {
            "source": "spc",
            "dataset": "noaa:spc-tornado-reports",
            "what it gave": f"tornado {int(spc_match['om'])}: same start time, "
            f"{spc_match['km_from_report']:.2f} km from the same point, "
            f"EF{int(spc_match['mag'])}, {spc_match['len']} mi by {int(spc_match['wid'])} yd",
        },
        {
            "source": "radar",
            "dataset": "noaa:nexrad-level2",
            "what it gave": f"KTLX sweep starting {scan_offset:+.0f} s from the report, "
            f"beam {beam_km:.2f} km up, "
            f"{max_dbz:.1f} dBZ maximum within {RADIUS_KM:.0f} km",
        },
        {
            "source": "rotation",
            "dataset": "noaa:mrms",
            "what it gave": f"mid-level shear peak {rotation.peak_in_25km.iloc[-1]:.0f} "
            f"({rotation_units}) in the 04:12 to 04:42 window, "
            f"{rotation.km_from_report.iloc[-1]:.1f} km from the reported start",
        },
        {
            "source": "lightning",
            "dataset": "noaa:goes-glm",
            "what it gave": f"{len(nearby):,} flashes within {RADIUS_KM:.0f} km in 20 minutes, "
            f"{before.mean():.1f} per minute before the report and {after.mean():.1f} after",
        },
        {
            "source": "environment",
            "dataset": "noaa:hrrr",
            "what it gave": f"04:00 UTC analysis, {lead.total_seconds() / 60:.0f} minutes early: "
            f"{environment_rows[0]['at_report_cell']:.0f} J/kg surface CAPE and "
            f"{environment_rows[1]['at_report_cell']:.0f} m2/s2 0-3 km helicity at the report cell",
        },
    ]
)
pd.set_option("display.max_colwidth", 150)
answer

assert verify(manifest) == []
usdata_cli = Path(sys.executable).with_name("usdata")
for style in ("text", "bibtex"):
    finished = subprocess.run(
        [str(usdata_cli), "cite", str(manifest), "--format", style],
        capture_output=True,
        text=True,
        check=True,
    )
    print(f"$ usdata cite {manifest} --format {style}")
    print(finished.stdout.rstrip())
    print()

again = pull(manifest)
assert again.from_lockfile and all(item.from_cache for item in again.fetched)
assert verify(manifest) == []
print("Restored from the lockfile:", again.from_lockfile)
print(f"Notebook seconds: {time.perf_counter() - notebook_started:.0f}")