How high was Hurricane Helene's storm surge at Cedar Key?

import math
from datetime import UTC, datetime
from pathlib import Path

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import pandas as pd
from IPython.display import display

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)
items = {item.asset.dataset_id: item for item in result.fetched}
assert len(items) == len(result.fetched) == 3, "expected exactly one asset per source"


def series(dataset_id: str) -> pd.DataFrame:
    """Open one CO-OPS CSV with stripped column names and a UTC time index."""
    frame = items[dataset_id].open().rename(columns=str.strip)
    frame["Date Time"] = pd.to_datetime(frame["Date Time"], utc=True)
    return frame.set_index("Date Time")


observed = series("noaa:coops-water-levels")
predicted = series("noaa:coops-tide-predictions")
assert observed.index.equals(predicted.index), "observed and predicted rows must share timestamps"
print(f"Rows: {len(observed)} observed, {len(predicted)} predicted")
print(f"Observed quality flags: {observed['Quality'].value_counts().to_dict()}")
display(observed[["Water Level", "Sigma", "Quality"]].head(3))
display(predicted[["Prediction"]].head(3))

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()

levels = pd.DataFrame({"observed": observed["Water Level"], "predicted": predicted["Prediction"]})
levels["residual"] = levels["observed"] - levels["predicted"]
assert not levels.isna().any().any(), "gaps in the record would need explicit handling"
peak_residual_at = levels["residual"].idxmax()
peak_level_at = levels["observed"].idxmax()
peak = levels.loc[peak_residual_at]
highest = levels.loc[peak_level_at]
print(f"Peak residual: {peak['residual']:.2f} m at {peak_residual_at:%Y-%m-%d %H:%M} UTC")
print(
    f"Peak observed level: {highest['observed']:.2f} m MLLW at {peak_level_at:%Y-%m-%d %H:%M} UTC"
)
print(f"Predicted tide at the peak residual: {peak['predicted']:.2f} m MLLW")
around = slice(
    peak_residual_at - pd.Timedelta(minutes=18), peak_residual_at + pd.Timedelta(minutes=18)
)
levels.loc[around].round(3)

fig, (water, surge) = plt.subplots(2, 1, sharex=True, layout="constrained")
water.plot(levels.index, levels["observed"], color="#b45631", linewidth=1.4, label="Observed")
water.plot(
    levels.index, levels["predicted"], color="#2563a6", linewidth=1.1, label="Predicted tide"
)
water.set_title("Cedar Key, Florida (NOAA station 8727520), Hurricane Helene")
water.set_ylabel("Water level (m MLLW)")
water.legend(loc="upper left", frameon=False)
surge.fill_between(levels.index, 0, levels["residual"], color="#b45631", alpha=0.35)
surge.plot(levels.index, levels["residual"], color="#b45631", linewidth=1.2)
surge.axvline(peak_residual_at, color="#444444", linewidth=0.8, linestyle="--")
surge.set_title("Residual: observed minus predicted")
surge.set_ylabel("Residual (m)")
surge.set_xlabel("UTC")
surge.xaxis.set_major_formatter(mdates.DateFormatter("%d %b\n%H:%M", tz=UTC))
plt.show()

GAUGE = (29.135, -83.032)


def distance_km(lat: float, lon: float) -> float:
    """Great-circle distance from the gauge to a point, in kilometres."""
    lat1, lon1, lat2, lon2 = map(math.radians, (GAUGE[0], GAUGE[1], lat, lon))
    half_dlat, half_dlon = (lat2 - lat1) / 2, (lon2 - lon1) / 2
    a = math.sin(half_dlat) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(half_dlon) ** 2
    return 2 * 6371.0 * math.asin(math.sqrt(a))


tracks = items["noaa:hurdat2"].open()
helene = tracks[(tracks["name"] == "HELENE") & (tracks["time"].dt.year == 2024)].copy()
assert len(helene) > 0, "the HURDAT2 revision should include the 2024 season"
helene["distance_km"] = [
    distance_km(lat, lon) for lat, lon in zip(helene["latitude"], helene["longitude"], strict=True)
]
near = (helene["time"] >= peak_residual_at - pd.Timedelta(hours=18)) & (
    helene["time"] <= peak_residual_at + pd.Timedelta(hours=6)
)
columns = ["time", "record_identifier", "status", "latitude", "longitude"]
columns += ["max_wind_kt", "min_pressure_mb", "distance_km"]
window = helene[near]
display(window[columns].assign(distance_km=window["distance_km"].round(1)))

landfalls = helene[helene["record_identifier"] == "L"]
landfall = landfalls.loc[landfalls["distance_km"].idxmin()]
before = helene[helene["time"] <= peak_residual_at].iloc[-1]
after = helene[helene["time"] >= peak_residual_at].iloc[0]
span = after["time"] - before["time"]
fraction = 0.0 if span == pd.Timedelta(0) else (peak_residual_at - before["time"]) / span
storm_lat = before["latitude"] + fraction * (after["latitude"] - before["latitude"])
storm_lon = before["longitude"] + fraction * (after["longitude"] - before["longitude"])
lag_hours = (peak_residual_at - landfall["time"]).total_seconds() / 3600
side = "east" if GAUGE[1] > storm_lon else "west"
winds = "onshore" if GAUGE[1] > storm_lon else "offshore"
print(
    f"Landfall record: {landfall['time']:%Y-%m-%d %H:%M} UTC at "
    f"{landfall['latitude']:.1f}N {abs(landfall['longitude']):.1f}W, "
    f"{landfall['max_wind_kt']:.0f} kt, {landfall['distance_km']:.0f} km from the gauge"
)
print(f"Peak residual came {lag_hours:+.1f} h relative to landfall")
print(
    f"Interpolated storm position at the peak: {storm_lat:.2f}N {abs(storm_lon):.2f}W, "
    f"{distance_km(storm_lat, storm_lon):.0f} km from the gauge"
)
print(f"The gauge lay {side} of the track, on the {winds}-wind side of a northward-moving storm")

print(
    f"Helene added {peak['residual']:.2f} m above the predicted tide at Cedar Key, "
    f"peaking at {peak_residual_at:%H:%M} UTC on {peak_residual_at:%d %B %Y}."
)
print(
    f"The water stood {highest['observed']:.2f} m above MLLW at its highest, "
    f"{highest['residual']:.2f} m of which was residual."
)
print(
    f"The peak residual came {lag_hours:.1f} h after the best-track landfall point, "
    f"{landfall['distance_km']:.0f} km away."
)
print("Scope: one verified NOAA gauge, one harmonic prediction, one best-track revision.")

assert verify(manifest) == []
again = pull(manifest)
print(f"Restored from lockfile: {again.from_lockfile}")
print(f"All assets came from cache: {all(item.from_cache for item in again.fetched)}")
print()
for citation in cite_lockfile(manifest):
    print(citation.as_text())