🌡️ Downloading ECMWF S2S Ensemble Temperature Forecasts¶

#!/usr/bin/env python
"""
Download ECMWF S2S realtime daily-averaged 2m temperature (T2M)
for multiple ensemble members and compute the ensemble mean.

Example:
python download_ecmwf_s2s_t2m_ensemble_dailymean.py \
  --outdir data/s2s_ecmwf \
  --outfile s2s_ecmwf_t2m_ensmean_2025-12-01_ea.nc \
  --date 2025-12-01 \
  --lead-days 30 \
  --members 10 \
  --lat-min 3 --lat-max 15 --lon-min 33 --lon-max 48
"""

import argparse
import os
from datetime import datetime
from ecmwfapi import ECMWFDataServer
import xarray as xr
import numpy as np


def build_daily_step_string(lead_days: int) -> str:
    """Build ECMWF S2S daily step string: '0-24/24-48/48-72/...'"""
    periods = [f"{i*24}-{(i+1)*24}" for i in range(lead_days)]
    return "/".join(periods)


def retrieve_member_t2m_daily(
    date_str: str,
    lead_days: int,
    lat_min: float,
    lat_max: float,
    lon_min: float,
    lon_max: float,
    member: int,
    out_path: str,
    grid: str | None = None,
    fmt: str = "netcdf",
) -> bool:
    """
    Retrieve ECMWF S2S daily-averaged 2m temperature for one ensemble member.

    Returns True if successful, False otherwise.
    """
    server = ECMWFDataServer()
    if lead_days < 1 or lead_days > 46:
        raise ValueError("lead_days must be between 1 and 46 for ECMWF S2S.")

    steps = build_daily_step_string(lead_days)
    area = f"{lat_max}/{lon_min}/{lat_min}/{lon_max}"

    req = {
        "class": "s2",
        "dataset": "s2s",
        "expver": "prod",
        "origin": "ecmf",
        "model": "glob",
        "levtype": "sfc",
        "stream": "enfo",
        "type": "pf",           # perturbed forecast
        "number": str(member),  # ensemble member id (1..N)
        "param": "2t",
        "date": date_str,
        "time": "00:00:00",
        "step": steps,
        "area": area,
        "format": fmt,
        "target": out_path,
        "expect": "any",
    }
    if grid is not None:
        req["grid"] = grid

    print(f"[info] Requesting member={member}...")

    try:
        server.retrieve(req)
    except Exception as exc:
        print(f"[warn] ECMWF request failed for member={member}: {exc}")
        return False

    # Check file validity
    if (not os.path.exists(out_path)) or (os.path.getsize(out_path) < 500):
        print(f"[warn] Output for member={member} looks empty or too small")
        return False

    print(f"[done] Member {member} → {out_path}")
    return True


def compute_ensemble_mean(member_files, out_path, to_celsius=True):
    """
    Compute ensemble mean across member files (preserving lead days).

    Parameters
    ----------
    member_files : list
        List of paths to member NetCDF files
    out_path : str
        Output path for ensemble mean
    to_celsius : bool
        Convert from Kelvin to Celsius (default True)
    """
    print(f"[info] Merging {len(member_files)} members...")

    datasets = []
    valid_members = []

    for i, f in enumerate(member_files, 1):
        if not os.path.exists(f):
            print(f"[warn] Member file missing: {f}")
            continue
        try:
            ds = xr.open_dataset(f)
            # Handle different variable names
            if 't2m' in ds:
                da = ds['t2m']
            elif '2t' in ds:
                da = ds['2t']
            else:
                print(f"[warn] No temperature variable in {f}")
                continue
            datasets.append(da)
            valid_members.append(i)
        except Exception as exc:
            print(f"[warn] Failed to open {f}: {exc}")
            continue

    if not datasets:
        raise RuntimeError("No valid member files found")

    # Stack along member dimension
    da_stack = xr.concat(
        datasets,
        dim=xr.DataArray(valid_members, dims="member"),
    )

    # Compute ensemble mean
    ens_mean = da_stack.mean("member")

    # Convert to Celsius if needed
    if to_celsius and float(ens_mean.max()) > 200:  # Likely Kelvin
        ens_mean = ens_mean - 273.15
        ens_mean.attrs["units"] = "degC"
        print("[info] Converted temperature from Kelvin to Celsius")

    ens_mean.attrs["long_name"] = "Ensemble mean daily 2m temperature"
    ens_mean.name = "t2m"

    # Save to NetCDF
    ds_out = ens_mean.to_dataset(name="t2m")
    ds_out.attrs["title"] = "ECMWF S2S Ensemble Mean 2m Temperature"
    ds_out.attrs["source"] = "ECMWF S2S (param=2t, perturbed forecasts)"
    ds_out.attrs["n_members"] = len(valid_members)

    ds_out.to_netcdf(out_path)
    print(f"[done] Ensemble mean saved → {out_path}")


