#!/usr/bin/env python3
"""
Station-by-station regression analysis for NOAA and STOFS NAVD water levels.

Input is a single table with at least:
    station_id, time, noaa_navd_m, stofs_navd_m

Outputs are written under OUTPUT_DIR:
    figures/raw/{station_id}_raw_regression.png
    figures/detrended/{station_id}_detrended_regression.png
    station_regression_summary.csv
    processing_log.csv
    station_regression_report.html
"""

from __future__ import annotations

from dataclasses import dataclass
from html import escape
from pathlib import Path
import sys
import traceback

# Optional project-local dependencies, used only when this directory exists.
LOCAL_DEPS = Path(__file__).resolve().parent / ".python_deps"
if LOCAL_DEPS.exists():
    sys.path.insert(0, str(LOCAL_DEPS))

import matplotlib

matplotlib.use("Agg")
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pylib import ReadNC, get_stat
import statsmodels.api as sm


# ---------------------------------------------------------------------------
# User settings
# ---------------------------------------------------------------------------

ROOT = Path("/sciclone/schism10/sgao08")

# Change this to the merged hourly CSV containing the columns below.
INPUT_FILE = ROOT / "STOFS_single/regression/noaa_stofs_hourly_navd.csv"

# If INPUT_FILE does not exist, the script can use the existing paired NetCDF
# hourly products in this workspace and compute the same station-month table.
USE_NETCDF_IF_CSV_MISSING = True
NOAA_NC_FILE = ROOT / "STOFS_single/v1.1/obs_navd_hourly_2000_2025.nc"
STOFS_NC_FILE = ROOT / "STOFS_single/v1.1/stofs_navd_hourly_2000_2025.nc"
NETCDF_WATER_VAR = "water_level_navd_m"
NETCDF_STATION_ID_VAR = "station_id"
NETCDF_TIME_VAR = "time"

OUTPUT_DIR = ROOT / "STOFS_single/regression"
RAW_FIG_DIR = OUTPUT_DIR / "figures/raw"
DETRENDED_FIG_DIR = OUTPUT_DIR / "figures/detrended"
SUMMARY_CSV = OUTPUT_DIR / "station_regression_summary.csv"
PROCESSING_LOG_CSV = OUTPUT_DIR / "processing_log.csv"
REPORT_HTML = OUTPUT_DIR / "station_regression_report.html"

STATION_COL = "station_id"
TIME_COL = "time"
NOAA_COL = "noaa_navd_m"
STOFS_COL = "stofs_navd_m"

MIN_COMMON_MONTHS = 120
FIG_DPI = 200
TREND_HAC_MAXLAGS = 12

# Require adequate coverage from both NOAA and STOFS for each monthly mean.
MIN_VALID_HOURS_PER_MONTH = 300


@dataclass
class TrendResult:
    slope_m_per_year: float
    slope_mm_per_year: float
    intercept_m: float
    slope_ci_low_mm_per_year: float
    slope_ci_high_mm_per_year: float
    p_value: float
    r_squared: float
    fitted: np.ndarray


@dataclass
class PairRegressionResult:
    slope: float
    intercept: float
    p_value: float
    r_squared: float
    fitted: np.ndarray


def ensure_output_dirs() -> None:
    OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
    RAW_FIG_DIR.mkdir(parents=True, exist_ok=True)
    DETRENDED_FIG_DIR.mkdir(parents=True, exist_ok=True)


def safe_station_id(value: object) -> str:
    text = str(value).strip()
    return "".join(ch if ch.isalnum() or ch in ("-", "_", ".") else "_" for ch in text)


def read_input(input_file: Path) -> pd.DataFrame:
    if not input_file.exists():
        raise FileNotFoundError(
            f"Input file not found: {input_file}. Edit INPUT_FILE at the top of this script."
        )

    df = pd.read_csv(input_file)
    required = {STATION_COL, TIME_COL, NOAA_COL, STOFS_COL}
    missing = sorted(required.difference(df.columns))
    if missing:
        raise ValueError(f"Input file is missing required columns: {missing}")

    df = df[[STATION_COL, TIME_COL, NOAA_COL, STOFS_COL]].copy()
    df[STATION_COL] = df[STATION_COL].astype(str).str.strip()
    df[TIME_COL] = pd.to_datetime(df[TIME_COL], errors="coerce")
    df[NOAA_COL] = pd.to_numeric(df[NOAA_COL], errors="coerce")
    df[STOFS_COL] = pd.to_numeric(df[STOFS_COL], errors="coerce")

    df = df[df[STATION_COL].ne("") & df[TIME_COL].notna()].copy()
    if df.empty:
        raise ValueError("No rows remain after station/time parsing.")

    # Average duplicate station-time records instead of arbitrarily keeping one.
    df = (
        df.groupby([STATION_COL, TIME_COL], as_index=False)[[NOAA_COL, STOFS_COL]]
        .mean()
        .sort_values([STATION_COL, TIME_COL])
    )
    return df


def hourly_to_common_monthly(df: pd.DataFrame) -> pd.DataFrame:
    work = df.copy()
    work["month"] = work[TIME_COL].dt.to_period("M").dt.to_timestamp()

    monthly = (
        work.groupby([STATION_COL, "month"])
        .agg(
            noaa_monthly_m=(NOAA_COL, "mean"),
            stofs_monthly_m=(STOFS_COL, "mean"),
            noaa_valid_hours=(NOAA_COL, "count"),
            stofs_valid_hours=(STOFS_COL, "count"),
        )
        .reset_index()
    )

    monthly = monthly[
        monthly["noaa_monthly_m"].notna()
        & monthly["stofs_monthly_m"].notna()
        & (monthly["noaa_valid_hours"] >= MIN_VALID_HOURS_PER_MONTH)
        & (monthly["stofs_valid_hours"] >= MIN_VALID_HOURS_PER_MONTH)
    ].copy()

    return monthly.sort_values([STATION_COL, "month"])


