🌡️ Downloading ECMWF S2S Temperature Forecasts¶

#!/usr/bin/env python
"""
Download ECMWF S2S realtime **daily-averaged 2m temperature (T2M)**.

Example:
python download_ecmwf_s2s_t2m_daily.py \
  --outdir data/s2s_ecmwf \
  --outfile s2s_ecmwf_daily_t2m_2025-11-01_ea.nc \
  --date 2025-11-01 \
  --lead-days 5 \
  --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 ECMWF S2S *daily* step string with intervals:
    "0-24/24-48/48-72/..."

    According to the S2S docs, this is the way to request
    daily-mean/daily-accumulated products.
    """
    periods = []
    for i in range(lead_days):
        start = i * 24
        end = (i + 1) * 24
        periods.append(f"{start}-{end}")
    return "/".join(periods)


def retrieve_s2s_t2m_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-averaged
    2m temperature (T2M, param=2t), control forecast.
    """
    server = ECMWFDataServer()

    # S2S ECMWF daily-averaged typically 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.")

    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": "2t",      # ONLY 2m temperature
        "date": date_str,   # YYYY-MM-DD
        "time": "00:00:00",
        "step": steps,      # "0-24/24-48/..."
        "area": area,
        "format": fmt,
        "target": out_path,
        # Allow partial retrieval instead of failing with "Expected N, got M"
        "expect": "any",
    }

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

    print("[info] Submitting ECMWF S2S T2M daily request…")
    print("[info] Request:", request)
    server.retrieve(request)
    print(f"[info] Download finished → {out_path}")


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

def parse_args() -> argparse.Namespace:
    p = argparse.ArgumentParser(
        description="Download ECMWF S2S realtime daily-averaged 2m temperature (T2M)."
    )
    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 averages 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 dir 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_t2m_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 get_recent_s2s_dates(n=4):
    """Get the most recent n valid S2S dates (Mondays and Thursdays)."""
    today = datetime.now()
    dates = []

    # Go back up to 30 days to find valid dates
    for i in range(30):
        check_date = today - timedelta(days=i)
        if check_date.weekday() in [0, 3]:  # Monday=0, Thursday=3
            dates.append(check_date.strftime("%Y-%m-%d"))
            if len(dates) >= n:
                break

    return dates

print("Recent S2S dates:", get_recent_s2s_dates())

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_t2m_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 temperature variable (may be 't2m' or '2t')
temp_var = 't2m' if 't2m' in ds.data_vars else '2t' if '2t' in ds.data_vars else list(ds.data_vars)[0]
temp = ds[temp_var]

# Convert from Kelvin to Celsius if needed
if temp.max() > 200:  # Likely in Kelvin
    temp_c = temp - 273.15
    units = '°C'
else:
    temp_c = temp
    units = temp.attrs.get('units', '°C')

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

# Plot Week 1 mean temperature
fig, ax = plt.subplots(figsize=(10, 8))
week1_mean = temp_c.isel(time=slice(0, 7)).mean(dim='time')
week1_mean.plot(ax=ax, cmap='RdYlBu_r', vmin=15, vmax=35)
ax.set_title('S2S Week 1 Mean Daily Temperature')
plt.savefig('s2s_week1_temp.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(temp_c.time):
        weekly_mean = temp_c.isel(time=slice(start, end)).mean(dim='time')
        weekly_mean.plot(ax=axes[i], cmap='RdYlBu_r', vmin=15, vmax=35, 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_temp_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 = temp_c.sel(latitude=lat_point, longitude=lon_point, method='nearest')
point_data.plot(marker='o', figsize=(12, 4), color='orangered')
plt.title(f'S2S 30-Day Temperature Forecast 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_temp_timeseries.png', dpi=150, bbox_inches='tight')
plt.show()

# Calculate temperature anomaly (if climatology available)
# climatology = ...  # Load your climatology
# anomaly = temp_c - climatology

# Dataset attributes (example)
{
    'Conventions': 'CF-1.6',
    'history': 'Retrieved from ECMWF S2S',
    'institution': 'ECMWF'
}

# Variable attributes
{
    'units': 'K',
    'long_name': '2 metre temperature',
    'standard_name': 'air_temperature'
}

import xarray as xr

# Load both datasets
ds_temp = xr.open_dataset('data/s2s_ecmwf/s2s_ecmwf_t2m_ethiopia_30day.nc')
ds_precip = xr.open_dataset('data/s2s_ecmwf/s2s_ecmwf_tp_ethiopia_30day.nc')

# Rename variables for consistency
ds_temp = ds_temp.rename({'2t': 't2m'} if '2t' in ds_temp else {})
ds_precip = ds_precip.rename({'tp': 'precip'} if 'tp' in ds_precip else {})

# Convert temperature to Celsius
if ds_temp.t2m.max() > 200:
    ds_temp['t2m'] = ds_temp['t2m'] - 273.15
    ds_temp['t2m'].attrs['units'] = 'degC'

# Convert precipitation to mm/day if in meters
if ds_precip.precip.max() < 1:
    ds_precip['precip'] = ds_precip['precip'] * 1000
    ds_precip['precip'].attrs['units'] = 'mm/day'

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

# Verify
print(ds_combined)

# Save combined file
ds_combined.to_netcdf('data/s2s_ecmwf/s2s_combined_ethiopia_30day.nc')
print("Saved combined temperature and precipitation dataset")