scripts/data_acquisition/fetch_ecmwf_forecast.py

"""
Fetches real, physics-based short-to-medium-range forecast data from
Open-Meteo's dedicated ECMWF endpoint -- built specifically for genuine
flood early-warning lead time, which (per the recent design discussion)
structurally requires a real precipitation forecast, not just discharge
history extrapolated forward.

CONFIRMED real API details (not guessed): base endpoint
https://api.open-meteo.com/v1/ecmwf, HTTP GET, no authentication, CC BY
4.0 licence, up to 10 real forecast days. Also confirmed to exist (not
yet used here): an archived-individual-model-run-by-init-time feature
explicitly described by Open-Meteo as built for bias-correction
pipelines without look-ahead leakage -- exactly what a later step would
need to calibrate this forecast against our own real climatology (see
climate_climatology.py). The exact endpoint for that archive wasn't
directly confirmed in this pass -- flagged for a follow-up check, not
guessed at here.

Usage:
    python -m scripts.fetch_ecmwf_forecast --check
    python -m scripts.fetch_ecmwf_forecast --stations-path datasets/station_elevations.csv
"""
import argparse
from pathlib import Path

import pandas as pd
import requests

BASE_URL = "https://api.open-meteo.com/v1/ecmwf"

# Daily aggregates, matching this project's existing daily-resolution
# pipeline convention (ERA5 climatology etc.) rather than pulling hourly
# data that would need its own separate aggregation step.
DAILY_VARIABLES = [
    "precipitation_sum", "temperature_2m_max", "temperature_2m_min",
    "wind_speed_10m_max", "shortwave_radiation_sum", "et0_fao_evapotranspiration",
]


def check_access() -> bool:
    """Single real request against a real coordinate in the study area
    (matches H602021010's real location) to confirm the endpoint and
    variable names actually work before looping over every station."""
    params = {
        "latitude": 48.748705218, "longitude": 0.58156195,
        "daily": ",".join(DAILY_VARIABLES), "forecast_days": 5,
        "timezone": "Europe/Paris",
    }
    print(f"Testing {BASE_URL} with a real study-area coordinate...")
    try:
        resp = requests.get(BASE_URL, params=params, timeout=30)
        resp.raise_for_status()
        data = resp.json()
    except Exception as e:
        print(f"FAILED: {e}")
        return False

    if "daily" not in data:
        print(f"FAILED: response has no 'daily' key. Full response: {data}")
        return False

    print("OK -- real forecast returned:")
    for var in DAILY_VARIABLES:
        if var in data["daily"]:
            print(f"  {var}: {data['daily'][var][:3]}... ({len(data['daily'][var])} real day(s))")
        else:
            print(f"  WARNING: '{var}' not found in response -- variable name may be wrong. "
                  f"Real keys returned: {list(data['daily'].keys())}")
    return True


def fetch_forecast(lat: float, lon: float, forecast_days: int = 10, max_retries: int = 3) -> pd.DataFrame:
    """
    Retries on read timeout -- confirmed as a real, not hypothetical,
    intermittent failure mode against the live API: 2 of 27 real
    stations hit a 30s read timeout on an otherwise-clean run, while
    the other 25 succeeded immediately. A short exponential backoff
    handles this automatically rather than requiring a manual re-run
    for whichever stations happen to be unlucky that run.
    """
    import time

    params = {
        "latitude": lat, "longitude": lon,
        "daily": ",".join(DAILY_VARIABLES), "forecast_days": forecast_days,
        "timezone": "Europe/Paris",
    }
    last_error = None
    for attempt in range(max_retries):
        try:
            resp = requests.get(BASE_URL, params=params, timeout=30)
            resp.raise_for_status()
            data = resp.json()
            break
        except requests.exceptions.RequestException as e:
            last_error = e
            if attempt < max_retries - 1:
                wait = 2 ** attempt  # 1s, 2s, 4s
                print(f"    attempt {attempt + 1}/{max_retries} failed ({e}), retrying in {wait}s...")
                time.sleep(wait)
    else:
        raise last_error

    if "daily" not in data:
        raise ValueError(f"No 'daily' data in response for ({lat}, {lon}): {data}")

    df = pd.DataFrame(data["daily"])
    df["latitude"] = lat
    df["longitude"] = lon
    return df


def main() -> None:
    parser = argparse.ArgumentParser(description="Fetch real ECMWF forecast data from Open-Meteo")
    parser.add_argument("--check", action="store_true")
    parser.add_argument("--stations-path", type=Path, default=Path("datasets/station_elevations.csv"))
    parser.add_argument("--forecast-days", type=int, default=10)
    parser.add_argument("--output", type=Path, default=Path("datasets/ecmwf_forecast.csv"))
    args = parser.parse_args()

    if args.check:
        check_access()
        return

    stations = pd.read_csv(args.stations_path)
    lat_col = "latitude" if "latitude" in stations.columns else "lat"
    lon_col = "longitude" if "longitude" in stations.columns else "lon"

    all_dfs = []
    for _, row in stations.iterrows():
        code = row.get("station_code", f"{row[lat_col]:.3f}_{row[lon_col]:.3f}")
        print(f"Fetching real forecast for {code}...")
        try:
            df = fetch_forecast(row[lat_col], row[lon_col], args.forecast_days)
            df["station_code"] = code
            all_dfs.append(df)
        except Exception as e:
            print(f"  FAILED for {code}: {e}")

    if not all_dfs:
        print("No real forecast data fetched at all -- check connectivity, or run --check first.")
        return

    combined = pd.concat(all_dfs, ignore_index=True)
    args.output.parent.mkdir(parents=True, exist_ok=True)
    combined.to_csv(args.output, index=False)
    print(f"\nSaved real forecast data for {combined['station_code'].nunique()} station(s), "
          f"{args.forecast_days} real day(s) each, to {args.output}")


if __name__ == "__main__":
    main()