def netcdf_time_to_datetime(time_var) -> pd.DatetimeIndex:
    units = getattr(time_var, "units", "")
    values = np.asarray(time_var[:], dtype=float)
    if units.startswith("seconds since "):
        origin_text = units.replace("seconds since ", "", 1).replace("+0000", "").strip()
        origin = pd.to_datetime(origin_text, utc=True).tz_localize(None)
        return pd.DatetimeIndex(origin + pd.to_timedelta(values, unit="s"))

    from netCDF4 import num2date

    dates = num2date(values, units=units)
    return pd.DatetimeIndex(pd.to_datetime([date.isoformat() for date in dates]))


def decode_netcdf_station_ids(station_id_var) -> list[str]:
    from netCDF4 import chartostring

    values = station_id_var[:]
    if values.dtype.kind in {"S", "U"} and values.ndim == 2:
        return [str(item).strip() for item in chartostring(values)]
    return [str(item).strip() for item in values]


def paired_netcdf_to_common_monthly(noaa_nc: Path, stofs_nc: Path) -> tuple[pd.DataFrame, list[str]]:
    if not noaa_nc.exists() or not stofs_nc.exists():
        raise FileNotFoundError(
            f"CSV input is missing and paired NetCDF files are not available: {noaa_nc}, {stofs_nc}"
        )

    rows = []
    noaa_ds = ReadNC(str(noaa_nc), fmt=1)
    stofs_ds = ReadNC(str(stofs_nc), fmt=1)
    try:
        for var_name in [NETCDF_TIME_VAR, NETCDF_STATION_ID_VAR, NETCDF_WATER_VAR]:
            if var_name not in noaa_ds.variables:
                raise ValueError(f"Missing {var_name} in {noaa_nc}")
            if var_name not in stofs_ds.variables:
                raise ValueError(f"Missing {var_name} in {stofs_nc}")

        noaa_times = netcdf_time_to_datetime(noaa_ds.variables[NETCDF_TIME_VAR])
        stofs_times = netcdf_time_to_datetime(stofs_ds.variables[NETCDF_TIME_VAR])
        if len(noaa_times) != len(stofs_times) or not np.array_equal(noaa_times.values, stofs_times.values):
            raise ValueError("NOAA and STOFS NetCDF time axes do not match.")

        noaa_ids = decode_netcdf_station_ids(noaa_ds.variables[NETCDF_STATION_ID_VAR])
        stofs_ids = decode_netcdf_station_ids(stofs_ds.variables[NETCDF_STATION_ID_VAR])
        stofs_index = {station_id: idx for idx, station_id in enumerate(stofs_ids)}
        station_ids = [station_id for station_id in noaa_ids if station_id in stofs_index]

        if not station_ids:
            raise ValueError("No common station_id values between paired NetCDF files.")

        month_index = noaa_times.to_period("M").to_timestamp()
        base = pd.DataFrame({"month": month_index})
        noaa_var = noaa_ds.variables[NETCDF_WATER_VAR]
        stofs_var = stofs_ds.variables[NETCDF_WATER_VAR]

        for noaa_idx, station_id in enumerate(noaa_ids):
            if station_id not in stofs_index:
                continue
            stofs_idx = stofs_index[station_id]
            noaa_values = np.ma.filled(noaa_var[:, noaa_idx], np.nan).astype(float)
            stofs_values = np.ma.filled(stofs_var[:, stofs_idx], np.nan).astype(float)
            station_frame = base.copy()
            station_frame[NOAA_COL] = noaa_values
            station_frame[STOFS_COL] = stofs_values
            monthly = (
                station_frame.groupby("month")
                .agg(
                    noaa_monthly_m=(NOAA_COL, "mean"),
                    stofs_monthly_m=(STOFS_COL, "mean"),
                    noaa_valid_hours=(NOAA_COL, "count"),
                    stofs_valid_hours=(STOFS_COL, "count"),
                )
                .reset_index()
            )
            monthly = monthly[
                monthly["noaa_monthly_m"].notna()
                & monthly["stofs_monthly_m"].notna()
                & (monthly["noaa_valid_hours"] >= MIN_VALID_HOURS_PER_MONTH)
                & (monthly["stofs_valid_hours"] >= MIN_VALID_HOURS_PER_MONTH)
            ].copy()
            monthly[STATION_COL] = station_id
            rows.append(monthly)
    finally:
        noaa_ds.close()
        stofs_ds.close()

    if not rows:
        return pd.DataFrame(columns=[STATION_COL, "month", "noaa_monthly_m", "stofs_monthly_m"]), station_ids

    monthly_all = pd.concat(rows, ignore_index=True)
    columns = [STATION_COL, "month", "noaa_monthly_m", "stofs_monthly_m", "noaa_valid_hours", "stofs_valid_hours"]
    return monthly_all[columns].sort_values([STATION_COL, "month"]), station_ids


def years_since_start(dates: pd.Series | pd.DatetimeIndex) -> np.ndarray:
    dt = pd.to_datetime(dates)
    start = dt.min()
    return ((dt - start) / pd.Timedelta(days=365.2425)).to_numpy(dtype=float)


