added gap validation.

gap_validation/phenology_offsets.py

@@ -0,0 +1,146 @@
+"""TIMESAT transition dates on gap-degraded fusion series vs PhenoCam reference."""
+
+from __future__ import annotations
+
+import argparse
+import json
+from datetime import datetime
+from pathlib import Path
+
+from phenology_timesat import (
+    build_yraw_three_years,
+    phenocam_phenology_path,
+    run_timesat_phenology_from_yraw,
+)
+
+from gap_validation.batch_spatial import (
+    PRIMARY_SEASON,
+    _best_bti_from_metrics,
+    _parse_scenario,
+    _site_positions,
+)
+from gap_validation.calendar import load_manifest, validation_dir
+from gap_validation.temporal_pc import _fusion_gcc_timeseries
+
+
+def _day_offset(iso_a: str | None, iso_b: str | None) -> int | None:
+    if not iso_a or not iso_b:
+        return None
+    try:
+        a = datetime.strptime(iso_a[:10], "%Y-%m-%d").date()
+        b = datetime.strptime(iso_b[:10], "%Y-%m-%d").date()
+        return abs((a - b).days)
+    except ValueError:
+        return None
+
+
+def _timesat_transitions(by_date: dict[str, float], season: int) -> dict[str, str | None]:
+    y1, y2, y3 = season - 1, season, season + 1
+    yraw, _mode = build_yraw_three_years(by_date, y1, y2, y3)
+    out = run_timesat_phenology_from_yraw(yraw, (y1, y2, y3))
+    return {
+        "green_up": out.get("green_up_50pct_date"),
+        "green_down": out.get("green_down_50pct_date"),
+    }
+
+
+def _temporal_fusion_dir(
+    site: str, season: int, gap_days: int, transition: str, scenario_key: str
+) -> Path:
+    strategy, sigma, mode = _parse_scenario(scenario_key)
+    sig = 30 if sigma == 30 else 20
+    return (
+        validation_dir(site, season)
+        / "temporal"
+        / f"gap_{gap_days}_{transition}"
+        / f"{strategy}_sigma{sig}_{mode}"
+    )
+
+
+def compute_offsets_for_site(
+    site: str,
+    season: int,
+    site_position: tuple[float, float],
+    *,
+    gap_days_list: tuple[int, ...] = (15, 30),
+) -> list[dict]:
+    base = Path(f"data/{site}/{season}")
+    metrics_path = base / "metrics.json"
+    scenario_key = _best_bti_from_metrics(metrics_path)
+    if not scenario_key:
+        return []
+    ref_path = phenocam_phenology_path(site, season)
+    reference = (
+        json.loads(ref_path.read_text(encoding="utf-8")) if ref_path.is_file() else {}
+    )
+    manifest = load_manifest(site, season)
+    rows: list[dict] = []
+    for entry in manifest["entries"]:
+        gd = entry.get("gap_days")
+        tr = entry.get("transition")
+        if gd not in gap_days_list or tr not in ("green_up", "green_down"):
+            continue
+        fusion_dir = _temporal_fusion_dir(site, season, gd, tr, scenario_key)
+        if not fusion_dir.is_dir():
+            continue
+        _, _, mode = _parse_scenario(scenario_key)
+        ts = _fusion_gcc_timeseries(fusion_dir, site_position, mode)
+        if len(ts) < 10:
+            continue
+        fused = _timesat_transitions(ts, season)
+        ref_key = (
+            "green_up_50pct_date"
+            if tr == "green_up"
+            else "green_down_50pct_date"
+        )
+        ref_date = reference.get(ref_key)
+        fused_date = fused.get("green_up" if tr == "green_up" else "green_down")
+        rows.append(
+            {
+                "site_name": site,
+                "season": season,
+                "transition": tr,
+                "gap_days": gd,
+                "scenario": scenario_key,
+                "reference_date": ref_date,
+                "fused_date": fused_date,
+                "abs_day_offset": _day_offset(fused_date, ref_date),
+                "window_start": entry.get("window_start"),
+                "window_end": entry.get("window_end"),
+            }
+        )
+    return rows
+
+
+def write_phenology_offsets(
+    site: str,
+    season: int,
+    site_position: tuple[float, float],
+    *,
+    gap_days_list: tuple[int, ...] = (15, 30),
+) -> Path:
+    rows = compute_offsets_for_site(
+        site, season, site_position, gap_days_list=gap_days_list
+    )
+    out = validation_dir(site, season) / "gap_phenology_offsets.json"
+    payload = {"site_name": site, "season": season, "records": rows}
+    out.write_text(json.dumps(payload, indent=2) + "\n", encoding="utf-8")
+    return out
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description="Gap fusion TIMESAT offsets vs PhenoCam.")
+    ap.add_argument("--data-dir", type=Path, default=Path("data"))
+    ap.add_argument("--sites-geojson", type=Path, default=Path("data/sites.geojson"))
+    args = ap.parse_args()
+    positions = _site_positions(args.sites_geojson)
+    for site, season in sorted(PRIMARY_SEASON.items()):
+        pos = positions.get(site)
+        if not pos:
+            continue
+        p = write_phenology_offsets(site, season, pos)
+        print(p)
+
+
+if __name__ == "__main__":
+    main()