What did radar and satellites show around a reported tornado?

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

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pyproj import CRS, Geod, Transformer

from usdata import build_query, cite_lockfile, get, pull, select_by_time, verify
from usdata.providers import load_adapter
from usdata.readers import RadarDecodeError

# 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")
print("Executed (UTC):", datetime.now(UTC).isoformat(timespec="seconds"))
print(
    "; ".join(
        f"{name} {version(name)}" for name in ("usdata", "pandas", "xradar", "xarray", "pyproj")
    )
)

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

started = time.perf_counter()
result = pull(manifest)
assert len(result.fetched) == 3
items = {item.asset.dataset_id: item for item in result.fetched}
assert len(items) == 3
print("\nPull seconds:", round(time.perf_counter() - started, 2))
print("Total source bytes:", sum(item.provenance.size for item in result.fetched))
pd.DataFrame(
    {
        "dataset": item.asset.dataset_id,
        "asset": item.asset.id,
        "bytes": item.provenance.size,
        "retrieved_utc": item.provenance.retrieved_at.isoformat(timespec="seconds"),
        "checksum": item.provenance.checksum[:23] + "...",
    }
    for item in result.fetched
)

columns = [
    "EVENT_ID",
    "EVENT_TYPE",
    "STATE",
    "CZ_NAME",
    "BEGIN_DATE_TIME",
    "END_DATE_TIME",
    "CZ_TIMEZONE",
    "BEGIN_LAT",
    "BEGIN_LON",
    "TOR_F_SCALE",
    "SOURCE",
]
frame = items["noaa:storm-events"].open_csv(usecols=columns)
rows = frame.loc[frame.EVENT_ID.eq("1184052")]
assert len(rows) == 1
report = rows.iloc[0]
if report.CZ_TIMEZONE != "CST-6":
    raise ValueError(f"Review source timezone before matching: {report.CZ_TIMEZONE}")
local_time = datetime.strptime(report.BEGIN_DATE_TIME, "%d-%b-%y %H:%M:%S")
event_time = local_time.replace(tzinfo=timezone(timedelta(hours=-6))).astimezone(UTC)
assert event_time == report.BEGIN_UTC.to_pydatetime(), "the reader derived a different UTC time"
lat, lon = float(report.BEGIN_LAT), float(report.BEGIN_LON)
assert np.isfinite([lat, lon]).all() and -90 <= lat <= 90 and -180 <= lon <= 180
for key, value in report.items():
    print(f"{key:16s} {value}")
print("\nReport start UTC:", event_time.isoformat())
del frame, rows
gc.collect();

matches = []
for dataset_id, params in [
    ("noaa:nexrad-level2", {"site": "KTLX"}),
    ("noaa:goes-abi", {"satellite": 16, "channel": 13}),
]:
    query = build_query(
        start=event_time - timedelta(minutes=5), end=event_time + timedelta(minutes=5), **params
    )
    with load_adapter(get(dataset_id)) as provider:
        candidates = provider.list_assets(query)
    selection = select_by_time(
        candidates, target=event_time, tolerance=timedelta(minutes=5), direction="nearest"
    )
    prior = select_by_time(
        candidates, target=event_time, tolerance=timedelta(minutes=5), direction="at_or_before"
    )
    assert selection.asset is not None, f"No scan start within tolerance: {dataset_id}"
    nearest = selection.asset
    delta = selection.offset_seconds
    assert nearest.id == items[dataset_id].asset.id, "Revisit the manifest selection"
    matches.append(
        {
            "dataset": dataset_id,
            "candidates": selection.candidate_count,
            "eligible": selection.eligible_count,
            "scan_start_utc": nearest.time.start.isoformat(),
            "scan_end_utc": nearest.time.end.isoformat(),
            "start_offset_seconds": delta,
            "prior_start_offset_seconds": prior.offset_seconds,
        }
    )
pd.DataFrame(matches)

radar_item = items["noaa:nexrad-level2"]
try:
    full_volume = radar_item.open()
except RadarDecodeError as error:
    print("Full-volume limitation:", error)
else:
    print("This decoder passed full-volume alignment checks:", full_volume.groups)
    del full_volume
radar = radar_item.open_nexrad(sweep=0)
sweep = radar["sweep_0"].to_dataset()
assert radar.attrs["usdata"]["sweeps"] == ["sweep_0"]
assert radar.attrs["usdata"]["provenance"]["checksum"] == radar_item.provenance.checksum
print("Sweep dimensions:", dict(sweep.sizes))
print("Elevation degrees:", float(sweep.sweep_fixed_angle))
print("Ray time range:", str(sweep.time.min().values), "to", str(sweep.time.max().values))
nearby = sweep.DBZH.sel(range=slice(None, 100_000))
assert nearby.attrs["units"] == "dBZ"
bearing, _, distance = Geod(ellps="WGS84").inv(
    float(radar.longitude), float(radar.latitude), lon, lat
)
print("Report ground distance from KTLX (km):", round(distance / 1000, 2))