def compute_ensemble_statistics(member_files, out_path, to_celsius=True):
    """
    Compute full ensemble statistics (mean, std, percentiles).

    Use this for probabilistic products.
    """
    print(f"[info] Computing ensemble statistics from {len(member_files)} members...")

    datasets = []
    valid_members = []

    for i, f in enumerate(member_files, 1):
        if not os.path.exists(f):
            continue
        try:
            ds = xr.open_dataset(f)
            da = ds['t2m'] if 't2m' in ds else ds['2t']
            datasets.append(da)
            valid_members.append(i)
        except:
            continue

    if not datasets:
        raise RuntimeError("No valid member files found")

    # Stack along member dimension
    da_stack = xr.concat(datasets, dim=xr.DataArray(valid_members, dims="member"))

    # Convert to Celsius if needed
    if to_celsius and float(da_stack.max()) > 200:
        da_stack = da_stack - 273.15

    # Compute statistics
    ens_mean = da_stack.mean("member")
    ens_std = da_stack.std("member")
    ens_min = da_stack.min("member")
    ens_max = da_stack.max("member")
    ens_p10 = da_stack.quantile(0.1, dim="member")
    ens_p25 = da_stack.quantile(0.25, dim="member")
    ens_median = da_stack.quantile(0.5, dim="member")
    ens_p75 = da_stack.quantile(0.75, dim="member")
    ens_p90 = da_stack.quantile(0.9, dim="member")

    # Create output dataset
    ds_out = xr.Dataset({
        't2m_mean': ens_mean,
        't2m_std': ens_std,
        't2m_min': ens_min,
        't2m_max': ens_max,
        't2m_p10': ens_p10,
        't2m_p25': ens_p25,
        't2m_median': ens_median,
        't2m_p75': ens_p75,
        't2m_p90': ens_p90,
    })

    # Add attributes
    units = "degC" if to_celsius else "K"
    for var in ds_out.data_vars:
        ds_out[var].attrs["units"] = units

    ds_out.attrs["title"] = "ECMWF S2S Ensemble Temperature Statistics"
    ds_out.attrs["n_members"] = len(valid_members)

    ds_out.to_netcdf(out_path)
    print(f"[done] Ensemble statistics saved → {out_path}")


def main():
    p = argparse.ArgumentParser(
        description="Download ECMWF S2S ensemble daily-mean 2m temperature (T2M)."
    )
    p.add_argument("--outdir", required=True, help="Output directory")
    p.add_argument("--outfile", required=True, help="Final ensemble-mean filename")
    p.add_argument("--date", required=True, help="Forecast date (YYYY-MM-DD)")
    p.add_argument("--lead-days", type=int, required=True, help="Lead days (1-46)")
    p.add_argument("--members", type=int, default=10, help="Number of members (1-50)")
    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="Output grid (e.g., '0.5/0.5')")
    p.add_argument("--fmt", default="netcdf", choices=["netcdf", "grib"])
    p.add_argument("--keep-kelvin", action="store_true", 
                   help="Keep temperature in Kelvin (default: convert to Celsius)")
    p.add_argument("--full-stats", action="store_true",
                   help="Compute full statistics (mean, std, percentiles)")
    args = p.parse_args()

    # Validate date
    try:
        datetime.strptime(args.date, "%Y-%m-%d")
    except ValueError:
        raise SystemExit(f"Invalid date format: {args.date}")

    os.makedirs(args.outdir, exist_ok=True)
    member_files = []

    # Download each member
    for m in range(1, args.members + 1):
        fpath = os.path.join(args.outdir, f"t2m_member{m:02d}.nc")
        ok = retrieve_member_t2m_daily(
            args.date,
            args.lead_days,
            args.lat_min, args.lat_max,
            args.lon_min, args.lon_max,
            member=m,
            out_path=fpath,
            grid=args.grid,
            fmt=args.fmt,
        )
        if ok:
            member_files.append(fpath)
        else:
            print(f"[warn] Skipping member={m}")

    if not member_files:
        raise SystemExit("No member files downloaded successfully")

    # Compute ensemble mean or full statistics
    out_nc = os.path.join(args.outdir, args.outfile)

    if args.full_stats:
        compute_ensemble_statistics(
            member_files, out_nc, 
            to_celsius=not args.keep_kelvin
        )
    else:
        compute_ensemble_mean(
            member_files, out_nc,
            to_celsius=not args.keep_kelvin
        )


if __name__ == "__main__":
    main()

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

# Open the ensemble mean file
ds = xr.open_dataset('data/s2s_ensemble/s2s_ensmean_t2m_ethiopia.nc')

# Display dataset information
print(ds)

# Check dimensions
print(f"Lead times: {len(ds.time) if 'time' in ds.dims else 'N/A'}")
print(f"Temperature units: {ds.t2m.attrs.get('units', 'unknown')}")

# Check temperature range
print(f"Temperature range: {float(ds.t2m.min()):.1f} to {float(ds.t2m.max()):.1f}")

