diff --git a/6-statistics-fusion-order.py b/6-statistics-fusion-order.py new file mode 100644 index 0000000..3a222cd --- /dev/null +++ b/6-statistics-fusion-order.py @@ -0,0 +1,252 @@ +"""Step 6: Paired ItB-vs-BtI significance test across the full sample. + +Inputs (``data/``, ``{year}`` = ``--evaluation-year``): + +- ``metrics/{year}/{site}/metrics.json`` — per-site validation metrics (Step 5) + +Outputs (``data/statistics_fusion_order/``): + +- ``{year}.json`` — paired Wilcoxon + t-test summary for NSE, RMSE, nRMSE, r + +CLI: + +- ``--evaluation-year`` (default 2025) +- ``--alpha`` (default 0.05) — significance threshold for ``better_order`` + +This step aggregates across all sites with Step 5 output; it does not accept +``--site`` (a single-site filter would not support a sample-level test). +""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path +from typing import Any + +import numpy as np +from scipy.stats import ttest_rel, wilcoxon + +# --------------------------------------------------------------------------- +# Constants +# --------------------------------------------------------------------------- + +DATA_DIR = Path("data") +DEFAULT_YEAR = 2025 +DEFAULT_ALPHA = 0.05 + +METRICS = ["nse", "rmse", "nrmse", "r"] +LOWER_IS_BETTER = {"rmse", "nrmse"} + +MIN_PAIRS_WARNING = 6 + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _r4(v: float | None) -> float | None: + return round(v, 4) if v is not None else None + + +def _load_site_metrics(year: int) -> list[dict[str, Any]]: + """Return parsed ``metrics.json`` payloads for every site under ``{year}``.""" + metrics_dir = DATA_DIR / "metrics" / str(year) + if not metrics_dir.is_dir(): + return [] + + payloads: list[dict[str, Any]] = [] + for site_dir in sorted(metrics_dir.iterdir()): + if not site_dir.is_dir(): + continue + path = site_dir / "metrics.json" + if not path.is_file(): + continue + payloads.append(json.loads(path.read_text())) + return payloads + + +def collect_pairs( + site_metrics: list[dict[str, Any]], metric: str +) -> tuple[list[float], list[float], int]: + """Return paired BtI / ItB values for ``metric`` and count of dropped sites.""" + bti_vals: list[float] = [] + itb_vals: list[float] = [] + n_dropped = 0 + + for payload in site_metrics: + bti = payload.get("bti") + itb = payload.get("itb") + if not isinstance(bti, dict) or not isinstance(itb, dict): + n_dropped += 1 + continue + + bti_v = bti.get(metric) + itb_v = itb.get(metric) + if bti_v is None or itb_v is None: + n_dropped += 1 + continue + + bti_vals.append(float(bti_v)) + itb_vals.append(float(itb_v)) + + return bti_vals, itb_vals, n_dropped + + +def _better_order( + bti_vals: list[float], + itb_vals: list[float], + metric: str, + p_value: float | None, + alpha: float, +) -> str: + """Name the better fusion order when Wilcoxon p < alpha, else no difference.""" + if p_value is None or p_value >= alpha: + return "no significant difference" + + mean_diff = float(np.mean(itb_vals) - np.mean(bti_vals)) + if metric in LOWER_IS_BETTER: + return "itb" if mean_diff < 0 else "bti" + return "itb" if mean_diff > 0 else "bti" + + +def paired_test( + bti_vals: list[float], + itb_vals: list[float], + metric: str, + alpha: float, +) -> dict[str, Any]: + """Run paired Wilcoxon (primary) and t-test; return summary dict.""" + n_pairs = len(bti_vals) + bti_arr = np.array(bti_vals, dtype=float) + itb_arr = np.array(itb_vals, dtype=float) + diffs = itb_arr - bti_arr + + result: dict[str, Any] = { + "n_pairs": n_pairs, + "bti_mean": _r4(float(bti_arr.mean())) if n_pairs else None, + "bti_median": _r4(float(np.median(bti_arr))) if n_pairs else None, + "itb_mean": _r4(float(itb_arr.mean())) if n_pairs else None, + "itb_median": _r4(float(np.median(itb_arr))) if n_pairs else None, + "mean_diff": _r4(float(diffs.mean())) if n_pairs else None, + "median_diff": _r4(float(np.median(diffs))) if n_pairs else None, + "wilcoxon": {"statistic": None, "p_value": None}, + "ttest": {"statistic": None, "p_value": None}, + "better_order": "insufficient data", + } + + if n_pairs < 2: + return result + + wilcoxon_stat: float | None = None + wilcoxon_p: float | None = None + if np.any(diffs != 0): + try: + w_stat, w_p = wilcoxon(itb_arr, bti_arr) + wilcoxon_stat = float(w_stat) + wilcoxon_p = float(w_p) + except ValueError: + pass + + t_stat, t_p = ttest_rel(itb_arr, bti_arr) + + result["wilcoxon"] = { + "statistic": _r4(wilcoxon_stat), + "p_value": _r4(wilcoxon_p), + } + result["ttest"] = { + "statistic": _r4(float(t_stat)), + "p_value": _r4(float(t_p)), + } + result["better_order"] = _better_order( + bti_vals, itb_vals, metric, wilcoxon_p, alpha + ) + return result + + +def _print_summary( + year: int, alpha: float, n_sites_total: int, metrics_out: dict +) -> None: + print(f"\nPaired ItB vs BtI test — {year} (alpha={alpha}, sites={n_sites_total})") + print( + f"{'metric':<8} {'n':>4} {'BtI mean':>10} {'ItB mean':>10} " + f"{'diff':>10} {'W p':>8} {'t p':>8} better" + ) + print("-" * 78) + for metric in METRICS: + m = metrics_out[metric] + bti_mean = m["bti_mean"] + itb_mean = m["itb_mean"] + mean_diff = m["mean_diff"] + w_p = m["wilcoxon"]["p_value"] + t_p = m["ttest"]["p_value"] + better = m["better_order"] + + def _fmt(v: float | None) -> str: + return f"{v:10.4f}" if v is not None else f"{'—':>10}" + + print( + f"{metric:<8} {m['n_pairs']:>4} {_fmt(bti_mean)} {_fmt(itb_mean)} " + f"{_fmt(mean_diff)} " + f"{w_p if w_p is not None else '—':>8} " + f"{t_p if t_p is not None else '—':>8} {better}" + ) + if 0 < m["n_pairs"] < MIN_PAIRS_WARNING: + print( + f" warning: only {m['n_pairs']} pair(s) for {metric}; " + "interpret p-values cautiously" + ) + + +# --------------------------------------------------------------------------- +# CLI +# --------------------------------------------------------------------------- + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--evaluation-year", type=int, default=DEFAULT_YEAR) + parser.add_argument( + "--alpha", + type=float, + default=DEFAULT_ALPHA, + help="Significance threshold for better_order (default 0.05)", + ) + args = parser.parse_args() + + year = args.evaluation_year + alpha = args.alpha + + site_metrics = _load_site_metrics(year) + n_sites_total = len(site_metrics) + if n_sites_total == 0: + raise SystemExit( + f"No Step 5 metrics found under {DATA_DIR / 'metrics' / str(year)}" + ) + + metrics_out: dict[str, Any] = {} + for metric in METRICS: + bti_vals, itb_vals, n_dropped = collect_pairs(site_metrics, metric) + summary = paired_test(bti_vals, itb_vals, metric, alpha) + summary["n_dropped"] = n_dropped + metrics_out[metric] = summary + + payload = { + "year": year, + "alpha": alpha, + "n_sites_total": n_sites_total, + "metrics": metrics_out, + } + + out_dir = DATA_DIR / "statistics_fusion_order" + out_dir.mkdir(parents=True, exist_ok=True) + out_path = out_dir / f"{year}.json" + out_path.write_text(json.dumps(payload, separators=(",", ":"))) + + _print_summary(year, alpha, n_sites_total, metrics_out) + print(f"\nWritten → {out_path}") + + +if __name__ == "__main__": + main() diff --git a/README.md b/README.md index 3d64297..3724bd8 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # EFAST fusion with phenocam validation. -End-to-end pipeline from selecting sites from the global [PhenoCam Network](https://phenocam.nau.edu/) to run [EFAST](https://github.com/DHI-GRAS/efast) spatio-temporal fusion with Sentinel-2 / Sentinel-3 and validate GCCs accross sensors. The five numbered steps cover site selection, Sentinel data acquisition, different fusion orders, and accuracy metrics, all feeding a static web QA viewer. +End-to-end pipeline from selecting sites from the global [PhenoCam Network](https://phenocam.nau.edu/) to run [EFAST](https://github.com/DHI-GRAS/efast) spatio-temporal fusion with Sentinel-2 / Sentinel-3 and validate GCCs accross sensors. The numbered steps cover site selection, Sentinel data acquisition, different fusion orders, accuracy metrics, and sample-level statistics, all feeding a static web QA viewer. --- @@ -13,6 +13,7 @@ End-to-end pipeline from selecting sites from the global [PhenoCam Network](http | 3 | `3-sentinel-data.py` | Acquire S2 (Earth Search COG) and S3 OLCI SYN L2 (CDSE OpenEO); prepare REFL, DIST_CLOUD, and composite GeoTIFFs | | 4 | `4-fusion.py` | Run EFAST BtI (fuse reflectance → GCC) and ItB (fuse GCC directly) for each screened site | | 5 | `5-metrics.py` | Extract PhenoCam-matched timeseries, compute NSE/RMSE/r baselines and fusion metrics, emit per-site JSON and webapp manifest | +| 6 | `6-statistics-fusion-order.py` | Paired ItB-vs-BtI significance test (Wilcoxon + t-test) across all sites | --- @@ -41,9 +42,10 @@ uv run python 2-phenocam-screening.py --evaluation-year 2025 uv run python 3-sentinel-data.py --evaluation-year 2025 uv run python 4-fusion.py --evaluation-year 2025 uv run python 5-metrics.py --evaluation-year 2025 +uv run python 6-statistics-fusion-order.py --evaluation-year 2025 ``` -All steps accept `--evaluation-year` (default `2025`) and `--site` (optional, for single-site runs). Steps 3–5 are resumable — existing output files are skipped. +Steps 1–5 accept `--evaluation-year` (default `2025`) and `--site` (optional, for single-site runs). Step 6 is a full-sample aggregate and only accepts `--evaluation-year` and `--alpha` (default `0.05`). Steps 3–5 are resumable — existing output files are skipped. ```bash # single site @@ -70,6 +72,7 @@ Step 3 S3 download uses CDSE OpenEO (`SENTINEL3_SYN_L2_SYN`). Set `CDSE_USER` an | `fusion/{year}/{site}/bti/`, `.../itb/` | 4 | BtI fused reflectance + GCC; ItB fused GCC | | `metrics/{year}/{site}/` | 5 | Per-site timeseries, metrics, covariates JSON | | `metrics/manifest.json` | 5 | Webapp manifest (years + site metadata) | +| `statistics_fusion_order/{year}.json` | 6 | Paired ItB-vs-BtI test summary (NSE, RMSE, nRMSE, r) | --- diff --git a/index.html b/index.html index a461cbb..644cc54 100644 --- a/index.html +++ b/index.html @@ -147,14 +147,42 @@ body { margin: 0; font: 13px/1.4 system-ui, sans-serif; background: #f5f5f5; col display: flex; align-items: center; gap: 12px; padding: 10px 14px; background: #1a1a2e; color: #eee; flex-shrink: 0; } -#worldHeader h2 { margin: 0; font-size: 14px; font-weight: 600; color: #7eb8f7; flex: 1; } -#worldHeader .world-meta { font-size: 12px; color: #aaa; } +#overlayTabs { display: flex; gap: 4px; flex: 1; } +.otab { + padding: 4px 12px; border-radius: 4px; font-size: 13px; cursor: pointer; + border: 1px solid #555; background: transparent; color: #ccc; +} +.otab.active { background: #2a3f5f; color: #dceeff; border-color: #4a6fa5; } +.otab:hover:not(.active) { background: rgba(255, 255, 255, 0.07); } +#worldHeader .world-meta { font-size: 12px; color: #aaa; flex-shrink: 0; } #worldClose { font-size: 13px; padding: 4px 12px; border-radius: 4px; border: 1px solid #666; background: transparent; color: #ddd; cursor: pointer; + flex-shrink: 0; } #worldClose:hover { background: rgba(255, 255, 255, 0.08); } #worldMap { flex: 1; min-height: 0; } +#statsPanel { flex: 1; min-height: 0; overflow-y: auto; padding: 20px 24px; display: none; background: #f5f5f5; } +.stat-grid { display: grid; grid-template-columns: repeat(auto-fit, minmax(340px, 1fr)); gap: 14px; margin-top: 4px; } +.stat-card { background: #fff; border: 1px solid #e0e0e0; border-radius: 6px; padding: 14px 16px; } +.stat-card h3 { margin: 0 0 10px; font-size: 14px; color: #1a1a2e; display: flex; align-items: center; gap: 8px; flex-wrap: wrap; } +.stat-badge { + font-size: 11px; padding: 2px 7px; border-radius: 10px; font-weight: 600; + background: #e8f5e9; color: #1a6e2e; +} +.stat-badge.itb { background: #fff3e0; color: #c75c00; } +.stat-badge.bti { background: #e3f2fd; color: #0d47a1; } +.stat-badge.none { background: #f5f5f5; color: #777; font-weight: 400; } +.stat-badge.insuf { background: #fce4ec; color: #b71c1c; font-weight: 400; } +.stat-row-table { width: 100%; border-collapse: collapse; font-size: 12px; } +.stat-row-table td { padding: 3px 0; vertical-align: top; } +.stat-row-table .slabel { color: #888; width: 46%; } +.stat-row-table .sval { color: #222; font-variant-numeric: tabular-nums; } +.stat-pval { display: inline-block; font-family: monospace; } +.stat-pval.sig { color: #1a6e2e; font-weight: 600; } +.stat-divider { border: none; border-top: 1px solid #f0f0f0; margin: 8px 0; } +.stat-summary { font-size: 12px; color: #666; margin-bottom: 14px; } +.stat-nodata { color: #999; padding: 40px; text-align: center; font-size: 13px; } .world-popup { font-size: 12px; line-height: 1.35; } .world-popup b { display: block; margin-bottom: 2px; } .world-popup .veg { color: #2e7d32; font-size: 11px; } @@ -173,11 +201,15 @@ body { margin: 0; font: 13px/1.4 system-ui, sans-serif; background: #f5f5f5; col @@ -270,6 +302,9 @@ let fusionMode = "bti"; // bti | itb let miniMapInst = null, miniMarker = null; let worldMapInst = null, worldCluster = null; let worldOverlayOpen = false; +let overlayTab = "map"; +let statsData = null; +let statsYear = null; const maps3 = {}; // { s2, fusion, s3 } Leaflet instances const overlays3 = {}; // current ImageOverlay per map const markers3 = {}; // site dot markers per map @@ -298,6 +333,22 @@ const SERIES_LABELS = { s2: "S2 raw", s3: "S3 raw", }; +const METRIC_META = { + nse: { label: "NSE", full: "Nash–Sutcliffe Efficiency", better: "higher" }, + rmse: { label: "RMSE", full: "Root Mean Square Error", better: "lower" }, + nrmse: { label: "nRMSE", full: "Normalised RMSE", better: "lower" }, + r: { label: "r", full: "Pearson correlation", better: "higher" }, +}; + +const STAT_METRICS = ["nse", "rmse", "nrmse", "r"]; + +const BADGE_CLASS = { + itb: "itb", + bti: "bti", + "no significant difference": "none", + "insufficient data": "insuf", +}; + const INSPECTOR_SERIES = [ { key: "phenocam", label: "PhenoCam", cols: [{ h: "gcc_90", k: "gcc_90" }] }, { key: "bands_s2", label: "S2 reflectance", cols: ["B02","B03","B04"].map(b => ({ h: b, k: b })) }, @@ -503,6 +554,118 @@ function renderSitePanel(meta, cov) { `${rows.join("")}
${species}`; } +function fmtStat(v, decimals = 4) { + return v != null ? v.toFixed(decimals) : "—"; +} + +function fmtPval(p, alpha) { + if (p == null) return "—"; + const cls = p < alpha ? "stat-pval sig" : "stat-pval"; + return `${p.toFixed(4)}`; +} + +function betterOrderLabel(order) { + if (order === "itb") return "ItB better"; + if (order === "bti") return "BtI better"; + if (order === "no significant difference") return "No significant difference"; + if (order === "insufficient data") return "Insufficient data"; + return order; +} + +function updateWorldMeta() { + const sites = manifest?.sites?.[currentYear] || {}; + const n = Object.values(sites).filter(m => m.has_fusion).length; + if (overlayTab === "stats") { + const nPairs = statsData?.metrics?.nse?.n_pairs; + qs("#worldMeta").textContent = nPairs != null + ? `${nPairs} paired site${nPairs === 1 ? "" : "s"} · α=${statsData.alpha ?? 0.05} · ${currentYear}` + : `${n} fusion site${n === 1 ? "" : "s"} · ${currentYear}`; + } else { + qs("#worldMeta").textContent = + `${n} fusion site${n === 1 ? "" : "s"} · ${currentYear}`; + } +} + +function renderStatsPanel(data) { + const panel = qs("#statsPanel"); + const alpha = data.alpha ?? 0.05; + const cards = STAT_METRICS.map(key => { + const meta = METRIC_META[key]; + const m = data.metrics?.[key] || {}; + const badgeCls = BADGE_CLASS[m.better_order] || "none"; + const badge = `${betterOrderLabel(m.better_order)}`; + const row = (label, val) => + `${label}${val}`; + + return `
+

${meta.label} ${meta.full} ${badge}

+
${meta.better} is better
+ + ${row("BtI mean", fmtStat(m.bti_mean))} + ${row("BtI median", fmtStat(m.bti_median))} + ${row("ItB mean", fmtStat(m.itb_mean))} + ${row("ItB median", fmtStat(m.itb_median))} + ${row("Diff (ItB − BtI) mean", fmtStat(m.mean_diff))} + ${row("Diff (ItB − BtI) median", fmtStat(m.median_diff))} +
+
+ + ${row("Wilcoxon W", m.wilcoxon?.statistic ?? "—")} + ${row("Wilcoxon p", fmtPval(m.wilcoxon?.p_value, alpha))} + ${row("Paired t", m.ttest?.statistic ?? "—")} + ${row("Paired t p", fmtPval(m.ttest?.p_value, alpha))} +
+
+ + ${row("Paired sites", m.n_pairs ?? "—")} + ${row("Dropped sites", m.n_dropped ?? "—")} +
+
`; + }).join(""); + + panel.innerHTML = + `
Paired ItB vs BtI test across ${data.n_sites_total ?? "—"} site(s) with Step 5 metrics · significance α=${alpha}
` + + `
${cards}
`; + updateWorldMeta(); +} + +async function loadStatsPanel() { + const panel = qs("#statsPanel"); + panel.innerHTML = '
Loading…
'; + try { + const data = await fetch(`data/statistics_fusion_order/${currentYear}.json`) + .then(r => { if (!r.ok) throw new Error(); return r.json(); }); + statsData = data; + statsYear = currentYear; + renderStatsPanel(data); + } catch { + statsData = null; + statsYear = null; + panel.innerHTML = + '
No statistics file found — run 6-statistics-fusion-order.py first.
'; + updateWorldMeta(); + } +} + +function switchOverlayTab(tab, updateHash = true) { + overlayTab = tab; + document.querySelectorAll(".otab").forEach(btn => + btn.classList.toggle("active", btn.dataset.tab === tab)); + qs("#worldMap").style.display = tab === "map" ? "block" : "none"; + qs("#statsPanel").style.display = tab === "stats" ? "block" : "none"; + + if (tab === "stats") { + if (statsYear !== currentYear || !statsData) loadStatsPanel(); + else updateWorldMeta(); + if (updateHash) setHash("statistics"); + return; + } + + buildWorldMap(); + requestAnimationFrame(() => worldMapInst?.invalidateSize()); + if (updateHash) setHash("worldwide"); +} + // ── init ── async function init() { try { @@ -519,11 +682,18 @@ async function init() { yearSel.addEventListener("change", () => { currentYear = +yearSel.value; + statsData = null; + statsYear = null; buildSiteList(); - if (worldOverlayOpen) buildWorldMap(); + if (worldOverlayOpen) { + if (overlayTab === "stats") loadStatsPanel(); + else buildWorldMap(); + } }); qs("#worldMapBtn").addEventListener("click", () => openWorldOverlay()); + document.querySelectorAll(".otab").forEach(btn => + btn.addEventListener("click", () => switchOverlayTab(btn.dataset.tab))); qs("#worldClose").addEventListener("click", () => closeWorldOverlay()); qs("#worldOverlay").addEventListener("click", e => { if (e.target === qs("#worldOverlay")) closeWorldOverlay(); @@ -551,11 +721,12 @@ async function init() { onHashChange(); } -// ── hash routing (#worldwide, #2025/sitename) ── +// ── hash routing (#worldwide, #statistics, #2025/sitename) ── function parseHash() { const raw = location.hash.replace(/^#/, "").trim(); if (!raw) return { view: null, year: null, site: null }; if (raw === "worldwide") return { view: "worldwide", year: null, site: null }; + if (raw === "statistics") return { view: "statistics", year: null, site: null }; const parts = raw.split("/"); if (parts.length === 2 && /^\d{4}$/.test(parts[0])) return { view: "site", year: +parts[0], site: decodeURIComponent(parts[1]) }; @@ -565,6 +736,7 @@ function parseHash() { function setHash(view, year, site) { let hash = ""; if (view === "worldwide") hash = "worldwide"; + else if (view === "statistics") hash = "statistics"; else if (view === "site" && year && site) hash = `${year}/${encodeURIComponent(site)}`; const next = hash ? `#${hash}` : ""; @@ -574,7 +746,11 @@ function setHash(view, year, site) { function onHashChange() { const { view, year, site } = parseHash(); if (view === "worldwide") { - openWorldOverlay(false); + openWorldOverlay(false, "map"); + return; + } + if (view === "statistics") { + openWorldOverlay(false, "stats"); return; } if (view === "site" && year && site && manifest?.sites?.[year]?.[site]?.has_fusion) { @@ -592,15 +768,13 @@ function onHashChange() { } // ── worldwide map overlay ── -function openWorldOverlay(updateHash = true) { +function openWorldOverlay(updateHash = true, tab = "map") { if (!manifest) return; worldOverlayOpen = true; const overlay = qs("#worldOverlay"); overlay.classList.add("open"); overlay.setAttribute("aria-hidden", "false"); - if (updateHash) setHash("worldwide"); - buildWorldMap(); - requestAnimationFrame(() => worldMapInst?.invalidateSize()); + switchOverlayTab(tab, updateHash); } function closeWorldOverlay(updateHash = true) { @@ -608,7 +782,8 @@ function closeWorldOverlay(updateHash = true) { const overlay = qs("#worldOverlay"); overlay.classList.remove("open"); overlay.setAttribute("aria-hidden", "true"); - if (updateHash && parseHash().view === "worldwide") { + const view = parseHash().view; + if (updateHash && (view === "worldwide" || view === "statistics")) { if (currentSite) setHash("site", currentYear, currentSite); else history.replaceState(null, "", location.pathname + location.search); } @@ -712,7 +887,7 @@ function buildSiteList() { list.appendChild(li); } const h = parseHash(); - if (h.view === "worldwide") return; + if (h.view === "worldwide" || h.view === "statistics") return; if (h.view === "site" && h.year === currentYear && sites[h.site]?.has_fusion) { selectSite(h.site); return;