diff --git a/examples/era_lte_workflow.py b/examples/era_lte_workflow.py new file mode 100644 index 0000000..2cae6c6 --- /dev/null +++ b/examples/era_lte_workflow.py @@ -0,0 +1,264 @@ +"""ERA LTE workflow — drive the Climate Toolkit with ERA's own season windows. + +This is a *workflow built on the toolkit's public API*, not just a run over ERA +coordinates. For each ERA Long-Term Experiment (LTE) site it: + +1. translates the LTE's reported season window(s) into the toolkit's + ``fixed_season`` syntax (``"MM-DD:MM-DD"``, or two comma-separated windows), +2. runs :func:`climate_toolkit.analyze_climate_statistics` with that fixed + season (ERA's season, *not* auto-detection), and +3. optionally compares the toolkit's seasonal rainfall against ERA's reported + rainfall — the validation/benchmarking loop that motivates #110. + +Input CSV contract +------------------ +One row per LTE site-period, with columns (registry aliases in parentheses): + + site_id (Site.ID), lat (Latitude), lon (Longitude), + start_year (Year.start), end_year (Year.end), # required + s1_start, s1_end # optional season-1 window as month-day, "03-01" + s2_start, s2_end # optional second season + reported_rain_mm # optional ERA-reported seasonal rainfall + +Two modes, chosen per row: + +* **Fixed-season** (when ``s1_start``/``s1_end`` are present): the LTE's own + season window is translated to the toolkit's ``fixed_season`` and compared + against ``reported_rain_mm`` if given. Year-crossing seasons (end month before + start month) are supported — the toolkit wraps them to the next year. Season + windows come from ERA's ``Site.Start.S1`` / ``Site.End.S1`` (and S2) fields. +* **Auto** (no season columns): falls back to auto season detection, so ERA's + shipped registry ``data/unique_ltes.csv`` (lat/lon/years only) runs as-is. + +To get the fixed-season + comparison path, export the season windows and +reported rainfall from the ERA ``ERAg`` package tables (they are *not* in the +registry CSV) into the columns above; see https://github.com/ERAgriculture/LTEs. + +Usage +----- + # No credentials needed if --source nasa_power: + python examples/era_lte_workflow.py era_lte.csv --dry-run # show season mapping + python examples/era_lte_workflow.py era_lte.csv --source nasa_power --limit 5 + python examples/era_lte_workflow.py era_lte.csv --out era_lte_compare.csv +""" + +from __future__ import annotations + +import argparse +import re + +import pandas as pd + +from climate_toolkit.climate_statistics import analyze_climate_statistics + +COLUMN_ALIASES = { + "Site.ID": "site_id", + "Latitude": "lat", + "Longitude": "lon", + "Year.start": "start_year", + "Year.end": "end_year", + "Site.Start.S1": "s1_start", + "Site.End.S1": "s1_end", + "Site.Start.S2": "s2_start", + "Site.End.S2": "s2_end", +} + + +def parse_month_day(value) -> str: + """Normalize an ERA season boundary to zero-padded ``"MM-DD"``. + + Accepts ``"MM-DD"``, ``"M-D"``, ``"MM/DD"``, ``"MM.DD"``, 4-digit ``"MMDD"``, + and date-like objects. Raises ``ValueError`` for anything unparseable or out + of range, so bad ERA rows fail loudly rather than silently mis-seasoning. + """ + if value is None or (isinstance(value, float) and pd.isna(value)): + raise ValueError("season boundary is missing") + if hasattr(value, "month") and hasattr(value, "day"): + month, day = int(value.month), int(value.day) + else: + token = str(value).strip() + parts = re.split(r"[-/.]", token) + if len(parts) == 2: + month, day = int(parts[0]), int(parts[1]) + elif re.fullmatch(r"\d{4}", token): + month, day = int(token[:2]), int(token[2:]) + else: + raise ValueError(f"unrecognized month-day value: {value!r}") + if not (1 <= month <= 12 and 1 <= day <= 31): + raise ValueError(f"month-day out of range: {value!r} -> {month:02d}-{day:02d}") + return f"{month:02d}-{day:02d}" + + +def lte_to_fixed_season(row) -> str: + """Build a toolkit ``fixed_season`` string from an LTE row's season windows. + + One season -> ``"MM-DD:MM-DD"``; two seasons -> ``"MM-DD:MM-DD,MM-DD:MM-DD"``. + """ + s1 = f"{parse_month_day(row['s1_start'])}:{parse_month_day(row['s1_end'])}" + s2_start, s2_end = row.get("s2_start"), row.get("s2_end") + has_s2 = not ( + s2_start is None + or s2_end is None + or (isinstance(s2_start, float) and pd.isna(s2_start)) + or str(s2_start).strip() == "" + ) + if has_s2: + s2 = f"{parse_month_day(s2_start)}:{parse_month_day(s2_end)}" + return f"{s1},{s2}" + return s1 + + +def has_season_windows(row) -> bool: + """True if the row carries an S1 season window (enables fixed-season mode). + + ERA's shipped ``unique_ltes.csv`` is registry-only (no season windows), so + the workflow falls back to auto season detection for it; an enriched export + with ``s1_start``/``s1_end`` unlocks the fixed-season + comparison path. + """ + def _present(key): + val = row.get(key) + return not ( + val is None + or (isinstance(val, float) and pd.isna(val)) + or str(val).strip() == "" + ) + + return _present("s1_start") and _present("s1_end") + + +def _year(value): + """Parse an ERA ``Year.start`` / ``Year.end`` cell into an int, or None.""" + try: + return int(float(str(value).strip()[:4])) + except (ValueError, TypeError): + return None + + +def load_sites(csv_path) -> pd.DataFrame: + """Load an ERA LTE export or the shipped registry, applying name aliases. + + Only lat/lon/years are required; season windows are optional (see + :func:`has_season_windows`). ERA's ``unique_ltes.csv`` is cp1252-encoded, so + fall back to that if UTF-8 decoding fails. + """ + try: + df = pd.read_csv(csv_path) + except UnicodeDecodeError: + df = pd.read_csv(csv_path, encoding="cp1252") + df = df.rename(columns=COLUMN_ALIASES) + df["start_year"] = df["start_year"].map(_year) + df["end_year"] = df["end_year"].map(_year) + df = df.dropna(subset=["lat", "lon"]) + df = df[df["start_year"].notna() & df["end_year"].notna()] + return df.reset_index(drop=True) + + +def seasonal_rows(result) -> list[dict]: + """Flatten ``season_statistics`` into per (year, season) toolkit metrics.""" + rows = [] + for s in result.get("season_statistics", []): + precip = s.get("precipitation") or {} + wb = s.get("water_balance") or {} + length = s.get("length_days") + rainy = precip.get("rainy_days") + rows.append( + { + "year": s.get("year"), + "season_number": s.get("season_number"), + "tk_rain_total_mm": precip.get("total_mm"), + "tk_rainy_days": rainy, + "tk_dry_days": (length - rainy) + if (length is not None and rainy is not None) + else None, + "tk_NDWS": wb.get("NDWS"), + "tk_WRSI": wb.get("WRSI"), + } + ) + return rows + + +def run_site(row, source: str) -> list[dict]: + """Run the toolkit for one LTE site, return per-season metrics. + + Uses ERA's own season window when present (``fixed_season``), otherwise + falls back to auto season detection so the shipped registry CSV still runs. + """ + kwargs = dict( + location_coord=(float(row["lat"]), float(row["lon"])), + start_year=int(row["start_year"]), + end_year=int(row["end_year"]), + source=source, + verbose=False, + ) + if has_season_windows(row): + fixed = lte_to_fixed_season(row) + kwargs["fixed_season"] = fixed + else: + fixed = None + result = analyze_climate_statistics(**kwargs) + out = [] + for tk in seasonal_rows(result): + record = { + "site_id": row.get("site_id"), + "mode": "fixed" if fixed else "auto", + "fixed_season": fixed, + **tk, + } + if "reported_rain_mm" in row and not pd.isna(row.get("reported_rain_mm")): + reported = float(row["reported_rain_mm"]) + record["reported_rain_mm"] = reported + if tk["tk_rain_total_mm"] is not None: + record["rain_delta_mm"] = tk["tk_rain_total_mm"] - reported + out.append(record) + return out + + +def run(csv_path, source="nasa_power", out_csv="era_lte_compare.csv", limit=None, dry_run=False): + sites = load_sites(csv_path) + if limit: + sites = sites.head(limit) + print(f"{len(sites)} ERA LTE site-periods") + + if dry_run: + for _, row in sites.iterrows(): + if has_season_windows(row): + try: + print(f" {row.get('site_id')}: fixed_season={lte_to_fixed_season(row)!r}") + except ValueError as exc: + print(f" {row.get('site_id')}: [bad season window] {exc}") + else: + print(f" {row.get('site_id')}: auto season detection (no season window)") + return sites + + out_rows = [] + for _, row in sites.iterrows(): + try: + out_rows.extend(run_site(row, source)) + except Exception as exc: # one bad site shouldn't stop the workflow + print(f"[skip] {row.get('site_id')}: {exc}") + + frame = pd.DataFrame(out_rows) + frame.to_csv(out_csv, index=False) + print(f"wrote {out_csv} ({len(frame)} rows)") + if "rain_delta_mm" in frame and frame["rain_delta_mm"].notna().any(): + deltas = frame["rain_delta_mm"].dropna() + print( + "toolkit vs ERA-reported rainfall — " + f"bias={deltas.mean():.1f} mm | MAE={deltas.abs().mean():.1f} mm | n={len(deltas)}" + ) + return frame + + +def main(): + parser = argparse.ArgumentParser(description="ERA LTE workflow on the Climate Toolkit.") + parser.add_argument("csv", help="ERA LTE export CSV (see module docstring for columns)") + parser.add_argument("--source", default="nasa_power", help="Climate source (default: nasa_power)") + parser.add_argument("--out", default="era_lte_compare.csv", help="Output CSV path") + parser.add_argument("--limit", type=int, default=None, help="Only process the first N sites") + parser.add_argument("--dry-run", action="store_true", help="Show season mapping without fetching") + args = parser.parse_args() + run(args.csv, args.source, args.out, args.limit, args.dry_run) + + +if __name__ == "__main__": + main() diff --git a/tests/test_era_lte_workflow.py b/tests/test_era_lte_workflow.py new file mode 100644 index 0000000..321abd0 --- /dev/null +++ b/tests/test_era_lte_workflow.py @@ -0,0 +1,79 @@ +"""Unit tests for the ERA LTE workflow's season-window mapper. + +Covers the credential-free core: translating ERA season boundaries into the +toolkit's ``fixed_season`` syntax. The actual ``analyze_climate_statistics`` +calls are exercised elsewhere / by hand — these tests pin the pure translation. +""" + +import os +import sys +import unittest +from datetime import date + +REPO_ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +if REPO_ROOT not in sys.path: + sys.path.insert(0, REPO_ROOT) + +from examples.era_lte_workflow import ( + has_season_windows, + lte_to_fixed_season, + parse_month_day, +) + + +class ParseMonthDayTests(unittest.TestCase): + def test_accepted_formats_normalize_to_mm_dd(self): + cases = { + "03-01": "03-01", + "3-1": "03-01", + "03/01": "03-01", + "03.01": "03-01", + "0301": "03-01", + "12-31": "12-31", + } + for raw, expected in cases.items(): + with self.subTest(raw=raw): + self.assertEqual(parse_month_day(raw), expected) + + def test_date_like_object(self): + self.assertEqual(parse_month_day(date(2020, 3, 1)), "03-01") + + def test_invalid_values_raise(self): + for bad in ("13-01", "03-40", "abc", "", None): + with self.subTest(bad=bad): + with self.assertRaises(ValueError): + parse_month_day(bad) + + +class LteToFixedSeasonTests(unittest.TestCase): + def test_single_season(self): + row = {"s1_start": "03-01", "s1_end": "05-31"} + self.assertEqual(lte_to_fixed_season(row), "03-01:05-31") + + def test_two_seasons(self): + row = {"s1_start": "03-01", "s1_end": "05-31", "s2_start": "10-01", "s2_end": "12-15"} + self.assertEqual(lte_to_fixed_season(row), "03-01:05-31,10-01:12-15") + + def test_year_crossing_passes_through(self): + # cessation month before onset month; the toolkit wraps to next year. + row = {"s1_start": "10-01", "s1_end": "03-31"} + self.assertEqual(lte_to_fixed_season(row), "10-01:03-31") + + def test_empty_s2_ignored(self): + row = {"s1_start": "03-01", "s1_end": "05-31", "s2_start": "", "s2_end": ""} + self.assertEqual(lte_to_fixed_season(row), "03-01:05-31") + + +class ModeSelectionTests(unittest.TestCase): + def test_season_windows_present_enables_fixed_mode(self): + self.assertTrue(has_season_windows({"s1_start": "03-01", "s1_end": "05-31"})) + + def test_registry_row_falls_back_to_auto(self): + # ERA's shipped unique_ltes.csv has no season columns. + for row in ({}, {"s1_start": "", "s1_end": ""}, {"s1_start": None, "s1_end": None}): + with self.subTest(row=row): + self.assertFalse(has_season_windows(row)) + + +if __name__ == "__main__": + unittest.main()