# Plot Week 1 ensemble mean temperature
fig, ax = plt.subplots(figsize=(10, 8))
week1_mean = ds.t2m.isel(time=slice(0, 7)).mean(dim='time')
week1_mean.plot(ax=ax, cmap='RdYlBu_r', vmin=15, vmax=35)
ax.set_title('S2S Ensemble Mean: Week 1 Daily Temperature')
plt.savefig('s2s_t2m_ensemble_week1.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 = ds.t2m.sel(latitude=lat_point, longitude=lon_point, method='nearest')
point_data.plot(marker='o', figsize=(12, 4), color='orangered')
plt.title(f'S2S Ensemble Mean Temperature for Addis Ababa')
plt.ylabel('Temperature (°C)')
plt.xlabel('Lead Time')
plt.grid(True, alpha=0.3)
plt.axhline(y=point_data.mean(), color='gray', linestyle='--', label='Mean')
plt.legend()
plt.savefig('s2s_t2m_ensemble_timeseries.png', dpi=150, bbox_inches='tight')
plt.show()

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

# Load all member files
member_files = sorted(glob.glob('data/s2s_ensemble/t2m_member*.nc'))
print(f"Found {len(member_files)} member files")

# Stack all members
datasets = []
for i, f in enumerate(member_files, 1):
    ds = xr.open_dataset(f)
    da = ds['t2m'] if 't2m' in ds else ds['2t']
    # Convert to Celsius if needed
    if da.max() > 200:
        da = da - 273.15
    ds_exp = da.expand_dims({'member': [i]})
    datasets.append(ds_exp)

# Combine all members
t2m_all = xr.concat(datasets, dim='member')
print(t2m_all)

# Compute ensemble statistics
ens_mean = t2m_all.mean(dim='member')
ens_std = t2m_all.std(dim='member')
ens_p10 = t2m_all.quantile(0.1, dim='member')
ens_p90 = t2m_all.quantile(0.9, dim='member')

# Compute probability of warm anomaly (e.g., > 25°C)
threshold = 25  # °C
prob_warm = (t2m_all > threshold).mean(dim='member') * 100

# Create output dataset
ds_out = xr.Dataset({
    't2m_mean': ens_mean,
    't2m_std': ens_std,
    't2m_p10': ens_p10,
    't2m_p90': ens_p90,
    'prob_above_25C': prob_warm,
})

# Save
ds_out.to_netcdf('data/s2s_ensemble/s2s_t2m_probabilistic.nc')
print("Saved probabilistic temperature products")

# Visualize ensemble spread
fig, axes = plt.subplots(1, 4, figsize=(16, 4))

for i, week in enumerate([0, 7, 14, 21]):
    if week + 7 <= len(ens_std.time):
        weekly_std = ens_std.isel(time=slice(week, week+7)).mean(dim='time')
        im = weekly_std.plot(ax=axes[i], cmap='YlOrRd', vmin=0, vmax=3, add_colorbar=False)
        axes[i].set_title(f'Week {i+1}')
        axes[i].set_xlabel('')
        axes[i].set_ylabel('')

plt.suptitle('Ensemble Spread (Standard Deviation) in °C')
plt.tight_layout()
plt.savefig('s2s_t2m_spread.png', dpi=150, bbox_inches='tight')
plt.show()

# Spaghetti plot for a single point
lat_point, lon_point = 9.0, 38.7
point_all = t2m_all.sel(latitude=lat_point, longitude=lon_point, method='nearest')

plt.figure(figsize=(12, 5))
for m in range(len(member_files)):
    plt.plot(point_all.isel(member=m), color='gray', alpha=0.3, linewidth=0.5)
plt.plot(ens_mean.sel(latitude=lat_point, longitude=lon_point, method='nearest'), 
         color='red', linewidth=2, label='Ensemble Mean')
plt.fill_between(
    range(len(ens_mean.time)),
    ens_p10.sel(latitude=lat_point, longitude=lon_point, method='nearest'),
    ens_p90.sel(latitude=lat_point, longitude=lon_point, method='nearest'),
    color='red', alpha=0.2, label='10th-90th percentile'
)
plt.xlabel('Lead Day')
plt.ylabel('Temperature (°C)')
plt.title(f'S2S Temperature Ensemble for Addis Ababa')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('s2s_t2m_spaghetti.png', dpi=150, bbox_inches='tight')
plt.show()

import xarray as xr

# Load ensemble means
ds_temp = xr.open_dataset('data/s2s_ensemble/s2s_ensmean_t2m_ethiopia.nc')
ds_precip = xr.open_dataset('data/s2s_ensemble/s2s_ensmean_tp_ethiopia.nc')

# Merge datasets
ds_combined = xr.merge([ds_temp, ds_precip])

# Verify
print(ds_combined)
# Variables: t2m, tp

# Save combined file
ds_combined.to_netcdf('data/s2s_ensemble/s2s_combined_ethiopia.nc')
print("Saved combined ensemble mean dataset")