def fit_time_trend(dates: pd.Series, values: pd.Series) -> TrendResult:
    y = pd.to_numeric(values, errors="coerce").to_numpy(dtype=float)
    x = years_since_start(dates)
    mask = np.isfinite(x) & np.isfinite(y)
    if mask.sum() < 3:
        raise ValueError("Need at least 3 valid points for time trend.")

    x_valid = x[mask]
    y_valid = y[mask]
    model = sm.OLS(y_valid, sm.add_constant(x_valid)).fit(
        cov_type="HAC",
        cov_kwds={"maxlags": TREND_HAC_MAXLAGS},
    )
    fitted = np.full_like(y, np.nan, dtype=float)
    fitted[mask] = model.predict(sm.add_constant(x_valid))

    ci = model.conf_int(alpha=0.05)
    slope = float(model.params[1])
    return TrendResult(
        slope_m_per_year=slope,
        slope_mm_per_year=slope * 1000.0,
        intercept_m=float(model.params[0]),
        slope_ci_low_mm_per_year=float(ci[1, 0] * 1000.0),
        slope_ci_high_mm_per_year=float(ci[1, 1] * 1000.0),
        p_value=float(model.pvalues[1]),
        r_squared=float(model.rsquared),
        fitted=fitted,
    )


def fit_pair_regression(x_values: pd.Series | np.ndarray, y_values: pd.Series | np.ndarray) -> PairRegressionResult:
    x = pd.to_numeric(pd.Series(x_values), errors="coerce").to_numpy(dtype=float)
    y = pd.to_numeric(pd.Series(y_values), errors="coerce").to_numpy(dtype=float)
    mask = np.isfinite(x) & np.isfinite(y)
    if mask.sum() < 3:
        raise ValueError("Need at least 3 valid points for pair regression.")

    x_valid = x[mask]
    y_valid = y[mask]
    model = sm.OLS(y_valid, sm.add_constant(x_valid)).fit()
    fitted = np.full_like(y, np.nan, dtype=float)
    fitted[mask] = model.predict(sm.add_constant(x_valid))

    return PairRegressionResult(
        slope=float(model.params[1]),
        intercept=float(model.params[0]),
        p_value=float(model.pvalues[1]),
        r_squared=float(model.rsquared),
        fitted=fitted,
    )


def pair_metrics(noaa: pd.Series | np.ndarray, stofs: pd.Series | np.ndarray) -> dict[str, float]:
    obs = pd.to_numeric(pd.Series(noaa), errors="coerce").to_numpy(dtype=float)
    mod = pd.to_numeric(pd.Series(stofs), errors="coerce").to_numpy(dtype=float)
    mask = np.isfinite(obs) & np.isfinite(mod)
    if mask.sum() < 2:
        return {
            "pearson_r": np.nan,
            "bias_m": np.nan,
            "rmse_m": np.nan,
            "mae_m": np.nan,
        }

    metrics = get_stat(mod[mask], obs[mask])
    return {
        "pearson_r": float(metrics.R),
        "bias_m": float(metrics.ME),
        "rmse_m": float(metrics.RMSD),
        "mae_m": float(metrics.MAE),
    }


def format_float(value: float, digits: int = 3) -> str:
    if value is None or not np.isfinite(value):
        return "NA"
    return f"{value:.{digits}f}"


def format_p(value: float) -> str:
    if value is None or not np.isfinite(value):
        return "NA"
    if value < 1e-4:
        return f"{value:.2e}"
    return f"{value:.4f}"


def set_month_axis(ax: plt.Axes) -> None:
    ax.xaxis.set_major_locator(mdates.YearLocator(base=2))
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
    ax.grid(True, alpha=0.25)


def equal_limits(a: np.ndarray, b: np.ndarray, pad_fraction: float = 0.06) -> tuple[float, float]:
    values = np.concatenate([np.asarray(a, dtype=float), np.asarray(b, dtype=float)])
    values = values[np.isfinite(values)]
    if values.size == 0:
        return -1.0, 1.0
    lo = float(np.min(values))
    hi = float(np.max(values))
    if np.isclose(lo, hi):
        pad = max(abs(lo) * 0.1, 0.1)
    else:
        pad = (hi - lo) * pad_fraction
    return lo - pad, hi + pad


def stats_text_raw(station_id: str, monthly: pd.DataFrame, metrics: dict[str, float]) -> str:
    lines = [
        f"station_id: {station_id}",
        f"period: {monthly['month'].min():%Y-%m} to {monthly['month'].max():%Y-%m}",
        f"common months: {len(monthly)}",
        "",
        "NOAA trend:",
        (
            f"  {format_float(metrics['raw_noaa_trend_mm_per_year'])} mm/yr "
            f"[{format_float(metrics['raw_noaa_trend_ci_low_mm_per_year'])}, "
            f"{format_float(metrics['raw_noaa_trend_ci_high_mm_per_year'])}], "
            f"p={format_p(metrics['raw_noaa_trend_p_value'])}, "
            f"R2={format_float(metrics['raw_noaa_trend_r_squared'])}"
        ),
        "STOFS trend:",
        (
            f"  {format_float(metrics['raw_stofs_trend_mm_per_year'])} mm/yr "
            f"[{format_float(metrics['raw_stofs_trend_ci_low_mm_per_year'])}, "
            f"{format_float(metrics['raw_stofs_trend_ci_high_mm_per_year'])}], "
            f"p={format_p(metrics['raw_stofs_trend_p_value'])}, "
            f"R2={format_float(metrics['raw_stofs_trend_r_squared'])}"
        ),
        (
            "difference trend: "
            f"{format_float(metrics['raw_difference_trend_mm_per_year'])} mm/yr, "
            f"p={format_p(metrics['raw_difference_trend_p_value'])}"
        ),
        "",
        (
            "STOFS = intercept + slope x NOAA: "
            f"slope={format_float(metrics['raw_pair_slope'])}, "
            f"intercept={format_float(metrics['raw_pair_intercept_m'])} m, "
            f"R2={format_float(metrics['raw_pair_r_squared'])}"
        ),
        (
            f"Pearson r={format_float(metrics['raw_pearson_r'])}, "
            f"bias={format_float(metrics['raw_bias_m'])} m, "
            f"RMSE={format_float(metrics['raw_rmse_m'])} m, "
            f"MAE={format_float(metrics['raw_mae_m'])} m"
        ),
    ]
    return "\n".join(lines)


