🌧️ Downloading ECMWF S2S Precipitation Forecasts¶

#!/usr/bin/env python
"""
Download ECMWF S2S realtime **daily total precipitation** (24-hour accumulations).

Example:
python download_ecmwf_s2s_tp_daily.py \
  --outdir data/s2s_ecmwf \
  --outfile s2s_ecmwf_daily_tp_2025-11-01_ea.nc \
  --date 2025-11-01 \
  --lead-days 30 \
  --lat-min 3 --lat-max 15 --lon-min 33 --lon-max 48
"""

import argparse
import os
from datetime import datetime
from ecmwfapi import ECMWFDataServer


# --------------------------------------------------------------------------- #
# Helpers
# --------------------------------------------------------------------------- #

def build_daily_step_string(lead_days: int) -> str:
    """
    Build a MARS step string for 24-hour lead times.

    For lead_days = 5 → "24/48/72/96/120"
    (i.e. end of each 24-hour accumulation period).
    """
    if lead_days < 1:
        raise ValueError("lead_days must be >= 1")

    steps = [str(24 * (i + 1)) for i in range(lead_days)]
    return "/".join(steps)


def retrieve_s2s_tp_daily(
    date_str: str,
    lead_days: int,
    lat_min: float,
    lat_max: float,
    lon_min: float,
    lon_max: float,
    out_path: str,
    grid: str | None = None,
    fmt: str = "netcdf",
) -> None:
    """
    Submit an ECMWF S2S request for **daily total precipitation** (control forecast).

    The request uses:
      class  = s2
      stream = enfo
      type   = cf  (control forecast)
      origin = ecmf
      param  = tp  (internally mapped to paramId 228228 = 24h total precip)
    """
    # Safety: ECMWF S2S daily product goes up to 46 days
    if lead_days < 1 or lead_days > 46:
        raise ValueError("lead_days must be between 1 and 46 for ECMWF S2S.")

    server = ECMWFDataServer()

    steps = build_daily_step_string(lead_days)

    # ECMWF area string is N/W/S/E
    area = f"{lat_max}/{lon_min}/{lat_min}/{lon_max}"

    request = {
        "class": "s2",
        "dataset": "s2s",
        "expver": "prod",
        "origin": "ecmf",
        "model": "glob",
        "levtype": "sfc",
        "stream": "enfo",
        "type": "cf",          # control forecast
        "number": "0",         # control member
        "param": "tp",         # daily total precip (internally paramId 228228)
        "date": date_str,      # YYYY-MM-DD, forecast start date
        "time": "00:00:00",
        "step": steps,         # e.g. "24/48/72/..."
        "area": area,
        "format": fmt,
        "target": out_path,
        # Avoid failing with "Expected N, got M" if some steps are missing
        "expect": "any",
    }

    # Optional horizontal interpolation grid, e.g. "1.5/1.5" or "0.5/0.5"
    if grid is not None:
        request["grid"] = grid

    print("[info] Submitting S2S daily TP request to ECMWF...")
    print("[info] Request:", request)
    server.retrieve(request)
    print(f"[info] Download finished: {out_path}")
    print("[info] NetCDF will contain a variable like 'tp' with dims (time, lat, lon).")


# --------------------------------------------------------------------------- #
# CLI
# --------------------------------------------------------------------------- #

def parse_args() -> argparse.Namespace:
    p = argparse.ArgumentParser(
        description="Download ECMWF S2S realtime daily total precipitation (24h accum)."
    )
    p.add_argument("--outdir", required=True, help="Output directory")
    p.add_argument("--outfile", required=True, help="Output filename (NetCDF)")
    p.add_argument(
        "--date",
        required=True,
        help="Forecast initial date (YYYY-MM-DD) for S2S start",
    )
    p.add_argument(
        "--lead-days",
        type=int,
        required=True,
        help="Number of lead days (1–46) of daily totals to retrieve",
    )
    p.add_argument("--lat-min", type=float, required=True)
    p.add_argument("--lat-max", type=float, required=True)
    p.add_argument("--lon-min", type=float, required=True)
    p.add_argument("--lon-max", type=float, required=True)
    p.add_argument(
        "--grid",
        default=None,
        help="Optional output grid resolution 'lat/lon', e.g. '1.5/1.5' or '0.5/0.5'",
    )
    p.add_argument(
        "--fmt",
        default="netcdf",
        choices=["netcdf", "grib"],
        help="Output format (default: netcdf)",
    )
    return p.parse_args()


def main() -> None:
    args = parse_args()

    # Ensure output directory exists
    os.makedirs(args.outdir, exist_ok=True)
    out_path = os.path.join(args.outdir, args.outfile)

    # Basic date sanity check
    try:
        datetime.strptime(args.date, "%Y-%m-%d")
    except ValueError as exc:
        raise SystemExit(f"Invalid --date '{args.date}', expected YYYY-MM-DD") from exc

    retrieve_s2s_tp_daily(
        date_str=args.date,
        lead_days=args.lead_days,
        lat_min=args.lat_min,
        lat_max=args.lat_max,
        lon_min=args.lon_min,
        lon_max=args.lon_max,
        out_path=out_path,
        grid=args.grid,
        fmt=args.fmt,
    )


