"""Rebuild the draft's figures from retained research files. No downloads."""
from __future__ import annotations

import datetime as dt
import csv
import gzip
import hashlib
import json
import sys
from pathlib import Path
from zoneinfo import ZoneInfo

import matplotlib

matplotlib.use("Agg")
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Patch


ROOT = Path(__file__).resolve().parents[3]
sys.path.insert(0, str(ROOT / "scripts/research"))
import nov1_rain_long_windows as water

OUT = Path(__file__).resolve().parent
RESEARCH = ROOT / "data/research/nov1-rain-pulse-20261105"
SOURCE_NAMES = [
    "rain-pulse-analogues.json",
    "pulse-forecast-calibration.json",
    "tide-phase-screen.json",
    "long-window-study.json",
    "nov15-outlook-snapshot.json.gz",
    "nov6-afd-snapshot.json",
]
sources = {name: RESEARCH / name for name in SOURCE_NAMES}
comparison_path = RESEARCH / "local-gauge-check/comparison.json"
comparison = json.loads(comparison_path.read_text())
by_event = {row["event"]: row for row in comparison["events"]}
analogues = json.loads(sources["rain-pulse-analogues.json"].read_text())
calibration = json.loads(sources["pulse-forecast-calibration.json"].read_text())
phase = json.loads(sources["tide-phase-screen.json"].read_text())
long_window = json.loads(sources["long-window-study.json"].read_text())
PACIFIC = ZoneInfo("America/Los_Angeles")
TEAL, ORANGE, NAVY, LIGHT, GRAY = "#176c69", "#c05a21", "#102d3a", "#dfeaed", "#779096"
plt.rcParams.update({
    "font.family": "DejaVu Sans", "font.size": 10,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.edgecolor": "#9bb5b8", "axes.labelcolor": NAVY,
    "text.color": NAVY, "xtick.color": "#36535e", "ytick.color": "#36535e",
    "savefig.facecolor": "white",
})


def save(fig, stem):
    fig.savefig(OUT / f"{stem}.png", dpi=175, bbox_inches="tight")
    fig.savefig(OUT / f"{stem}.svg", bbox_inches="tight")
    plt.close(fig)


# Conceptual timing diagram: observed rain may start as the tide falls and continue
# across low tide. The pulse shape illustrates the three shorter analogues;
# it is deliberately not drawn as a measured November 1 time series.
x = np.linspace(-8, 4, 241)
rain_x = np.array([-7.5, -6.8, -6.1, -5.4, -4.7, -4, -3.3, -2.6, -1.9, -1.2])
rain_y = np.array([.18, .28, .12, .38, .68, .86, .72, .53, .27, .12])
tide = np.cos(np.pi * x / 6)
rise_progress = np.clip((x + 4) / 2.7, 0, 1)
fall_progress = np.clip((2.3 - x) / 3.6, 0, 1)
smooth = lambda value: value * value * (3 - 2 * value)
pulse = np.where(x <= -1.3, smooth(rise_progress), smooth(fall_progress))
fig, (rain_ax, tide_ax, pulse_ax) = plt.subplots(
    3, 1, figsize=(10.8, 6.8), sharex=True,
    gridspec_kw={"height_ratios": [1.0, 1.2, 1.15], "hspace": .18},
)
fig.subplots_adjust(top=.86, bottom=.13, left=.13, right=.97)
fig.suptitle("Rain can cross the tide turn before the extra rise peaks", x=.13, y=.97,
             ha="left", fontsize=16, weight="bold")
fig.text(.13, .905, "Timing hypothesis  ·  schematic, not a measured storm", fontsize=10, color=GRAY)
rain_ax.bar(rain_x, rain_y, width=.58, color=ORANGE, alpha=.88)
rain_ax.set(ylabel="Local rain", ylim=(0, 1.15), yticks=[])
rain_ax.text(-4.3, 1.03, "Rain continues as tide rises", color=ORANGE, fontsize=9)
tide_ax.plot(x, tide, color=NAVY, lw=2.5)
tide_ax.set(ylabel="Predicted tide", ylim=(-1.18, 1.35), yticks=[])
tide_ax.annotate("Incoming tide", xy=(-2.7, .18), xytext=(-5.2, .88),
                 arrowprops={"arrowstyle": "->", "color": NAVY}, fontsize=9)