def stats_text_detrended(station_id: str, monthly: pd.DataFrame, metrics: dict[str, float]) -> str:
    lines = [
        f"station_id: {station_id}",
        f"common months: {len(monthly)}",
        "",
        f"NOAA detrended anomaly SD: {format_float(metrics['detrended_noaa_sd_m'])} m",
        f"STOFS detrended anomaly SD: {format_float(metrics['detrended_stofs_sd_m'])} m",
        (
            "regression: "
            f"slope={format_float(metrics['detrended_pair_slope'])}, "
            f"intercept={format_float(metrics['detrended_pair_intercept_m'])} m"
        ),
        f"Pearson r: {format_float(metrics['detrended_pearson_r'])}",
        f"R2: {format_float(metrics['detrended_pair_r_squared'])}",
        f"RMSE: {format_float(metrics['detrended_rmse_m'])} m",
        f"MAE: {format_float(metrics['detrended_mae_m'])} m",
    ]
    return "\n".join(lines)


def plot_raw_station(
    station_id: str,
    monthly: pd.DataFrame,
    noaa_trend: TrendResult,
    stofs_trend: TrendResult,
    diff_trend: TrendResult,
    pair_reg: PairRegressionResult,
    metrics: dict[str, float],
    out_png: Path,
) -> None:
    dates = monthly["month"]
    noaa = monthly["noaa_monthly_m"].to_numpy(dtype=float)
    stofs = monthly["stofs_monthly_m"].to_numpy(dtype=float)
    diff = stofs - noaa

    fig, axes = plt.subplots(2, 2, figsize=(13, 9), constrained_layout=True)
    fig.suptitle(f"{station_id} Raw Monthly Mean Regression", fontsize=15)

    ax = axes[0, 0]
    ax.plot(dates, noaa, color="#1f77b4", lw=1.1, label="NOAA monthly mean")
    ax.plot(dates, stofs, color="#d62728", lw=1.1, label="STOFS monthly mean")
    ax.plot(dates, noaa_trend.fitted, color="#1f77b4", lw=2.0, ls="--", label="NOAA trend")
    ax.plot(dates, stofs_trend.fitted, color="#d62728", lw=2.0, ls="--", label="STOFS trend")
    ax.set_title("A. Monthly Mean Time Series")
    ax.set_ylabel("Water level (m NAVD88)")
    ax.legend(fontsize=8)
    set_month_axis(ax)

    ax = axes[0, 1]
    ax.plot(dates, diff, color="#2ca02c", lw=1.1, label="STOFS - NOAA")
    ax.plot(dates, diff_trend.fitted, color="black", lw=2.0, ls="--", label="difference trend")
    ax.axhline(0.0, color="0.45", lw=0.8)
    ax.set_title("B. Difference Time Series")
    ax.set_ylabel("STOFS - NOAA (m)")
    ax.legend(fontsize=8)
    set_month_axis(ax)

    ax = axes[1, 0]
    ax.scatter(noaa, stofs, s=15, alpha=0.65, color="#4c78a8", edgecolors="none", label="common months")
    lo, hi = equal_limits(noaa, stofs)
    ax.plot([lo, hi], [lo, hi], color="0.35", lw=1.2, label="1:1 line")
    order = np.argsort(noaa)
    ax.plot(noaa[order], pair_reg.fitted[order], color="#d62728", lw=2.0, label="regression")
    ax.set_xlim(lo, hi)
    ax.set_ylim(lo, hi)
    ax.set_aspect("equal", adjustable="box")
    ax.set_title("C. STOFS vs NOAA")
    ax.set_xlabel("NOAA monthly mean (m NAVD88)")
    ax.set_ylabel("STOFS monthly mean (m NAVD88)")
    ax.grid(True, alpha=0.25)
    ax.legend(fontsize=8)

    ax = axes[1, 1]
    ax.axis("off")
    ax.set_title("D. Station Statistics")
    ax.text(
        0.02,
        0.98,
        stats_text_raw(station_id, monthly, metrics),
        va="top",
        ha="left",
        family="monospace",
        fontsize=8.5,
        transform=ax.transAxes,
    )

    fig.savefig(out_png, dpi=FIG_DPI)
    plt.close(fig)


def compute_detrended_anomalies(monthly: pd.DataFrame) -> tuple[pd.DataFrame, TrendResult, TrendResult]:
    det = monthly.copy()
    det["calendar_month"] = det["month"].dt.month

    noaa_clim = det.groupby("calendar_month")["noaa_monthly_m"].transform("mean")
    stofs_clim = det.groupby("calendar_month")["stofs_monthly_m"].transform("mean")
    det["noaa_anomaly_m"] = det["noaa_monthly_m"] - noaa_clim
    det["stofs_anomaly_m"] = det["stofs_monthly_m"] - stofs_clim

    noaa_anom_trend = fit_time_trend(det["month"], det["noaa_anomaly_m"])
    stofs_anom_trend = fit_time_trend(det["month"], det["stofs_anomaly_m"])
    det["noaa_detrended_anomaly_m"] = det["noaa_anomaly_m"] - noaa_anom_trend.fitted
    det["stofs_detrended_anomaly_m"] = det["stofs_anomaly_m"] - stofs_anom_trend.fitted
    return det, noaa_anom_trend, stofs_anom_trend