scene = items["noaa:goes-abi"].open()
projection = scene.goes_imager_projection.attrs
transform = Transformer.from_crs("EPSG:4326", CRS.from_cf(projection), always_xy=True)
px, py = transform.transform(lon, lat)
height = projection["perspective_point_height"]
sx, sy = px / height, py / height
assert np.isfinite([sx, sy]).all()
assert float(scene.x.min()) <= sx <= float(scene.x.max())
assert float(scene.y.min()) <= sy <= float(scene.y.max())
ix = int(np.abs(scene.x.values - sx).argmin())
iy = int(np.abs(scene.y.values - sy).argmin())
point = scene.isel(x=ix, y=iy)
print(
    "Ground-projected nearest pixel:",
    {
        "x": ix,
        "y": iy,
        "DQF": int(point.DQF),
        "accepted_CMI_K": round(float(point.CMI), 1) if int(point.DQF) == 0 else None,
    },
)
half_window = 120_000 / height
roi = scene.isel(
    x=np.flatnonzero(abs(scene.x.values - sx) < half_window),
    y=np.flatnonzero(abs(scene.y.values - sy) < half_window),
)
brightness = roi.CMI.where(roi.DQF == 0)
assert brightness.attrs["units"] == "K" and int(brightness.count()) > 0
print("Crop dimensions:", dict(brightness.sizes), "DQF=0 pixels:", int(brightness.count()))
print(f"Coldest accepted pixel in the crop: {float(brightness.min()):.1f} K")

fig = plt.figure(layout="constrained")
ax_radar = fig.add_subplot(121, projection="polar")
ax_goes = fig.add_subplot(122)
mesh = ax_radar.pcolormesh(
    np.deg2rad(nearby.azimuth.values),
    nearby.range.values / 1000,
    nearby.values.T,
    cmap="turbo",
    shading="auto",
    vmin=-10,
    vmax=65,
)
ax_radar.scatter(
    np.deg2rad(bearing), distance / 1000, marker="*", s=150, c="white", edgecolors="black"
)
ax_radar.set_theta_zero_location("N")
ax_radar.set_theta_direction(-1)
ax_radar.set_ylim(0, 100)
ax_radar.tick_params(labelsize=7)
ax_radar.set_title("KTLX first sweep, slant range (km)\n04:40:53 UTC start", fontsize=9)
fig.colorbar(mesh, ax=ax_radar, label="Reflectivity (dBZ)", shrink=0.7)
mesh = ax_goes.pcolormesh(
    brightness.x, brightness.y, brightness, cmap="magma", shading="auto", vmin=190, vmax=290
)
ax_goes.scatter(sx, sy, marker="*", s=150, c="cyan", edgecolors="black")
ax_goes.grid(False)
ax_goes.tick_params(labelsize=7)
ax_goes.set(
    xlabel="East-west scan angle (rad)",
    ylabel="North-south scan angle (rad)",
)
ax_goes.set_title("GOES-16 C13, DQF = 0\n04:41:17 UTC start", fontsize=9)
fig.colorbar(mesh, ax=ax_goes, label="Brightness temperature (K)", shrink=0.7)
fig.suptitle("Tornado report 1184052 begins 7 May 2024 04:39 UTC (stars)")
plt.show()
plt.close(fig)
del radar, sweep, nearby, scene, roi, brightness, point
gc.collect();

assert verify(manifest) == []
for citation in cite_lockfile(manifest):
    print(citation.as_text())
    print()

again = pull(manifest)
assert again.from_lockfile and all(item.from_cache for item in again.fetched)
assert verify(manifest) == []
with tempfile.TemporaryDirectory(prefix="usdata-event-restore-") as directory:
    started = time.perf_counter()
    restored = pull(manifest, root=Path(directory))
    assert restored.from_lockfile and not any(item.from_cache for item in restored.fetched)
    assert verify(manifest, root=Path(directory)) == []
    assert [item.provenance.checksum for item in restored.fetched] == [
        item.provenance.checksum for item in result.fetched
    ]
    print("Empty-cache restoration seconds:", round(time.perf_counter() - started, 2))
print("All three inputs verified; cached and empty-cache locked restoration passed.")