tide_ax.annotate("Outgoing tide", xy=(2, .35), xytext=(1, -.25),
                 arrowprops={"arrowstyle": "->", "color": NAVY}, fontsize=9)
pulse_ax.plot(x, pulse, color=TEAL, lw=2.5)
pulse_ax.fill_between(x, pulse, color=TEAL, alpha=.15)
pulse_ax.scatter([-1.3], [1], color=TEAL, s=50, zorder=3)
pulse_ax.annotate("Extra-water peak\n1–2h before high", xy=(-1.3, 1), xytext=(-4.3, .92),
                  arrowprops={"arrowstyle": "->", "color": TEAL}, fontsize=9, color=TEAL)
pulse_ax.text(.45, .63, "Short pulses fade\nas tide falls", fontsize=9, color=TEAL)
pulse_ax.set(ylabel="Extra residual", ylim=(0, 1.35), yticks=[])
for axis in (rain_ax, tide_ax, pulse_ax):
    axis.axvline(-6, color=GRAY, lw=1, ls=":")
    axis.axvline(0, color=GRAY, lw=1, ls=":")
    axis.set_xlim(-8, 4)
    axis.grid(axis="x", color=LIGHT, lw=.7)
pulse_ax.set_xticks([-8, -6, -4, -2, 0, 2, 4],
                   ["−8", "Low tide\n−6", "−4", "−2", "High tide\n0", "+2", "+4"])
fig.text(.13, .035, "Three shorter pulses faded 1–2½h after high. The November 15, 2025 pulse lasted much longer.",
         fontsize=8.5, color="#36535e")
save(fig, "rain-across-tide-turn")


# Fixed elapsed-hour bins are deliberately distinguished from the moving one-hour peak.
nov1 = next(x for x in analogues["profiles"] if x["event"] == "2026-11-01")
hours = list(range(len(nov1["oneHourMeanChangeByElapsedHourInches"])))
fig, (ax, bars) = plt.subplots(2, 1, figsize=(10.8, 7), sharex=True,
                               gridspec_kw={"height_ratios": [3.2, 1.25], "hspace": .12})
ax.plot(hours, nov1["oneHourMeanChangeByElapsedHourInches"], marker="o", ms=3,
        color=TEAL, lw=2.1, label="Fixed one-hour bins")
ax.axhline(2, color=ORANGE, ls="--", lw=1.2, label="+2 in watch threshold")
ax.axhline(0, color=GRAY, lw=.8)
ax.axvspan(6.88, 11.38, color=ORANGE, alpha=.1, label="Sustained pulse")
ax.axvline(7.8, color=ORANGE, lw=1.3)
ax.scatter([7.8], [4.501], color=ORANGE, zorder=5)
ax.annotate("Rolling 1h peak: +4.5 in\n2:55 a.m. PST", (7.8, 4.501),
            xytext=(10.3, 5.3), arrowprops={"arrowstyle": "-", "color": ORANGE},
            color=ORANGE, fontsize=9)
ax.axvline(8.97, color=NAVY, ls=":", lw=1.3)
ax.text(9.2, -2.4, "4:05 a.m. PST\npredicted high tide", fontsize=9, va="bottom")
ax.set(xlim=(-.4, 25), ylim=(-3.1, 6.2), ylabel="Residual change from pre-rain 12h mean (in)",
       title="November 1: a short pulse inside a modest daily change")
ax.grid(axis="y", color=LIGHT)
ax.legend(frameon=False, loc="upper right", fontsize=8.5)
nov1_window = next(x for x in long_window["events"] if x["event"] == "2026-11-01")
windows = nov1_window["windows"]
vals = [windows[k]["changeVsPre12hInches"] for k in ["post0to12h", "post12to24h", "post0to24h"]]
bars.bar([6, 18, 28], vals, width=[9, 9, 7], color=[TEAL, GRAY, NAVY])
for x, val in zip([6, 18, 28], vals):
    bars.text(x, val + (.15 if val >= 0 else -.45), f"{val:+.2f} in", ha="center", fontsize=9)