def plot_detrended_station(
    station_id: str,
    det: pd.DataFrame,
    pair_reg: PairRegressionResult,
    metrics: dict[str, float],
    out_png: Path,
) -> None:
    dates = det["month"]
    noaa = det["noaa_detrended_anomaly_m"].to_numpy(dtype=float)
    stofs = det["stofs_detrended_anomaly_m"].to_numpy(dtype=float)
    diff = stofs - noaa
    diff_mean = float(np.nanmean(diff))

    fig, axes = plt.subplots(2, 2, figsize=(13, 9), constrained_layout=True)
    fig.suptitle(f"{station_id} Detrended Monthly Anomaly Regression", fontsize=15)

    ax = axes[0, 0]
    ax.plot(dates, noaa, color="#1f77b4", lw=1.0, label="NOAA detrended anomaly")
    ax.plot(dates, stofs, color="#d62728", lw=1.0, label="STOFS detrended anomaly")
    ax.axhline(0.0, color="0.45", lw=0.8)
    ax.set_title("A. Detrended Monthly Anomalies")
    ax.set_ylabel("Detrended anomaly (m)")
    ax.legend(fontsize=8)
    set_month_axis(ax)

    ax = axes[0, 1]
    ax.plot(dates, diff, color="#2ca02c", lw=1.0, label="detrended STOFS - NOAA")
    ax.axhline(diff_mean, color="black", lw=1.8, ls="--", label=f"mean={diff_mean:.3f} m")
    ax.axhline(0.0, color="0.45", lw=0.8)
    ax.set_title("B. Detrended Difference")
    ax.set_ylabel("Difference (m)")
    ax.legend(fontsize=8)
    set_month_axis(ax)

    ax = axes[1, 0]
    ax.scatter(noaa, stofs, s=15, alpha=0.65, color="#4c78a8", edgecolors="none", label="common months")
    lo, hi = equal_limits(noaa, stofs)
    ax.plot([lo, hi], [lo, hi], color="0.35", lw=1.2, label="1:1 line")
    order = np.argsort(noaa)
    ax.plot(noaa[order], pair_reg.fitted[order], color="#d62728", lw=2.0, label="regression")
    ax.set_xlim(lo, hi)
    ax.set_ylim(lo, hi)
    ax.set_aspect("equal", adjustable="box")
    ax.set_title("C. Detrended STOFS vs NOAA")
    ax.set_xlabel("NOAA detrended anomaly (m)")
    ax.set_ylabel("STOFS detrended anomaly (m)")
    ax.grid(True, alpha=0.25)
    ax.legend(fontsize=8)

    ax = axes[1, 1]
    ax.axis("off")
    ax.set_title("D. Detrended Statistics")
    ax.text(
        0.02,
        0.98,
        stats_text_detrended(station_id, det, metrics),
        va="top",
        ha="left",
        family="monospace",
        fontsize=9,
        transform=ax.transAxes,
    )

    fig.savefig(out_png, dpi=FIG_DPI)
    plt.close(fig)