if __name__ == "__main__":
    main()

from datetime import datetime, timedelta

def next_s2s_date(from_date=None):
    """Find the next valid S2S forecast date (Monday or Thursday)."""
    if from_date is None:
        from_date = datetime.now()

    weekday = from_date.weekday()  # 0=Monday, 3=Thursday

    # Days until next Monday or Thursday
    days_to_monday = (7 - weekday) % 7 if weekday != 0 else 0
    days_to_thursday = (3 - weekday) % 7 if weekday != 3 else 0

    if days_to_monday == 0:
        return from_date.strftime("%Y-%m-%d")
    elif days_to_thursday == 0:
        return from_date.strftime("%Y-%m-%d")
    elif days_to_monday < days_to_thursday:
        return (from_date + timedelta(days=days_to_monday)).strftime("%Y-%m-%d")
    else:
        return (from_date + timedelta(days=days_to_thursday)).strftime("%Y-%m-%d")

print(f"Next S2S date: {next_s2s_date()}")

import xarray as xr
import matplotlib.pyplot as plt
import numpy as np

# Open the forecast file
ds = xr.open_dataset('data/s2s_ecmwf/s2s_ecmwf_tp_ethiopia_30day.nc')

# Display dataset information
print(ds)

# Check dimensions
print(f"Lead times: {len(ds.time) if 'time' in ds.dims else len(ds.step)}")
print(f"Latitude range: {float(ds.latitude.min()):.2f} to {float(ds.latitude.max()):.2f}")
print(f"Longitude range: {float(ds.longitude.min()):.2f} to {float(ds.longitude.max()):.2f}")

# Get precipitation variable (may be 'tp' or 'tprate')
precip_var = 'tp' if 'tp' in ds.data_vars else list(ds.data_vars)[0]
precip = ds[precip_var]

# Check units - S2S tp is typically in meters, convert to mm
if precip.max() < 1:  # Likely in meters
    precip_mm = precip * 1000
    units = 'mm/day'
else:
    precip_mm = precip
    units = precip.attrs.get('units', 'unknown')

print(f"Precipitation range: {float(precip_mm.min()):.1f} to {float(precip_mm.max()):.1f} {units}")

# Plot Week 1 mean precipitation
fig, ax = plt.subplots(figsize=(10, 8))
week1_mean = precip_mm.isel(time=slice(0, 7)).mean(dim='time')
week1_mean.plot(ax=ax, cmap='Blues', vmin=0, vmax=20)
ax.set_title('S2S Week 1 Mean Daily Precipitation')
plt.savefig('s2s_week1_precip.png', dpi=150, bbox_inches='tight')
plt.show()

# Plot weekly evolution
weeks = [0, 7, 14, 21, 28]
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
for i, (start, end) in enumerate(zip(weeks[:-1], weeks[1:])):
    if end <= len(precip_mm.time):
        weekly_mean = precip_mm.isel(time=slice(start, end)).mean(dim='time')
        weekly_mean.plot(ax=axes[i], cmap='Blues', vmin=0, vmax=20, add_colorbar=False)
        axes[i].set_title(f'Week {i+1}')
        axes[i].set_xlabel('')
        axes[i].set_ylabel('')
plt.tight_layout()
plt.savefig('s2s_weekly_evolution.png', dpi=150, bbox_inches='tight')
plt.show()

# Time series for a point
lat_point, lon_point = 9.0, 38.7  # Addis Ababa
point_data = precip_mm.sel(latitude=lat_point, longitude=lon_point, method='nearest')
point_data.plot(marker='o', figsize=(12, 4))
plt.title(f'S2S 30-Day Precipitation Forecast for Addis Ababa')
plt.ylabel('Precipitation (mm/day)')
plt.xlabel('Lead Time')
plt.grid(True, alpha=0.3)
plt.savefig('s2s_timeseries.png', dpi=150, bbox_inches='tight')
plt.show()

# To request ensemble members, modify the request:
request = {
    ...
    "type": "pf",          # perturbed forecast (ensemble)
    "number": "1/2/3/.../50",  # all 50 members
    ...
}

request = {
    "class": "s2",           # S2S project
    "dataset": "s2s",        # S2S database
    "origin": "ecmf",        # ECMWF model
    "stream": "enfo",        # Ensemble forecast
    "type": "cf",            # Control forecast
    "param": "tp",           # Total precipitation
    "step": "24/48/72/...",  # Daily accumulation steps
    "area": "N/W/S/E",       # Bounding box
    "format": "netcdf",      # Output format
}