bars.set(xlim=(-.4, 32), ylim=(-1.0, 2.4), ylabel="Window mean change (in)",
         xlabel="Elapsed hours after first creek rain report, 8:07 p.m. PDT Oct 31")
bars.set_xticks([0, 6, 12, 18, 24, 28], ["0", "6", "12", "18", "24", "24h mean"])
bars.grid(axis="y", color=LIGHT)
fig.text(.07, .02, "One-hour curve uses fixed elapsed bins; peak marker uses a moving one-hour mean. Nov 1 crossed the DST change.", fontsize=8, color="#36535e")
save(fig, "november-1-pulse")


# The airport's small November 1 total is real at SNA, but it does not
# describe the rainfall at the nearby Corona del Mar and upper-bay gauges.
def cumulative_counter(path, field, origin, end):
    with gzip.open(path, "rt") as stream:
        records = []
        for row in csv.DictReader(stream):
            if not row.get(field):
                continue
            time = dt.datetime.strptime(row["utc_valid"], "%Y-%m-%d %H:%M:%S").replace(tzinfo=dt.timezone.utc)
            if origin-dt.timedelta(hours=2) <= time <= end:
                records.append((time, float(row[field])))
    records.sort()
    baseline = [v for t,v in records if t < origin][-1]
    running = baseline
    points = [(0,0)]
    for time,value in records:
        if time >= origin:
            running = max(running,value)
            points.append(((time-origin).total_seconds()/3600, running-baseline))
    return points


rain_origin = dt.datetime.fromisoformat(by_event["2026-11-01"]["coronaDelMarFirstNearbyIncrementUTC"])
rain_end = rain_origin + dt.timedelta(hours=12)
cdm_path = RESEARCH / "local-gauge-check/2026-11-01-cdmc1.csv.gz"
sanc_path = RESEARCH / "sanc1-20261020-1104.csv.gz"
sna_path = ROOT / "data/research/data-page-refresh-20261105/airport.gz"
cdm_points = cumulative_counter(cdm_path,"PCIRRZZ",rain_origin,rain_end)
sanc_points = cumulative_counter(sanc_path,"PCIRRZZ",rain_origin,rain_end)
sna_points = [(0,0)];sna_total=0
with gzip.open(sna_path,"rt") as stream:
    for row in csv.DictReader(stream):
        time=dt.datetime.strptime(row["valid"],"%Y-%m-%d %H:%M").replace(tzinfo=dt.timezone.utc)
        if rain_origin <= time <= rain_end and row["p01i"] not in ("","M","T"):
            sna_total += float(row["p01i"])
            sna_points.append(((time-rain_origin).total_seconds()/3600,sna_total))
fig,(ax,clock)=plt.subplots(1,2,figsize=(11.5,5),gridspec_kw={"width_ratios":[1.15,1],"wspace":.34})
for label,points,color in [("Corona del Mar",cdm_points,ORANGE),("SANC1, Campus Drive",sanc_points,TEAL),("John Wayne Airport",sna_points,GRAY)]:
    ax.step([x for x,y in points],[y for x,y in points],where="post",lw=2.2,color=color,label=label)
    ax.text(12.05,points[-1][1],f"{points[-1][1]:.2f} in",color=color,va="center",fontsize=9)
ax.set(xlim=(0,14),ylim=(-.03,.92),xlabel="Hours since Corona del Mar first increment",
       ylabel="New measured rain (inches)",title="November 1 local rain")
ax.grid(axis="y",color=LIGHT);ax.legend(frameon=False,loc="upper left",fontsize=8)
positive=[by_event[name] for name in ["2025-11-15","2025-11-21","2025-12-24","2026-11-01"]]
for y,row in enumerate(positive[::-1]):
    cdm=dt.datetime.fromisoformat(row["coronaDelMarFirstNearbyIncrementUTC"])
    peak=dt.datetime.fromisoformat(row["newportPeakCenterUTC"])
    sna=dt.datetime.fromisoformat(row["snaFirstMeasurableHourEndUTC"])
    clock.scatter([(peak-cdm).total_seconds()/3600],[y],marker="D",s=74,color=TEAL,zorder=3)
    clock.scatter([(sna-cdm).total_seconds()/3600],[y],marker="o",s=70,facecolor="white",edgecolor=ORANGE,lw=2,zorder=3)