def analyze_station(station_id: str, monthly: pd.DataFrame) -> dict[str, object]:
    monthly = monthly.sort_values("month").reset_index(drop=True)
    station_safe = safe_station_id(station_id)
    raw_png = RAW_FIG_DIR / f"{station_safe}_raw_regression.png"
    det_png = DETRENDED_FIG_DIR / f"{station_safe}_detrended_regression.png"

    noaa = monthly["noaa_monthly_m"]
    stofs = monthly["stofs_monthly_m"]
    diff = stofs - noaa

    noaa_trend = fit_time_trend(monthly["month"], noaa)
    stofs_trend = fit_time_trend(monthly["month"], stofs)
    diff_trend = fit_time_trend(monthly["month"], diff)
    raw_pair_reg = fit_pair_regression(noaa, stofs)
    raw_pair_metrics = pair_metrics(noaa, stofs)

    metrics: dict[str, object] = {
        "station_id": station_id,
        "first_common_month": monthly["month"].min().strftime("%Y-%m"),
        "last_common_month": monthly["month"].max().strftime("%Y-%m"),
        "common_months": int(len(monthly)),
        "raw_figure": raw_png.relative_to(OUTPUT_DIR).as_posix(),
        "detrended_figure": det_png.relative_to(OUTPUT_DIR).as_posix(),
        "raw_noaa_trend_mm_per_year": noaa_trend.slope_mm_per_year,
        "raw_noaa_trend_ci_low_mm_per_year": noaa_trend.slope_ci_low_mm_per_year,
        "raw_noaa_trend_ci_high_mm_per_year": noaa_trend.slope_ci_high_mm_per_year,
        "raw_noaa_trend_p_value": noaa_trend.p_value,
        "raw_noaa_trend_r_squared": noaa_trend.r_squared,
        "raw_stofs_trend_mm_per_year": stofs_trend.slope_mm_per_year,
        "raw_stofs_trend_ci_low_mm_per_year": stofs_trend.slope_ci_low_mm_per_year,
        "raw_stofs_trend_ci_high_mm_per_year": stofs_trend.slope_ci_high_mm_per_year,
        "raw_stofs_trend_p_value": stofs_trend.p_value,
        "raw_stofs_trend_r_squared": stofs_trend.r_squared,
        "raw_difference_trend_mm_per_year": diff_trend.slope_mm_per_year,
        "raw_difference_trend_ci_low_mm_per_year": diff_trend.slope_ci_low_mm_per_year,
        "raw_difference_trend_ci_high_mm_per_year": diff_trend.slope_ci_high_mm_per_year,
        "raw_difference_trend_p_value": diff_trend.p_value,
        "raw_difference_trend_r_squared": diff_trend.r_squared,
        "raw_pair_slope": raw_pair_reg.slope,
        "raw_pair_intercept_m": raw_pair_reg.intercept,
        "raw_pair_p_value": raw_pair_reg.p_value,
        "raw_pair_r_squared": raw_pair_reg.r_squared,
        "raw_pearson_r": raw_pair_metrics["pearson_r"],
        "raw_bias_m": raw_pair_metrics["bias_m"],
        "raw_rmse_m": raw_pair_metrics["rmse_m"],
        "raw_mae_m": raw_pair_metrics["mae_m"],
    }

    plot_raw_station(
        station_id,
        monthly,
        noaa_trend,
        stofs_trend,
        diff_trend,
        raw_pair_reg,
        metrics,
        raw_png,
    )

    det, _, _ = compute_detrended_anomalies(monthly)
    det_noaa = det["noaa_detrended_anomaly_m"]
    det_stofs = det["stofs_detrended_anomaly_m"]
    det_pair_reg = fit_pair_regression(det_noaa, det_stofs)
    det_pair_metrics = pair_metrics(det_noaa, det_stofs)

    metrics.update(
        {
            "detrended_noaa_sd_m": float(np.nanstd(det_noaa, ddof=1)),
            "detrended_stofs_sd_m": float(np.nanstd(det_stofs, ddof=1)),
            "detrended_pair_slope": det_pair_reg.slope,
            "detrended_pair_intercept_m": det_pair_reg.intercept,
            "detrended_pair_p_value": det_pair_reg.p_value,
            "detrended_pair_r_squared": det_pair_reg.r_squared,
            "detrended_pearson_r": det_pair_metrics["pearson_r"],
            "detrended_rmse_m": det_pair_metrics["rmse_m"],
            "detrended_mae_m": det_pair_metrics["mae_m"],
        }
    )

    plot_detrended_station(
        station_id,
        det,
        det_pair_reg,
        metrics,
        det_png,
    )
    return metrics


def write_processing_log(log_rows: list[dict[str, object]]) -> None:
    columns = ["station_id", "status", "reason", "common_months", "first_month", "last_month"]
    pd.DataFrame(log_rows).reindex(columns=columns).to_csv(PROCESSING_LOG_CSV, index=False)


def write_summary(summary_rows: list[dict[str, object]]) -> pd.DataFrame:
    summary = pd.DataFrame(summary_rows)
    if not summary.empty:
        summary = summary.sort_values("station_id")
    summary.to_csv(SUMMARY_CSV, index=False)
    return summary


def html_number(value: object, digits: int = 3, is_p: bool = False) -> str:
    try:
        val = float(value)
    except (TypeError, ValueError):
        return ""
    if not np.isfinite(val):
        return ""
    return escape(format_p(val) if is_p else format_float(val, digits))


def html_sort_value(value: object) -> str:
    try:
        val = float(value)
    except (TypeError, ValueError):
        return ""
    if not np.isfinite(val):
        return "NaN"
    return f"{val:.12g}"


def html_num_cell(value: object, digits: int = 3, is_p: bool = False) -> str:
    return f"<td data-sort='{html_sort_value(value)}'>{html_number(value, digits=digits, is_p=is_p)}</td>"


def trend_label(slope: object, p_value: object) -> tuple[str, str]:
    try:
        slope_val = float(slope)
        p_val = float(p_value)
    except (TypeError, ValueError):
        return "not significant", "not-significant"
    if not np.isfinite(slope_val) or not np.isfinite(p_val):
        return "not significant", "not-significant"
    if p_val >= 0.05:
        return "not significant", "not-significant"
    if slope_val > 0:
        return "significant increase", "significant-increase"
    return "significant decrease", "significant-decrease"


