#!/usr/bin/env python3 """ Download ERA5 hourly temperature variables from CDS and compute daily summaries. Supported Variables ------------------- - 2m_temperature (t2m) -> Daily Mean - minimum_2m_temperature_since_previous_post_processing (mn2t) -> Daily Min - maximum_2m_temperature_since_previous_post_processing (mx2t) -> Daily Max Key features ------------ - Monthly hourly download to avoid CDS "cost limits exceeded" - Computes daily aggregations (mean, min, or max) per month - Region bounding box support (N/W/S/E) - Optional unit conversion to Celsius - Optional cleanup of hourly files - Optional merge of all daily files into one NetCDF """ import argparse import os import glob from datetime import datetime import numpy as np import xarray as xr import cdsapi # Mapping from friendly names to CDS variable names and NetCDF short names VAR_MAP = { "t2m": { "cds_name": "2m_temperature", "short_name": "t2m", "long_name": "2m temperature", "agg": "mean" }, "tmin": { "cds_name": "minimum_2m_temperature_since_previous_post_processing", "short_name": "mn2t", "long_name": "minimum 2m temperature", "agg": "min" }, "tmax": { "cds_name": "maximum_2m_temperature_since_previous_post_processing", "short_name": "mx2t", "long_name": "maximum 2m temperature", "agg": "max" } } # --------------------------------------------------------------------------- # # Helpers: time handling # --------------------------------------------------------------------------- # def find_time_dim(ds: xr.Dataset) -> str: for cand in ("time", "valid_time"): if cand in ds.coords or cand in ds.dims: return cand for name, coord in ds.coords.items(): if "time" in name.lower(): return name try: if np.issubdtype(coord.dtype, np.datetime64): return name except Exception: pass raise KeyError("Could not find a time dimension/coordinate.") def standardise_time_for_resample(ds: xr.Dataset) -> xr.Dataset: time_dim = find_time_dim(ds) if time_dim in ds.dims and time_dim not in ds.coords: if time_dim in ds.variables: ds = ds.assign_coords({time_dim: ds[time_dim]}) if time_dim != "time": ds = ds.rename({time_dim: "time"}) if "time" in ds.coords: if not np.issubdtype(ds["time"].dtype, np.datetime64): ds = xr.decode_cf(ds) return ds # --------------------------------------------------------------------------- # # ERA5 retrieval and processing # --------------------------------------------------------------------------- # def build_monthly_request(year: int, month: int, area: list[float], variables: list[str]) -> dict: year_str = f"{year:04d}" month_str = f"{month:02d}" days = [f"{d:02d}" for d in range(1, 32)] times = [f"{h:02d}:00" for h in range(0, 24)] cds_vars = [VAR_MAP[v]["cds_name"] for v in variables] return { "product_type": "reanalysis", "variable": cds_vars, "year": year_str, "month": month_str, "day": days, "time": times, "area": area, # [N, W, S, E] "format": "netcdf", } def retrieve_hourly_month(year: int, month: int, out_path: str, area: list[float], variables: list[str]) -> None: client = cdsapi.Client() request = build_monthly_request(year, month, area, variables) print(f"[info] Requesting ERA5 hourly for {year:04d}-{month:02d} (Vars: {variables})...") result = client.retrieve("reanalysis-era5-single-levels", request) result.download(out_path) def compute_daily_summaries(hourly_path: str, daily_path: str, variables: list[str], to_celsius: bool) -> None: ds = xr.open_dataset(hourly_path) ds = standardise_time_for_resample(ds) daily_results = [] for v in variables: short_name = VAR_MAP[v]["short_name"] long_name = VAR_MAP[v]["long_name"] agg_type = VAR_MAP[v]["agg"] if short_name not in ds: # Try to find by long name if short name mismatch found = False for dv in ds.data_vars: if ds[dv].attrs.get("long_name", "").lower() == long_name.lower(): ds = ds.rename({dv: short_name}) found = True break if not found: print(f"[warning] Variable {short_name} not found in {hourly_path}. Skipping.") continue # Aggregate if agg_type == "mean": da_daily = ds[short_name].resample(time="1D").mean() elif agg_type == "min": da_daily = ds[short_name].resample(time="1D").min() elif agg_type == "max": da_daily = ds[short_name].resample(time="1D").max() # Unit conversion if to_celsius: da_daily = da_daily - 273.15 da_daily.attrs["units"] = "degC" else: da_daily.attrs["units"] = "K" da_daily.attrs["long_name"] = f"Daily {agg_type} {long_name}" daily_results.append(da_daily.to_dataset(name=short_name)) if not daily_results: ds.close() return ds_daily = xr.merge(daily_results) ds_daily.attrs["source"] = "ERA5 reanalysis (CDS)" ds_daily.attrs["processing"] = "Hourly download; daily aggregation computed with xarray" encoding = {v: {"zlib": True, "complevel": 4} for v in ds_daily.data_vars} ds_daily.to_netcdf(daily_path, encoding=encoding) ds.close() ds_daily.close() def merge_daily_files(daily_files: list[str], out_path: str) -> None: if not daily_files: return ds = xr.open_mfdataset(daily_files, combine="by_coords") if "time" in ds.coords: ds = ds.sortby("time") encoding = {v: {"zlib": True, "complevel": 4} for v in ds.data_vars} ds.to_netcdf(out_path, encoding=encoding) ds.close() # --------------------------------------------------------------------------- # # CLI # --------------------------------------------------------------------------- # def main() -> None: p = argparse.ArgumentParser(description="Download ERA5 hourly temperature variables and compute daily summaries.") p.add_argument("--start-year", type=int, required=True) p.add_argument("--end-year", type=int, required=True) 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("--outdir", required=True) p.add_argument("--variables", nargs="+", default=["t2m", "tmin", "tmax"], help="Vars to download: t2m, tmin, tmax") p.add_argument("--merge-outfile", default=None) p.add_argument("--to-celsius", action="store_true") p.add_argument("--keep-hourly", action="store_true") p.add_argument("--delete-hourly", action="store_true") p.add_argument("--keep-monthly-daily", action="store_true") args = p.parse_args() os.makedirs(args.outdir, exist_ok=True) area = [args.lat_max, args.lon_min, args.lat_min, args.lon_max] hourly_dir = os.path.join(args.outdir, "hourly_monthly") daily_dir = os.path.join(args.outdir, "daily_monthly") os.makedirs(hourly_dir, exist_ok=True) os.makedirs(daily_dir, exist_ok=True) daily_files = [] for year in range(args.start_year, args.end_year + 1): for month in range(1, 13): hourly_path = os.path.join(hourly_dir, f"era5_temp_hourly_{year:04d}_{month:02d}.nc") daily_path = os.path.join(daily_dir, f"era5_temp_daily_{year:04d}_{month:02d}.nc") if not os.path.exists(hourly_path): retrieve_hourly_month(year, month, hourly_path, area, args.variables) if not os.path.exists(daily_path): print(f"[info] Computing daily summaries for {year:04d}-{month:02d}...") compute_daily_summaries(hourly_path, daily_path, args.variables, args.to_celsius) daily_files.append(daily_path) if args.delete_hourly or (not args.keep_hourly and not args.delete_hourly): try: os.remove(hourly_path) except OSError: pass if args.merge_outfile: merged_path = os.path.join(args.outdir, args.merge_outfile) merge_daily_files(daily_files, merged_path) if not args.keep_monthly_daily: for f in daily_files: try: os.remove(f) except OSError: pass print("[info] Done.") if __name__ == "__main__": main()