clock.axvline(0,color=NAVY,lw=1)
clock.set(xlim=(-1,27),ylim=(-.5,3.6),xlabel="Hours since Corona del Mar first increment",
          title="First airport report versus water peak")
clock.set_yticks(range(4),["Nov 1, 2026","Dec 24, 2025","Nov 21, 2025","Nov 15, 2025"])
clock.grid(axis="x",color=LIGHT)
clock.legend(handles=[plt.Line2D([],[],marker="D",ls="",color=TEAL,label="Newport one-hour peak"),
                      plt.Line2D([],[],marker="o",ls="",markerfacecolor="white",markeredgecolor=ORANGE,color=ORANGE,label="First measurable SNA hour end")],
             frameon=False,loc="upper left",fontsize=8)
fig.suptitle("The local rain gauges change the storm clock",y=1.04)
fig.text(.05,-.04,"SNA reports hourly rain ending at its time stamp. CDMC1 and SANC1 are cumulative counters; stale lower values are ignored. Elapsed time handles the clock change.",fontsize=8,color="#36535e")
save(fig,"local-rain-gauges")


# Each storm gets one row: duration and timing are visible without selecting a two-hour maximum.
wet = {name: by_event[name]["coronaAlignedPulse"] for name in ["2025-11-15","2025-11-21","2025-12-24","2026-11-01"]}
labels = ["Nov 15, 2025", "Nov 21, 2025", "Dec 24, 2025", "Nov 1, 2026"]
fig, ax = plt.subplots(figsize=(10.8, 5))
for y, (name, event) in enumerate(list(wet.items())[::-1]):
    onset = event["pulseStartLagHours"]
    end = onset + event["durationAbove2InchesHours"]
    ax.plot([onset, end], [y, y], lw=9, color=TEAL, solid_capstyle="butt", alpha=.75)
    ax.scatter([event["peakLagHours"]], [y], s=85, marker="D", color=ORANGE, zorder=4)
    high = next(x for x in phase["wetEvents"] if x["event"] == name)
    ax.scatter([event["peakLagHours"] + high["hoursFromPeakToNextHigh"]],
               [y], s=75, marker="|", color=NAVY, zorder=5)
    ax.text(end + .3, y, f"{event['durationAbove2InchesHours']:.1f}h", va="center", fontsize=9)
ax.set_yticks(range(4), labels[::-1])
ax.set(xlim=(0, 36), ylim=(-.6, 3.6), xlabel="Hours after first Corona del Mar rain increment",
       title="Four sustained pulses: onset, duration, peak and next high tide")
ax.grid(axis="x", color=LIGHT)
ax.legend(handles=[Patch(facecolor=TEAL, alpha=.75, label="1h mean above +2 in"),
                   plt.Line2D([], [], marker="D", ls="", color=ORANGE, label="1h peak"),
                   plt.Line2D([], [], marker="|", ls="", ms=12, color=NAVY, label="Predicted high after peak")],
          loc="upper left", frameon=False, fontsize=8.5)
fig.text(.1, .01, "Rebased on Corona del Mar rain onset: peak rise +4.6 to +5.9 inches. Six other episodes had no sustained pulse.", fontsize=8, color="#36535e")
save(fig, "four-pulse-timelines")


# Individual multi-day event charts: aligned tide, detided pulse and hourly
# Corona del Mar rain, all on the same local-time axis. No additional sources
# are fetched; the plotted rain totals are checked against comparison.json.
corona_sources = [
    RESEARCH / "local-gauge-check/2025-11-2026-04-cdmc1.csv.gz",
    RESEARCH / "local-gauge-check/2026-10-cdmc1.csv.gz",
    RESEARCH / "local-gauge-check/2026-11-01-cdmc1.csv.gz",
]
corona_records = {}
for path in corona_sources:
    with gzip.open(path, "rt") as stream:
        for row in csv.DictReader(stream):
            if row["station"] != "CDMC1" or not row["PCIRRZZ"]:
                continue
            when = dt.datetime.strptime(row["utc_valid"], "%Y-%m-%d %H:%M:%S").replace(tzinfo=dt.timezone.utc)
            value = float(row["PCIRRZZ"])
            old = corona_records.get(when)
            if old is not None and old != value:
                raise ValueError(f"Conflicting Corona del Mar counter at {when}")
            corona_records[when] = value