def write_html_report(summary: pd.DataFrame, log_rows: list[dict[str, object]]) -> None:
    ok_count = len(summary)
    skipped_count = sum(1 for row in log_rows if row.get("status") == "skipped")
    failed_count = sum(1 for row in log_rows if row.get("status") == "failed")

    rows = []
    for _, row in summary.iterrows():
        raw_link = escape(str(row.get("raw_figure", "")))
        det_link = escape(str(row.get("detrended_figure", "")))
        label, label_class = trend_label(
            row.get("raw_difference_trend_mm_per_year"),
            row.get("raw_difference_trend_p_value"),
        )
        rows.append(
            "<tr>"
            f"<td>{escape(str(row.get('station_id', '')))}</td>"
            f"<td>{escape(str(row.get('first_common_month', '')))}</td>"
            f"<td>{escape(str(row.get('last_common_month', '')))}</td>"
            f"{html_num_cell(row.get('common_months'), digits=0)}"
            f"{html_num_cell(row.get('raw_noaa_trend_mm_per_year'))}"
            f"{html_num_cell(row.get('raw_stofs_trend_mm_per_year'))}"
            f"{html_num_cell(row.get('raw_difference_trend_mm_per_year'))}"
            f"{html_num_cell(row.get('raw_difference_trend_p_value'), is_p=True)}"
            f"<td><span class='tag {label_class}'>{escape(label)}</span></td>"
            f"{html_num_cell(row.get('raw_pair_slope'))}"
            f"{html_num_cell(row.get('raw_pair_r_squared'))}"
            f"{html_num_cell(row.get('raw_pearson_r'))}"
            f"{html_num_cell(row.get('raw_bias_m'))}"
            f"{html_num_cell(row.get('raw_rmse_m'))}"
            f"{html_num_cell(row.get('detrended_pair_slope'))}"
            f"{html_num_cell(row.get('detrended_pair_r_squared'))}"
            f"{html_num_cell(row.get('detrended_pearson_r'))}"
            f"{html_num_cell(row.get('detrended_rmse_m'))}"
            f"<td><a href='{raw_link}'><img src='{raw_link}' alt='raw regression'></a></td>"
            f"<td><a href='{det_link}'><img src='{det_link}' alt='detrended regression'></a></td>"
            "</tr>"
        )

    log_table_rows = []
    for log_row in log_rows:
        status = str(log_row.get("status", ""))
        status_class = {
            "success": "success",
            "skipped": "insufficient",
            "failed": "failed",
        }.get(status, "not-significant")
        log_table_rows.append(
            "<tr>"
            f"<td>{escape(str(log_row.get('station_id', '')))}</td>"
            f"<td><span class='tag {status_class}'>{escape(status)}</span></td>"
            f"<td>{escape(str(log_row.get('reason', '')))}</td>"
            f"{html_num_cell(log_row.get('common_months'), digits=0)}"
            f"<td>{escape(str(log_row.get('first_month', '')))}</td>"
            f"<td>{escape(str(log_row.get('last_month', '')))}</td>"
            "</tr>"
        )
    log_html = "".join(log_table_rows)

    if summary.empty:
        period_text = "No successful station rows"
    else:
        period_text = (
            f"{escape(str(summary['first_common_month'].min()))} to "
            f"{escape(str(summary['last_common_month'].max()))} common monthly NOAA/STOFS NAVD water levels"
        )

    html = f"""<!doctype html>
<html lang="en">
<head>
  <meta charset="utf-8">
  <meta name="viewport" content="width=device-width, initial-scale=1">
  <title>NOAA-STOFS Station Regression Report</title>
  <style>
    body {{ font-family: Arial, sans-serif; margin: 18px; color: #20242a; background: #f7f8f4; }}
    h1 {{ font-size: 22px; margin: 0 0 8px; }}
    h2 {{ font-size: 17px; margin: 18px 0 8px; }}
    .meta {{ color: #555; margin-bottom: 12px; line-height: 1.4; }}
    input {{ width: min(520px, 95vw); padding: 8px 10px; border: 1px solid #9aa3a8; border-radius: 4px; margin: 8px 0 14px; background: white; }}
    table {{ border-collapse: collapse; width: 100%; background: white; font-size: 13px; }}
    th, td {{ border-bottom: 1px solid #d8ddd9; padding: 7px 8px; text-align: left; vertical-align: middle; white-space: nowrap; }}
    th {{ position: sticky; top: 0; background: #e8eee9; cursor: pointer; z-index: 1; }}
    tbody tr:hover {{ background: #f1f6f4; }}
    img {{ width: 168px; height: 110px; object-fit: cover; border: 1px solid #ccd2ce; }}
    a {{ color: #145a7a; }}
    .tag {{ display: inline-block; min-width: 118px; text-align: center; padding: 3px 6px; border-radius: 999px; font-size: 12px; }}
    .significant-increase {{ background: #d7efe0; color: #155c32; }}
    .significant-decrease {{ background: #f3d9d4; color: #8c241c; }}
    .not-significant {{ background: #e6e8eb; color: #3d444b; }}
    .insufficient {{ background: #fff2c7; color: #6a4a00; }}
    .success {{ background: #d7efe0; color: #155c32; }}
    .failed {{ background: #f3d9d4; color: #8c241c; }}
    .table-wrap {{ overflow-x: auto; }}
    .links a {{ margin-right: 12px; }}
  </style>
</head>
<body>
  <h1>NOAA-STOFS Station Regression Report</h1>
  <div class="meta">
    {period_text}; raw monthly mean and detrended anomaly regressions.
    Successful stations: {ok_count}; skipped: {skipped_count}; failed: {failed_count};
    minimum common months: {MIN_COMMON_MONTHS}; minimum valid hours per source-month:
    {MIN_VALID_HOURS_PER_MONTH}.
    <div class="links">
      <a href="{escape(SUMMARY_CSV.name)}">summary CSV</a>
      <a href="{escape(PROCESSING_LOG_CSV.name)}">processing log CSV</a>
    </div>
  </div>
  <input id="search" type="search" placeholder="Search station id, metric value, or trend label">
  <h2>Station Summary</h2>
  <div class="table-wrap">
  <table id="stationTable" class="sortable">
    <thead>
      <tr>
        <th>station_id</th>
        <th>first</th>
        <th>last</th>
        <th>months</th>
        <th>NOAA trend<br>mm/yr</th>
        <th>STOFS trend<br>mm/yr</th>
        <th>diff trend<br>mm/yr</th>
        <th>diff p</th>
        <th>label</th>
        <th>raw slope</th>
        <th>raw R2</th>
        <th>raw r</th>
        <th>raw bias<br>m</th>
        <th>raw RMSE<br>m</th>
        <th>detr slope</th>
        <th>detr R2</th>
        <th>detr r</th>
        <th>detr RMSE<br>m</th>
        <th>raw fig</th>
        <th>detr fig</th>
      </tr>
    </thead>
    <tbody>
      {"".join(rows)}
    </tbody>
  </table>
  </div>
  <h2>Processing Log</h2>
  <div class="table-wrap">
    <table id="logTable" class="sortable">
      <thead>
        <tr>
          <th>station_id</th>
          <th>status</th>
          <th>reason</th>
          <th>common months</th>
          <th>first</th>
          <th>last</th>
        </tr>
      </thead>
      <tbody>
        {log_html}
      </tbody>
    </table>
  </div>
  <script>
    function cellValue(row, index) {{
      const cell = row.children[index];
      if (!cell) return '';
      if (cell.dataset.sort !== undefined) {{
        const numeric = Number(cell.dataset.sort);
        return Number.isFinite(numeric) ? numeric : cell.innerText.toLowerCase();
      }}
      return cell.innerText.toLowerCase();
    }}
    document.querySelectorAll('table.sortable').forEach(table => {{
      const tbody = table.tBodies[0];
      table.querySelectorAll('th').forEach((th, index) => {{
        th.addEventListener('click', () => {{
          const rows = Array.from(tbody.rows);
          const ascending = th.dataset.asc !== 'true';
          table.querySelectorAll('th').forEach(h => h.dataset.asc = '');
          th.dataset.asc = ascending ? 'true' : 'false';
          rows.sort((a, b) => {{
            const av = cellValue(a, index);
            const bv = cellValue(b, index);
            if (Number.isFinite(av) && Number.isFinite(bv)) return ascending ? av - bv : bv - av;
            return ascending ? String(av).localeCompare(String(bv)) : String(bv).localeCompare(String(av));
          }});
          rows.forEach(row => tbody.appendChild(row));
        }});
      }});
    }});
    const search = document.getElementById('search');
    const stationBody = document.getElementById('stationTable').tBodies[0];
    search.addEventListener('input', () => {{
      const q = search.value.toLowerCase();
      Array.from(stationBody.rows).forEach(row => {{
        row.style.display = row.innerText.toLowerCase().includes(q) ? '' : 'none';
      }});
    }});
  </script>
</body>
</html>
"""
    REPORT_HTML.write_text(html, encoding="utf-8")


