add Manning's n calibration and SRC comparison plot scripts

.gitignore

@@ -0,0 +1,42 @@
+# Python
+__pycache__/
+*.py[cod]
+*.pyo
+*.pyd
+.Python
+*.egg-info/
+*.egg
+dist/
+build/
+
+# Virtual environment
+.venv/
+venv/
+env/
+
+# Pytest / coverage
+.pytest_cache/
+.coverage
+htmlcov/
+
+# Data files (large, transferred separately)
+data/*.parquet
+data/*.xlsx
+data/*.csv
+data/*.tif
+data/*.geojson
+
+# Pipeline cache and intermediate outputs
+cache/
+
+# Jupyter
+.ipynb_checkpoints/
+*.ipynb
+
+# OS
+.DS_Store
+Thumbs.db
+
+# IDE
+.vscode/settings.json
+.idea/

scripts/README.md

@@ -0,0 +1,59 @@
+# FIMbox pipeline scripts
+
+Run everything with `.venv\Scripts\python.exe scripts/<name>.py` **from the
+repo root** (several scripts resolve `data/` and other paths relative to the
+working directory).
+
+## 1. Pipeline run order
+
+| # | Script | What it does |
+|---|---|---|
+| 0 | `s00_study_boundary.py` | Builds the study-area boundary from `data/study_area.xlsx` (run once) |
+| 1 | `s01_download_dem_fast.py` | Downloads the 10 m 3DEP DEM for the study area |
+| 2 | `s02_stage1_all_hucs.py` | Stage 1 — data download, per HUC8 |
+| 3 | `s03_stage2_all_hucs.py` | Stage 2 — HAND processing, produces SRCs + hydroTables per branch |
+| 3b | `s03b_fix_nonmonotonic_srcs.py` | **Run right after s03.** Forces every rating curve monotone (raster artifacts from HAND geometry can otherwise make discharge dip with rising stage). Fixed ~3,000 curves study-wide in the last full run. |
+| 4 | `s04_usgs_gage_crosswalk_all_hucs.py` | Crosswalks USGS gauges onto branches/HydroIDs; writes `usgs_elev_table.csv` per HUC |
+| 5 | `s05d_calibrate_B3_wedge_continuity.py` | **Manning's n calibration — the production method.** See §2. |
+| 6 | `s06_stage3_all_hucs.py` | Streamflow / feature-ID extraction for a forecast event |
+| 7 | `s07_stage4_all_hucs.py` | **Ablation baseline only** — OWP's own channel/overbank calibration. Overwrites the B3 hydroTables. Do not run this in the main sequence; only run it (on a restored baseline) when producing the L0/L1/L2 comparison arms against B3. |
+| 8 | `s08_stage5_all_hucs.py` | FIM generation — inundation depth/extent rasters |
+
+## 2. Calibration — how n is calibrated and why B3 was chosen
+
+- **`s05_ncalib_core.py`** — the shared engine. Not run directly; everything
+  below imports it. Handles the starting-point match (datum/bathymetry
+  correction), conveyance smoothing, and both n-application schemes
+  (whole-section and incremental).
+- **`s05d_calibrate_B3_wedge_continuity.py`** — **the method actually applied
+  in production.** Continuity-wedge starting-point match + incremental
+  (slice) n application. Writes both n definitions to every calibrated
+  hydroTable: `zonal_n_applied`/`n_{yr}yr` (slice n — drives `discharge_cms`,
+  what FIM consumes) and `whole_section_n`/`n_eff_{yr}yr` (an independently
+  calibrated whole-section curve, `discharge_cms_wholesection`, kept for
+  comparison with OWP's 0.06/0.12 defaults — not used by FIM).
+- **`s05a`/`s05b`/`s05c`** (methods A, B1, B2) — the three *rejected*
+  starting-point corrections. Kept only so the method comparison is
+  reproducible — **do not use for production runs.**
+- **`s05e_evaluate_methods.py`** — reruns A/B1/B2/B3 offline (reads baseline
+  backups, writes nothing to the pipeline) and produces the evidence for
+  choosing B3. Output: `E:/SI/out/calibration_analysis/method_evaluation/`.
+- **`s05f_evaluate_schemes.py`** — same idea for whole-section vs incremental
+  n application; justifies the incremental choice. Output:
+  `E:/SI/out/calibration_analysis/scheme_evaluation/`.
+- Full narrative + evidence: `docs/method_decision_deck.pptx`.
+
+## 3. Analysis & export
+
+| Script | Produces | Output folder |
+|---|---|---|
+| `plot_v2_diagnostics.py` | Datum-gap, clip-rate, n-profile, anchor-residual, hold-out validation, and starting-point figures + one RC panel per calibrated gauge | `E:/SI/out/calibration_analysis/v2_diagnostics/` |
+| `plot_calibration_analysis.py` | Study-area aggregate plots (n distributions, Q scatter, coverage, metrics) + per-branch RC plots | `E:/SI/out/calibration_analysis/` and `E:/SI/out/HUC{huc8}/src_plots/` |
+| `plot_src_ensemble.py` | SRC ensemble comparisons (side-by-side, overlaid, normalized) across all calibrated gauges | `E:/SI/out/calibration_analysis/` |
+| `export_rc_panels_and_videos.py` | Consolidates all per-gauge RC panels into one folder in three annotation variants (slice n / whole-section n / both), plus a slideshow video per variant | `E:/SI/out/calibration_analysis/rc_panels_export/` |
+
+## `legacy/`
+
+Retired scripts from earlier phases of the project (single-HUC prototypes,
+the old monolithic calibration script, earlier plotting tools) — kept for
+history, not maintained. See `legacy/README.md` for what superseded each one.

scripts/export_rc_panels_and_videos.py

@@ -0,0 +1,197 @@
+"""
+Consolidate per-gauge RC panels into one export folder and render them as
+slideshow videos — three variants of the SAME calibrated curve, differing
+only in which n definition is annotated in the title:
+
+    slice/          annotated with slice n (applied, incremental scheme)
+    wholesection/    annotated with whole-section n (effective, back-computed)
+    both/            annotated with both definitions stacked
+
+Outputs (one consolidated location):
+    E:/SI/out/calibration_analysis/rc_panels_export/
+        slice/{NNN}_{huc8}_{loc}.png
+        wholesection/{NNN}_{huc8}_{loc}.png
+        both/{NNN}_{huc8}_{loc}.png
+        video_slice_n.mp4
+        video_wholesection_n.mp4
+        video_both_n.mp4
+
+Reuses the data-assembly logic from plot_v2_diagnostics.py (same anchors,
+same baseline/calibrated hydroTable curves) — only the panel's annotation
+and output path differ. Panels are numbered in a fixed (huc8, location_id)
+sort order so the three videos show gauges in the same sequence.
+
+Run AFTER plot_v2_diagnostics.py (needs ncalib_diagnostics.csv per HUC):
+    .venv\\Scripts\\python.exe scripts/export_rc_panels_and_videos.py
+"""
+from __future__ import annotations
+
+import logging
+from pathlib import Path
+
+import matplotlib as mpl
+mpl.use("Agg")   # headless backend — TkAgg default is ~10-20x slower for batch savefig
+
+import imageio.v2 as imageio
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+import plot_v2_diagnostics as v2  # reuse load_shared / load_diagnostics / gauge_records
+
+mpl.rcParams.update({
+    "font.family": "DejaVu Sans", "font.size": 10.5,
+    "axes.labelsize": 12, "axes.titlesize": 10.5, "axes.titleweight": "bold",
+    "axes.spines.top": False, "axes.spines.right": False,
+    "axes.grid": True, "grid.alpha": 0.20, "grid.linestyle": ":",
+    "legend.fontsize": 9.5, "legend.framealpha": 0.92, "figure.dpi": 140,
+})
+
+C_ORIG, C_CALIB, C_OBS = "#2471a3", "#c0392b", "#111111"
+C_WSEC = "#7f8c8d"   # whole-section curve — same gray as the scheme evaluation
+RY = v2.RECURRENCE_YEARS
+
+EXPORT_DIR = Path("E:/SI/out/calibration_analysis/rc_panels_export")
+VARIANTS = ("slice", "wholesection", "both")
+FPS = 10
+HOLD_SECONDS = 1.5
+
+
+def n_line(d: pd.Series, col_prefix: str) -> str:
+    return "   ".join(f"n({yr}) = {d[f'{col_prefix}_{yr}yr']:.3f}"
+                      for yr in RY if pd.notna(d.get(f"{col_prefix}_{yr}yr")))
+
+
+def draw_panel(g: dict, anchors: list, s_o: pd.DataFrame, s_c: pd.DataFrame,
+               grc: pd.DataFrame, variant: str) -> plt.Figure:
+    fig, ax = plt.subplots(figsize=(9.5, 6.6), constrained_layout=True)
+
+    h_max = max(a[1] for a in anchors) * 1.6
+    q_max = max(a[2] for a in anchors) * 1.8
+
+    for i, (yr, h_r, _) in enumerate(anchors):
+        ax.axhspan(0 if i == 0 else anchors[i - 1][1], h_r,
+                   color="#187E86", alpha=0.045 if i % 2 == 0 else 0.09)
+    if not grc.empty:
+        ax.plot(grc["flow_cms"], grc["hand"], color=C_OBS, lw=2.0,
+                label="USGS observed RC", zorder=5)
+    ax.plot(s_o["Discharge (m3s-1)"], s_o["Stage"], color=C_ORIG, lw=1.8,
+            label="Baseline SRC (n = 0.06)", zorder=4)
+
+    # The two calibrated curves are genuinely different objects: the slice
+    # (incremental) curve is what FIM consumes; the whole-section curve is an
+    # independently calibrated diagnostic (discharge_cms_wholesection column).
+    has_wsec = ("discharge_cms_wholesection" in s_c.columns
+                and s_c["discharge_cms_wholesection"].notna().any())
+    if variant in ("slice", "both"):
+        ax.plot(s_c["Discharge (m3s-1)"], s_c["Stage"], color=C_CALIB, lw=2.6,
+                ls="--", label="Calibrated SRC — slice n (applied)", zorder=6)
+    if variant in ("wholesection", "both") and has_wsec:
+        ax.plot(s_c["discharge_cms_wholesection"], s_c["Stage"], color=C_WSEC,
+                lw=2.2, ls="-.", label="Calibrated SRC — whole-section n",
+                zorder=5)
+    elif variant == "wholesection" and not has_wsec:
+        # single-anchor / all-pinned gauges: the two schemes coincide
+        ax.plot(s_c["Discharge (m3s-1)"], s_c["Stage"], color=C_WSEC, lw=2.2,
+                ls="-.", label="Calibrated SRC — whole-section n (= slice)",
+                zorder=5)
+    for yr, h_r, Q_r in anchors:
+        ax.scatter([Q_r], [h_r], s=52, color=C_CALIB, edgecolors="white",
+                   linewidths=1.4, zorder=8)
+        ax.annotate(f"{yr}-yr", (Q_r, h_r), textcoords="offset points",
+                    xytext=(7, -3), fontsize=8.5, color=C_CALIB, fontweight="bold")
+
+    d = g["diag"]
+    mode_txt = {"wedge_continuity": f"continuity wedge  δ = {d['delta_m']:.2f} m",
+                "wedge_topwidth":   f"TopWidth wedge  δ = {d['delta_m']:.2f} m",
+                "q0_offset":        f"baseflow offset  Q₀ = {d['q0_cms']:.1f} m³/s"}.get(
+                    d["mode"], d["mode"])
+
+    if variant == "slice":
+        n_txt = "Slice n (applied):  " + n_line(d, "n")
+    elif variant == "wholesection":
+        n_txt = "Whole-section n (effective):  " + n_line(d, "n_eff")
+    else:  # both
+        n_txt = ("Slice n:          " + n_line(d, "n") + "\n" +
+                 "Whole-section n:  " + n_line(d, "n_eff"))
+
+    ax.set_title(f"Gauge {g['location_id']}  ·  HydroID {g['hydroid']}  ·  "
+                 f"HUC {g['huc8']}\n{mode_txt}\n{n_txt}", fontsize=9.5)
+    ax.set_xlim(0, q_max)
+    ax.set_ylim(0, h_max)
+    ax.set_xlabel("Discharge  (m³/s)")
+    ax.set_ylabel("HAND stage  (m)")
+    ax.legend(loc="lower right")
+    return fig
+
+
+def main():
+    for v in VARIANTS:
+        (EXPORT_DIR / v).mkdir(parents=True, exist_ok=True)
+
+    hucs = [str(int(c)).zfill(8)
+            for c in pd.read_excel(v2.EXCEL_PATH)[v2.HUC_CODE_COL]]
+    rc, recur = v2.load_shared()
+    diag = v2.load_diagnostics(hucs)
+    records = sorted(v2.gauge_records(diag, rc, recur),
+                     key=lambda r: (r[0]["huc8"], r[0]["location_id"]))
+    log.info("Rendering %d gauges × %d variants", len(records), len(VARIANTS))
+
+    paths = {v: [] for v in VARIANTS}
+    for i, (g, anchors, s_o, s_c, grc) in enumerate(records, 1):
+        fname = f"{i:03d}_{g['huc8']}_{g['location_id']}.png"
+        for variant in VARIANTS:
+            fig = draw_panel(g, anchors, s_o, s_c, grc, variant)
+            out = EXPORT_DIR / variant / fname
+            fig.savefig(out, dpi=140, bbox_inches="tight", facecolor="white")
+            plt.close(fig)
+            paths[variant].append(out)
+        if i % 20 == 0:
+            log.info("  %d/%d gauges rendered", i, len(records))
+    log.info("Rendered %d gauges → %s", len(records), EXPORT_DIR)
+
+    labels = {"slice": "Slice n (applied, incremental scheme)",
+              "wholesection": "Whole-section n (effective)",
+              "both": "Both n definitions"}
+    for variant in VARIANTS:
+        make_video(paths[variant], EXPORT_DIR / f"video_{variant}_n.mp4",
+                  labels[variant])
+
+
+def make_video(png_paths: list[Path], out_path: Path, label: str) -> None:
+    if not png_paths:
+        log.warning("No frames for %s — skipping video", out_path.name)
+        return
+    frames = [imageio.imread(p) for p in png_paths]
+    h, w = frames[0].shape[:2]
+    # matplotlib bbox_inches='tight' can vary output pixel size slightly per
+    # panel — pad/crop every frame to a common canvas so the video encoder
+    # gets a constant frame size.
+    canvas_h = max(f.shape[0] for f in frames)
+    canvas_w = max(f.shape[1] for f in frames)
+    canvas_h += canvas_h % 2   # even dimensions for the video codec
+    canvas_w += canvas_w % 2
+
+    hold_frames = max(1, round(FPS * HOLD_SECONDS))
+    writer = imageio.get_writer(str(out_path), fps=FPS, codec="libx264",
+                                quality=8, macro_block_size=None)
+    for f in frames:
+        padded = np.full((canvas_h, canvas_w, 3), 255, dtype=np.uint8)
+        fh, fw = f.shape[:2]
+        rgb = f[..., :3] if f.shape[-1] == 4 else f
+        padded[:fh, :fw] = rgb
+        for _ in range(hold_frames):
+            writer.append_data(padded)
+    writer.close()
+    log.info("Saved video (%s, %d gauges, %.0fs): %s",
+             label, len(png_paths), len(png_paths) * HOLD_SECONDS, out_path)
+
+
+logging.basicConfig(level=logging.INFO,
+                    format="%(asctime)s [%(levelname)s] %(message)s",
+                    datefmt="%H:%M:%S")
+log = logging.getLogger(__name__)
+
+
+if __name__ == "__main__":
+    main()

scripts/legacy/README.md

@@ -0,0 +1,26 @@
+# Legacy scripts
+
+These files are retired — kept for history only, not maintained, and not
+part of the active pipeline. Each is superseded by a file in `scripts/`
+(see `../README.md` for the active set):
+
+| Legacy file | Superseded by |
+|---|---|
+| `s01_download_dem.py` | `s01_download_dem_fast.py` (reads AOI window from cloud-optimized GeoTIFFs directly — ~8x faster, no national VRT parse) |
+| `s02_stage1_single_huc.py` | `s02_stage1_all_hucs.py` |
+| `s03_stage2_single_huc.py` | `s03_stage2_all_hucs.py` |
+| `s04_usgs_gage_crosswalk.py` | `s04_usgs_gage_crosswalk_all_hucs.py` |
+| `s05_calibrate_n_recurrence.py` | `s05_ncalib_core.py` + `s05a`–`s05d` |
+| `s05_revert_uncalibrated_srcs.py` | Not needed — fixed a failure mode (trunk-branch averaging propagating onto ungauged reaches) that cannot occur under the new pipeline, which only ever touches gauged HydroIDs. Was hardcoded to one HUC (`03020102`) and the pre-correction `D:/SI/out` drive path. |
+| `s06_stage3_single_huc.py` | `s06_stage3_all_hucs.py` |
+| `s07_stage4_single_huc.py` | `s07_stage4_all_hucs.py` |
+| `s08_stage5_single_huc.py` | `s08_stage5_all_hucs.py` |
+| `plot_calibration_slides.py` | The method-decision deck (`docs/method_decision_deck.pptx`) and `plot_v2_diagnostics.py` |
+| `plot_src_calibration.py`, `plot_src_calibration_inset.py`, `plot_srcs.py` | `plot_v2_diagnostics.py`, `export_rc_panels_and_videos.py` |
+| `plot_workflow_chart.py` | — (superseded by the pipeline order documented in `../README.md`) |
+| `open_parquet.py` | Scratch one-liner; use `pandas.read_parquet` directly |
+
+`scripts/s05_calibrate_n_recurrence_all_hucs.py` (the old monolithic v2
+calibration script) is not moved here — it has local uncommitted edits that
+are intentionally being kept out of git history. It is also fully superseded
+by `s05_ncalib_core.py` + `s05a`–`s05d`.

scripts/legacy/open_parquet.py

@@ -0,0 +1,4 @@
+# impport parquet
+import pyarrow.parquet as pq
+table = pq.read_table(r"C:\Users\Ali\OneDrive - CUNY\Desktop\SI\fimbox_SI26\data\usgs_rating_curves.parquet")
+print(table)