corona_increments = []
previous_time = running_max = None
for when, value in sorted(corona_records.items()):
    if previous_time is not None and when - previous_time > dt.timedelta(days=30):
        running_max = value
    if running_max is None:
        running_max = value
    if value > running_max + .005:
        corona_increments.append((when, round(value - running_max, 3)))
        running_max = value
    previous_time = when

timeline_metadata = []
rain_tide_overlap = []
def predicted_tide(when):
    stamp = int(when.timestamp() * 1000)
    if stamp >= water.study.recent_stamps[0]:
        return water.study.interpolated(water.study.recent_noaa, water.study.recent_stamps, stamp, 360_000)
    return water.study.interpolated(water.study.hourly_noaa, water.study.hourly_stamps, stamp, water.H)


def incoming_tide(when):
    before = predicted_tide(when - dt.timedelta(hours=1))
    after = predicted_tide(when + dt.timedelta(hours=1))
    return before is not None and after is not None and after > before


for name, label in zip(
    ["2025-11-15", "2025-11-21", "2025-12-24", "2026-11-01"],
    ["November 15, 2025", "November 21, 2025", "December 24, 2025", "November 1, 2026"],
):
    event = by_event[name]
    onset = dt.datetime.fromisoformat(event["coronaDelMarFirstNearbyIncrementUTC"])
    start = onset - dt.timedelta(hours=18)
    end = onset + dt.timedelta(hours=42)
    onset_ms = int(onset.timestamp() * 1000)
    pulse = water.longest_pulse(onset_ms, event["coronaAlignedPulse"]["pre12hResidualMeanInches"])
    pulse_start = dt.datetime.fromisoformat(pulse["approximateCenterStartUTC"])
    pulse_end = dt.datetime.fromisoformat(pulse["approximateCenterEndUTC"])
    peak = dt.datetime.fromisoformat(event["newportPeakCenterUTC"])
    high = dt.datetime.fromisoformat(next(p for p in phase["wetEvents"] if p["event"] == name)["nextPredictedHighUTC"])
    assert start < onset < pulse_start < peak < pulse_end < end

    rain = [(t, v) for t, v in corona_increments if onset <= t < end]
    rain6 = round(sum(v for t, v in rain if t < onset + dt.timedelta(hours=6)), 3)
    rain24 = round(sum(v for t, v in rain if t < onset + dt.timedelta(hours=24)), 3)
    assert rain6 == event["coronaSinceOnset6hInches"], (name, rain6)
    assert rain24 == event["coronaSinceOnset24hInches"], (name, rain24)
    hours = [sum(v for t, v in rain if k <= (t-onset).total_seconds()/3600 < k+1) for k in range(42)]
    rain_times = [onset + dt.timedelta(hours=k) for k in range(42)]
    rain_last24 = max((t for t, _ in rain if t < onset + dt.timedelta(hours=24)), default=onset)
    before_peak = [(t, v) for t, v in rain if t <= peak]
    incoming_before_peak = [(t, v) for t, v in before_peak if incoming_tide(t)]
    assert incoming_before_peak, name
    rain_tide_overlap.append({
        "event": name,
        "firstRainOnIncomingTide": incoming_tide(onset),
        "rainBeforePeakInches": round(sum(v for _, v in before_peak), 3),
        "rainOnIncomingTideBeforePeakInches": round(sum(v for _, v in incoming_before_peak), 3),
        "firstIncomingRainUTC": incoming_before_peak[0][0].isoformat(),
        "lastRainBeforePeakUTC": before_peak[-1][0].isoformat(),
        "pulseEndHoursAfterPredictedHigh": round((pulse_end - high).total_seconds() / 3600, 2),
    })

    grid_start = int(start.timestamp() * 1000) // water.STEP * water.STEP
    grid_end = int(end.timestamp() * 1000) // water.STEP * water.STEP
    tide_times, observed, predicted = [], [], []
    for stamp in range(grid_start, grid_end + 1, water.STEP):
        residual = water.study.residual(stamp)
        if residual is None:
            continue
        tide_times.append(dt.datetime.fromtimestamp(stamp / 1000, dt.timezone.utc))
        observed.append(water.study.uc[stamp])
        predicted.append(water.study.uc[stamp] - residual / 12)
    assert len(tide_times) > 600, (name, len(tide_times))
    rolling_times, rolling_rise = [], []
    baseline = event["coronaAlignedPulse"]["pre12hResidualMeanInches"]
    for stamp in range(grid_start, grid_end - water.H + 1, water.STEP):
        window = water.interval(stamp, 1)
        if window["coverage"] >= .9:
            rolling_times.append(dt.datetime.fromtimestamp((stamp + water.H//2) / 1000, dt.timezone.utc))
            rolling_rise.append(window["meanInches"] - baseline)
    assert rolling_times

    fig, (tide_ax, pulse_ax, rain_ax) = plt.subplots(
        3, 1, figsize=(13.2, 8.6), sharex=True,
        gridspec_kw={"height_ratios": [2.25, 1.4, 1.15], "hspace": .11},
    )
    fig.subplots_adjust(top=.84, bottom=.13, left=.1, right=.93)
    fig.suptitle(f"{label}: rain, tide and the later water-level pulse", y=.97, fontsize=16, weight="bold")
    onset_label = onset.astimezone(PACIFIC).strftime("%b %-d, %-I:%M %p %Z")
    fig.text(.1, .91, f"First Corona del Mar rain: {onset_label}   |   Rain: {rain6:.2f} in first 6h, {rain24:.2f} in first 24h   |   Peak: +{event['coronaAlignedPulse']['peakOneHourRiseInches']:.1f} in", fontsize=10, color=NAVY)

    tide_ax.plot(tide_times, predicted, color=GRAY, lw=1.4, ls="--", label="NOAA predicted tide")
    tide_ax.plot(tide_times, observed, color=NAVY, lw=1.7, label="UC Davis measured water")
    tide_ax.axvspan(pulse_start, pulse_end, color=TEAL, alpha=.13, label="Sustained pulse")
    tide_ax.scatter([peak], [water.study.uc[int(peak.timestamp()*1000)]], marker="D", s=60,
                    color=TEAL, edgecolor="white", lw=.8, zorder=5)
    tide_ax.set_ylabel("Water level (ft MLLW)")
    tide_ax.set_title("Measured water and predicted astronomical tide", loc="left", fontsize=10, pad=8)
    tide_ax.legend(frameon=False, loc="upper left", ncol=3, fontsize=8.5)
    tide_ax.grid(axis="y", color=LIGHT)

    pulse_ax.plot(rolling_times, rolling_rise, color=TEAL, lw=2)
    pulse_ax.axhline(0, color=GRAY, lw=.8)
    pulse_ax.axhline(2, color=TEAL, lw=.85, ls="--")
    pulse_ax.axvspan(pulse_start, pulse_end, color=TEAL, alpha=.15)
    peak_rise = event["coronaAlignedPulse"]["peakOneHourRiseInches"]
    pulse_ax.scatter([peak], [peak_rise], marker="D", s=64, color=TEAL, edgecolor="white", lw=.8, zorder=5)
    pulse_ax.text(.99, .96, f"Diamond: pulse peak +{peak_rise:.1f} in\n{event['coronaAlignedPulse']['peakLagHours']:.1f}h after first rain",
                  transform=pulse_ax.transAxes, ha="right", va="top", fontsize=8.5, color=TEAL,
                  bbox={"facecolor": "white", "edgecolor": "none", "alpha": .88, "pad": 3})
    pulse_ax.text(.01, .96, f"Shaded: +2 in run, {pulse['durationHours']:.1f}h", transform=pulse_ax.transAxes,
                  va="top", fontsize=8.5, color=TEAL)
    pulse_ax.set_ylabel("1h residual rise\nabove pre-rain 12h mean (in)")
    pulse_ax.grid(axis="y", color=LIGHT)

    rain_ax.bar(rain_times, hours, width=dt.timedelta(hours=1), align="edge", color=ORANGE, edgecolor="white", linewidth=.25)
    rain_ax.axvspan(onset, rain_last24, color=ORANGE, alpha=.08)
    rain_ax.scatter([onset], [0], marker="v", s=70, color=ORANGE, zorder=5, clip_on=False)
    rain_ax.set_ylabel("Corona del Mar\nhourly rain (in)")
    rain_ax.set_ylim(0, max(.12, max(hours) * 1.3))
    rain_ax.text(.01, .95, f"Orange bars: measured rain; {rain24:.2f} in in first 24h", transform=rain_ax.transAxes,
                 va="top", fontsize=8.5, color=ORANGE)
    rain_ax.grid(axis="y", color=LIGHT)

    for axis in (tide_ax, pulse_ax, rain_ax):
        axis.axvline(onset, color=ORANGE, lw=1.1, ls=":")
        axis.axvline(high, color=GRAY, lw=1, ls=":")
        axis.set_xlim(start, end)
    rain_ax.xaxis.set_major_locator(mdates.HourLocator(byhour=[0, 6, 12, 18], tz=PACIFIC))
    rain_ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %-d\n%-I %p %Z", tz=PACIFIC))
    rain_ax.tick_params(axis="x", labelsize=8)
    rain_ax.set_xlabel("Newport local time; 60-hour chart around first rain")
    fig.text(.1, .045, "Orange line/triangle: first rain. Gray dotted line: predicted high after the marked pulse peak. Residual = measured water − predicted tide; one-hour mean is compared with the preceding 12-hour mean. Times are Pacific.", fontsize=8, color="#36535e")
    stem = f"event-tide-rain-pulse-{name}"
    save(fig, stem)
    timeline_metadata.append({
        "chart": stem, "event": name, "rainOnsetUTC": onset.isoformat(),
        "rainFirst6hInches": rain6, "rainFirst24hInches": rain24,
        "pulseStartUTC": pulse_start.isoformat(), "pulsePeakUTC": peak.isoformat(),
        "pulseEndUTC": pulse_end.isoformat(), "pulsePeakRiseInches": peak_rise,
        "predictedHighUTC": high.isoformat(), "waterPairsPlotted": len(tide_times),
    })


# The completed 24-hour total is useful for retrospective comparison; some
# storms continued raining after the water peak, so it is not a forecast input.
fig, ax = plt.subplots(figsize=(9.6, 5.2))
for event in [r for r in comparison["events"] if r["coronaAlignedPulse"]]:
    rain24 = event["coronaSinceOnset24hInches"]
    rise = event["coronaAlignedPulse"]["peakOneHourRiseInches"]
    ax.scatter(rain24, rise, s=90,
               facecolor=ORANGE if event["sustainedPulse"] else "white",
               edgecolor=ORANGE if event["sustainedPulse"] else GRAY, lw=1.8, zorder=3)
    if event["event"] in ("2026-11-01", "2025-12-26"):
        label = "Nov 1, 2026" if event["event"] == "2026-11-01" else "Dec 26, 2025"
        ax.annotate(f"{label}: +{rise:.1f} in", (rain24, rise), xytext=(12, 5),
                    textcoords="offset points", fontsize=9, color=NAVY)
ax.axhline(4, color=ORANGE, ls="--", lw=.9, alpha=.75)
ax.set(xlim=(-.06, 2.05), ylim=(0, 6.9),
       xlabel="Corona del Mar rain in first 24 hours (in)",
       ylabel="Peak one-hour residual rise (in)",
       title="Rain amount alone does not select the pulse")
ax.grid(axis="y", color=LIGHT)
ax.legend(handles=[plt.Line2D([], [], marker="o", ls="", ms=8, markerfacecolor=ORANGE,
                              markeredgecolor=ORANGE, label="Strong pulse"),
                   plt.Line2D([], [], marker="o", ls="", ms=8, markerfacecolor="white",
                              markeredgecolor=GRAY, label="Other rain episode")],
          frameon=False, loc="upper left", fontsize=9)
save(fig, "rain-amount-screen")


# Keep the conditional Friday-rain scenario separate from the light Saturday model run.
def local(day, hour=0, minute=0):
    return dt.datetime(2026, 11, day, hour, minute, tzinfo=PACIFIC)


fig, ax = plt.subplots(figsize=(10.8, 4.4))
fig.subplots_adjust(left=.28, right=.97, top=.70, bottom=.22)
fig.suptitle("Saturday's high comes before the model's light rain", x=.04, y=.96,
             ha="left", fontsize=15, weight="bold")
fig.text(.04, .84, "NWS Nov 6: late-week rain uncertain; generally light if any",
         fontsize=11, color="#36535e")
watch_start, watch_end = local(14, 8, 50), local(14, 9, 50)
model_start, model_end = local(14, 11), local(14, 16)
high = local(14, 10, 52)
ax.broken_barh([(mdates.date2num(watch_start), mdates.date2num(watch_end) - mdates.date2num(watch_start))],
               (2.05, .35), facecolors=ORANGE)
ax.text(local(14, 9, 20), 2.53, "8:50–9:50 a.m.", ha="center", fontsize=12, weight="bold", color=ORANGE)
ax.scatter([high], [1.35], marker="^", color=NAVY, s=120, zorder=4)
ax.text(high, 1.68, "10:52 a.m. · 4.89 ft", ha="center", fontsize=12, weight="bold")
ax.broken_barh([(mdates.date2num(model_start), mdates.date2num(model_end) - mdates.date2num(model_start))],
               (.25, .35), facecolors="#b5cace")
ax.text(local(14, 13, 30), .72, "11 a.m.–4 p.m. · ~0.19 in coast", ha="center", fontsize=12)
ax.set(xlim=(local(14, 7), local(14, 17)), ylim=(-.05, 2.85),
       xlabel="Saturday, November 14 · Newport local time")
ax.xaxis.label.set_size(12)
ax.set_yticks([2.22, 1.35, .42],
              ["If rain starts Friday evening:\nconditional peak watch",
               "NOAA predicted high", "One model: light rain"])
ax.tick_params(axis="y", length=0, pad=12, labelsize=12)
ax.set_xticks([local(14, hour) for hour in (8, 10, 12, 14, 16)],
              ["8 a.m.", "10 a.m.", "noon", "2 p.m.", "4 p.m."])
ax.tick_params(axis="x", labelsize=11)
ax.grid(axis="x", color=LIGHT)
save(fig, "forecast-tide-calendar")

provenance = {
    "purpose": "Draft blog figures only; no website build or publication",
    "sources": {str(path.relative_to(ROOT)): hashlib.sha256(path.read_bytes()).hexdigest()
                for path in [*sources.values(),comparison_path,*corona_sources,sanc_path,sna_path,
                             ROOT / "data/research/kelvin-improve-v4-20261024/audit/water-600-days.json.gz",
                             ROOT / "data/research/data-page-refresh-20261105/live-uc-davis.gz",
                             ROOT / "data/research/kelvin-continued-20261103/tidal-prediction-features/noaa-9410580-hourly-20251104-20261104.json.gz",
                             ROOT / "data/research/data-page-refresh-20261105/live-noaa.gz",
                             ROOT / "public/data/tides/daily-highs.json"]},
    "charts": ["rain-across-tide-turn", "november-1-pulse", "local-rain-gauges", "four-pulse-timelines", "rain-amount-screen", "forecast-tide-calendar",
               *[row["chart"] for row in timeline_metadata]],
    "eventTimelines": timeline_metadata,
    "rainTideOverlap": rain_tide_overlap,
    "rainTideOverlapMethod": "Corona del Mar positive cumulative-gauge increments from onset through pulse peak. Incoming tide means NOAA predicted water is higher one hour after the increment than one hour before it; pulse end uses the sustained moving-one-hour rise above +2 inches. Rain across tide turn is conceptual and not a causal attribution.",
    "notes": "NWS outlook is as of 2026-11-06 4:27 a.m. PST. Forecast calendar uses a retained ECMWF-based run, not a consensus forecast. All times Pacific, accounting for DST. Corona del Mar onset and rain amount sensitivity are in the comparison artifact; complete rain/tide sequences are in local-gauge-check/rain-tide-sequence.json.",
}
(OUT / "provenance.json").write_text(json.dumps(provenance, indent=2) + "\n")
print("Wrote 10 chart pairs and provenance in", OUT)