def run_analysis() -> None:
    ensure_output_dirs()

    try:
        if INPUT_FILE.exists():
            hourly = read_input(INPUT_FILE)
            monthly = hourly_to_common_monthly(hourly)
            all_station_ids = sorted(hourly[STATION_COL].unique())
        elif USE_NETCDF_IF_CSV_MISSING:
            print(f"Input CSV not found, using paired NetCDF files: {NOAA_NC_FILE}, {STOFS_NC_FILE}")
            monthly, all_station_ids = paired_netcdf_to_common_monthly(NOAA_NC_FILE, STOFS_NC_FILE)
        else:
            hourly = read_input(INPUT_FILE)
            monthly = hourly_to_common_monthly(hourly)
            all_station_ids = sorted(hourly[STATION_COL].unique())
    except Exception as exc:
        log_rows = [
            {
                "station_id": "__all__",
                "status": "failed",
                "reason": str(exc),
                "common_months": 0,
                "first_month": "",
                "last_month": "",
            }
        ]
        write_processing_log(log_rows)
        write_summary([])
        write_html_report(pd.DataFrame(), log_rows)
        raise

    summary_rows: list[dict[str, object]] = []
    log_rows: list[dict[str, object]] = []

    for station_id in sorted(all_station_ids):
        station_monthly = monthly[monthly[STATION_COL] == station_id].copy()
        common_months = int(len(station_monthly))
        first_month = station_monthly["month"].min().strftime("%Y-%m") if common_months else ""
        last_month = station_monthly["month"].max().strftime("%Y-%m") if common_months else ""

        if common_months < MIN_COMMON_MONTHS:
            log_rows.append(
                {
                    "station_id": station_id,
                    "status": "skipped",
                    "reason": f"only {common_months} common months; need {MIN_COMMON_MONTHS}",
                    "common_months": common_months,
                    "first_month": first_month,
                    "last_month": last_month,
                }
            )
            continue

        try:
            summary_rows.append(analyze_station(station_id, station_monthly))
            log_rows.append(
                {
                    "station_id": station_id,
                    "status": "success",
                    "reason": "",
                    "common_months": common_months,
                    "first_month": first_month,
                    "last_month": last_month,
                }
            )
        except Exception as exc:
            log_rows.append(
                {
                    "station_id": station_id,
                    "status": "failed",
                    "reason": f"{exc}; {traceback.format_exc(limit=1).strip()}",
                    "common_months": common_months,
                    "first_month": first_month,
                    "last_month": last_month,
                }
            )

    summary = write_summary(summary_rows)
    write_processing_log(log_rows)
    write_html_report(summary, log_rows)

    print(f"Wrote summary: {SUMMARY_CSV}")
    print(f"Wrote processing log: {PROCESSING_LOG_CSV}")
    print(f"Wrote HTML report: {REPORT_HTML}")
    print(f"Wrote figures under: {OUTPUT_DIR / 'figures'}")
    print(f"Successful stations: {len(summary_rows)}")
    print(f"Skipped stations: {sum(1 for row in log_rows if row['status'] == 'skipped')}")
    print(f"Failed stations: {sum(1 for row in log_rows if row['status'] == 'failed')}")


if __name__ == "__main__":
    run_analysis()