How warm was Oklahoma City on 6 May 2024 compared with its hourly normals?

from datetime import UTC, datetime
from pathlib import Path
from tempfile import TemporaryDirectory

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import pandas as pd

import usdata
from usdata import cite_lockfile, 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")
print("Executed (UTC):", datetime.now(UTC).isoformat(timespec="seconds"))
print(f"usdata {usdata.__version__}; pandas {pd.__version__}")

print(manifest.read_text())

result = pull(manifest)
observations = result.one("observations").open_csv(dtype={"DATE": "string"})
normals = result.one("normals").open_csv(dtype={"DATE": "string"})
print(f"Fetched {len(observations)} LCD reports and {len(normals)} hourly normals.")
print(f"Lockfile pins {len(result.lockfile.assets)} assets.")
print()
for item in result.fetched:
    print(item.asset.dataset_id)
    print(f"Bytes: {item.provenance.size}; cache hit: {item.from_cache}")
    print("Source:", item.provenance.source_url)
    print("Retrieved (UTC):", item.provenance.retrieved_at.isoformat(timespec="seconds"))
    print("Checksum:", item.provenance.checksum)
    print()

assert set(observations["STATION"]) == {"72353013967"}
assert set(normals["STATION"]) == {"USW00013967"}
reports = observations.loc[observations["REPORT_TYPE"].str.strip().eq("FM-15")].copy()
reports["observed_at"] = pd.to_datetime(reports["DATE"], format="%Y-%m-%dT%H:%M:%S")
# Flagged strings stay invalid; do not silently strip quality suffixes.
reports["observed_c"] = pd.to_numeric(reports["HourlyDryBulbTemperature"], errors="coerce")
reports["observed_c"] = reports["observed_c"].mask(reports["observed_c"].eq(-9999))
normals["normal_at"] = pd.to_datetime("2024-" + normals["DATE"], format="%Y-%m-%dT%H:%M:%S")
normals["normal_c"] = pd.to_numeric(normals["HLY-TEMP-NORMAL"], errors="coerce")
normals["normal_c"] = normals["normal_c"].mask(normals["normal_c"].eq(-9999))
assert not normals["normal_at"].duplicated().any()
comparison = pd.merge_asof(
    reports[["observed_at", "observed_c"]].sort_values("observed_at"),
    normals[["normal_at", "normal_c"]].sort_values("normal_at"),
    left_on="observed_at",
    right_on="normal_at",
    direction="nearest",
    tolerance=pd.Timedelta(minutes=10),
)
comparison["departure_c"] = comparison["observed_c"] - comparison["normal_c"]
comparison["offset_minutes"] = (
    comparison["normal_at"] - comparison["observed_at"]
).dt.total_seconds() / 60
usable = comparison.dropna(subset=["departure_c"])
assert len(usable) > 0, "No usable matched temperatures; inspect the source tables."
print(f"Routine reports: {len(reports)}; other report types: {len(observations) - len(reports)}")
print(f"Invalid observations: {comparison.observed_c.isna().sum()}")
print(f"Invalid normals: {normals.normal_c.isna().sum()}")
print(f"Unmatched reports: {comparison.normal_at.isna().sum()}; usable pairs: {len(usable)}")
comparison.tail()

fig, (temps, departures) = plt.subplots(2, 1, sharex=True, layout="constrained")
temps.plot(
    usable["observed_at"],
    usable["observed_c"],
    marker=".",
    color="#b45631",
    label="Routine observation",
)
temps.plot(
    usable["observed_at"],
    usable["normal_c"],
    color="#2563a6",
    label="Nearest hourly normal (within 10 min)",
)
temps.set_ylabel("Temperature (°C)")
temps.set_title("Will Rogers World Airport, Oklahoma City, 6 May 2024")
temps.legend(loc="upper left", frameon=False)
departures.bar(usable["observed_at"], usable["departure_c"], width=0.03, color="#b45631")
departures.axhline(0, color="#444444", linewidth=0.8)
departures.set_ylabel("Departure (°C)")
departures.set_xlabel("Report time, local standard time (CST, UTC-6)")
departures.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
plt.show()

warm = usable["departure_c"].gt(0).sum()
peak = usable.loc[usable["departure_c"].idxmax()]
print(f"{warm} of {len(usable)} usable reports were warmer than the matched normal.")
print(f"Mean departure across matched reports: {usable.departure_c.mean():.2f} C")
print(f"Largest warm departure: {peak.departure_c:.1f} C at {peak.observed_at:%Y-%m-%d %H:%M} CST")
print(
    f"Matched normal: {peak.normal_at:%Y-%m-%d %H:%M} CST ({peak.offset_minutes:.0f} minutes later)"
)

assert verify(manifest) == []
with TemporaryDirectory(prefix="usdata-hourly-restore-") as fresh:
    restored = pull(manifest, root=fresh)
    assert restored.from_lockfile
    assert all(not item.from_cache for item in restored.fetched)
    assert restored.lockfile == result.lockfile
    assert verify(manifest, root=fresh) == []
    for name in ("observations", "normals"):
        assert restored.one(name).path.read_bytes() == result.one(name).path.read_bytes()
print("Both inputs restored byte-for-byte into an empty cache and passed verification.")
print()
for citation in cite_lockfile(manifest):
    print(citation.as_text())