From 424f105ab8deb48e71cbbf158559cd2027c79c33 Mon Sep 17 00:00:00 2001 From: jome1 Date: Tue, 13 Oct 2026 11:31:52 +0200 Subject: [PATCH 1/4] add access to bias-corrected ERA5 data from Copernicus --- RELEASE_NOTES.rst | 4 + atlite/data.py | 31 ++- atlite/datasets/__init__.py | 3 +- atlite/datasets/era5_corrected.py | 322 ++++++++++++++++++++++++++++++ atlite/wind.py | 4 +- 5 files changed, 356 insertions(+), 8 deletions(-) create mode 100644 atlite/datasets/era5_corrected.py diff --git a/RELEASE_NOTES.rst b/RELEASE_NOTES.rst index 2f48f25e..dd9c94d9 100644 --- a/RELEASE_NOTES.rst +++ b/RELEASE_NOTES.rst @@ -20,6 +20,10 @@ Upcoming Release **Features** * Add support for Python 3.14. +* Add the ``era5-corrected`` module providing the bias-corrected C3S Energy ERA5 + 100 m wind speed ([``sis-energy-global-reanalysis``](https://cds.climate.copernicus.eu/datasets/sis-energy-global-reanalysis?tab=overview)) as ``wnd100m_corrected`` + and ``wnd100m_corrected_extrapolated``. Combine it with ``era5`` to get the shear exponent in order to extrapolate the corrected wind speed. Assign + the preferred variable to ``cutout.data["wnd100m"]`` to use it in ``cutout.wind``. **Bug fixes** diff --git a/atlite/data.py b/atlite/data.py index 2064a8dd..c7541495 100644 --- a/atlite/data.py +++ b/atlite/data.py @@ -5,8 +5,10 @@ from __future__ import annotations +import gc import logging import os +import time from functools import wraps from pathlib import Path from shutil import rmtree @@ -32,6 +34,19 @@ logger = logging.getLogger(__name__) +def _rename_with_retry(source: Path, target: Path, attempts: int = 20) -> None: + """Rename a file, retrying while Windows still holds handles on it.""" + for attempt in range(attempts): + try: + source.rename(target) + return + except PermissionError: + if attempt == attempts - 1: + raise + gc.collect() + time.sleep(0.25) + + def get_features( cutout: Cutout, module: str, @@ -233,12 +248,12 @@ def cutout_prepare( dask_kwargs : dict, default {} Dict with keyword arguments passed to `dask.compute`. monthly_requests : bool, optional - If True, the data is requested on a monthly basis in ERA5. This is useful for - large cutouts, where the data is requested in smaller chunks. The - default is False + If True, ERA5 and ERA5-corrected data are requested on a monthly basis. + This is useful for large cutouts, where the data is requested in smaller + chunks. The default is False. concurrent_requests : bool, optional - If True, the monthly data requests are posted concurrently. - Only has an effect if `monthly_requests` is True. The default is False. + If True, the time-chunked data requests are posted concurrently. The + default is False. Returns ------- @@ -319,10 +334,14 @@ def cutout_prepare( write_job.compute(**dask_kwargs) else: write_job.compute(**dask_kwargs) + del write_job if cutout.path.exists(): cutout.data.close() + # release lazily opened handles of the old file (Windows cannot + # replace files that are still open) + gc.collect() cutout.path.unlink() - Path(tmp).rename(cutout.path) + _rename_with_retry(Path(tmp), cutout.path) cutout.data = xr.open_dataset(cutout.path, chunks=cutout.chunks) diff --git a/atlite/datasets/__init__.py b/atlite/datasets/__init__.py index f678e3d2..322f795e 100644 --- a/atlite/datasets/__init__.py +++ b/atlite/datasets/__init__.py @@ -4,10 +4,11 @@ """atlite datasets.""" -from atlite.datasets import era5, era5_edh, gebco, glofas, sarah +from atlite.datasets import era5, era5_corrected, era5_edh, gebco, glofas, sarah modules = { "era5": era5, + "era5-corrected": era5_corrected, "era5-edh": era5_edh, "gebco": gebco, "glofas": glofas, diff --git a/atlite/datasets/era5_corrected.py b/atlite/datasets/era5_corrected.py new file mode 100644 index 00000000..0680848d --- /dev/null +++ b/atlite/datasets/era5_corrected.py @@ -0,0 +1,322 @@ +# SPDX-FileCopyrightText: Contributors to atlite +# +# SPDX-License-Identifier: MIT +""" +Module for downloading data from the C3S Energy bias-corrected ERA5 reanalysis. + +The data is retrieved from the Climate Data Store dataset +``sis-energy-global-reanalysis``, see +https://cds.climate.copernicus.eu/datasets/sis-energy-global-reanalysis + +Only the 100 m wind speed is available, in two variants which are both stored +in the cutout (feature ``wind_corrected``): + +* ``wnd100m_corrected``: bias-corrected 100 m wind speed + (CDS variable ``100m_wind_speed``) +* ``wnd100m_corrected_extrapolated``: wind speed extrapolated to 100 m with a + power law from the bias-corrected lower level wind speed + (CDS variable ``100m_wind_speed_extrapolated``) + +Combine the module with ``era5`` to get the roughness needed for hub heights +other than 100 m. Afterwards choose which 100 m wind speed is used by +:py:meth:`atlite.Cutout.wind` by assigning it to ``wnd100m``. + +Example +------- +>>> cutout = atlite.Cutout( +... "corrected.nc", +... module=["era5", "era5-corrected"], +... x=slice(-12, 19), +... y=slice(36, 52), +... time="2010-01-01", +... ) +>>> cutout.prepare() +>>> cutout.data["wnd100m"] = cutout.data["wnd100m_corrected_extrapolated"] +""" + +from __future__ import annotations + +import logging +import os +import tempfile +import zipfile +from contextlib import nullcontext +from typing import TYPE_CHECKING, Any + +import numpy as np +import pandas as pd +import xarray as xr +from dask import compute, delayed + +from atlite.datasets.cds_helper import _area, sanitize_chunks +from atlite.gis import maybe_swap_spatial_dims + +if TYPE_CHECKING: + from dask.utils import SerializableLock + + from atlite._types import PathLike + +logger = logging.getLogger(__name__) + +# Model and CRS Settings +crs = 4326 + +features = {"wind_corrected": ["wnd100m_corrected", "wnd100m_corrected_extrapolated"]} + +DATASET = "sis-energy-global-reanalysis" +# CDS variable -> (tag in the delivered file name, name in the cutout) +WIND_VARIABLES = { + "100m_wind_speed": ("_mbc_", "wnd100m_corrected"), + "100m_wind_speed_extrapolated": ("_pwl_", "wnd100m_corrected_extrapolated"), +} + + +def _rename_and_clean_coords(ds: xr.Dataset, add_lon_lat: bool = True) -> xr.Dataset: + """ + Standardize coordinates to atlite conventions. + + The dataset uses longitudes in [0, 360), which are shifted to [-180, 180) + and sorted, and are renamed to x/y. + + Parameters + ---------- + ds : xr.Dataset + Raw dataset with 'longitude' and 'latitude' coordinates. + add_lon_lat : bool, optional + Whether to add 'lon'/'lat' as coordinate aliases. Default True. + + Returns + ------- + xr.Dataset + Dataset with standardized coordinates. + """ + ds = ds.rename({"longitude": "x", "latitude": "y"}) + x = ((ds.x.astype(float) + 180) % 360) - 180 + ds = ds.assign_coords(x=np.round(x, 5), y=np.round(ds.y.astype(float), 5)) + ds = ds.sortby("x") + ds = maybe_swap_spatial_dims(ds) + if add_lon_lat: + ds = ds.assign_coords(lon=ds.coords["x"], lat=ds.coords["y"]) + return ds + + +def retrieve_data( + chunks: dict[str, int] | None = None, + tmpdir: PathLike | None = None, + lock: SerializableLock | None = None, + **updates: Any, +) -> xr.Dataset: + """ + Download data from the CDS API and return it as an xarray Dataset. + + If you want to track the state of your request go to + https://cds.climate.copernicus.eu/requests?tab=all + + Note that the CDS delivers complete months, irrespective of the requested + days. + + Parameters + ---------- + chunks : dict[str, int] or None, optional + Dask chunk specification of the returned dataset. + tmpdir : PathLike or None, optional + Directory for temporary download files. + lock : SerializableLock or None, optional + Lock for thread-safe file creation. + **updates : Any + Additional CDS API request parameters. Must include at least + 'variable', 'year', 'month' and 'day'. + + Returns + ------- + xr.Dataset + Downloaded data. + """ + import cdsapi + + request: dict[str, Any] = { + "spatial_coverage": "global", + "spatial_resolution": ["0_25_degree"], + "temporal_resolution": ["1_hour"], + "version": "1_00", + "file_format": "netcdf_4", + } + request.update(updates) + + assert {"year", "month", "day", "variable"}.issubset(request), ( + "Need to specify at least 'variable', 'year', 'month' and 'day'" + ) + + logger.debug("Requesting %s with API request: %s", DATASET, request) + + client = cdsapi.Client( + info_callback=logger.debug, debug=logging.root.level <= logging.DEBUG + ) + result = client.retrieve(DATASET, request) + + cm: Any = nullcontext() if lock is None else lock + with cm: + fd, archive = tempfile.mkstemp(suffix=".zip", dir=tmpdir) + os.close(fd) + + year_str = ( + ",".join(request["year"]) + if isinstance(request["year"], list) + else str(request["year"]) + ) + month_str = ( + ",".join(request["month"]) + if isinstance(request["month"], list) + else str(request["month"]) + ) + timestr = f"{year_str}-{month_str}" + varstr = "\n\t".join(f"{v} ({timestr})" for v in request["variable"]) + logger.info("CDS: Downloading variables\n\t%s\n", varstr) + result.download(archive) + + # The files are read into memory so that they can be removed right away, + # which is required on Windows where open files cannot be deleted. + # The zip can hold several files; the file name tags bias-corrected (mbc) + # and power-law extrapolated (pwl) data. When multiple months are requested, + # the zip contains one NetCDF file per month. + try: + datasets = [] + with zipfile.ZipFile(archive) as zf: + names = list(dict.fromkeys(n for n in zf.namelist() if n.endswith(".nc"))) + for variable in request["variable"]: + tag, new_name = WIND_VARIABLES[variable] + members = sorted(n for n in names if tag in n) + assert members, ( + f"Expected at least one NetCDF file for {variable}, got none" + ) + var_datasets = [] + for member in members: + fd, target = tempfile.mkstemp(suffix=".nc", dir=tmpdir) + try: + with os.fdopen(fd, "wb") as out: + out.write(zf.read(member)) + with xr.open_dataset(target) as raw: + var_datasets.append(raw.load().rename({"ws100": new_name})) + finally: + os.remove(target) + if len(var_datasets) == 1: + datasets.append(var_datasets[0]) + else: + combined = xr.concat(var_datasets, dim="time") + combined = combined.sortby("time") + datasets.append(combined) + finally: + os.remove(archive) + + ds = xr.merge(datasets) + ds = _rename_and_clean_coords(ds) + chunks = sanitize_chunks(chunks, time="time", x="x", y="y") + return ds.chunk(chunks) if isinstance(chunks, dict) else ds + + +def get_data( + cutout: Any, + feature: str, + tmpdir: PathLike, + lock: SerializableLock | None = None, + monthly_requests: bool = False, + concurrent_requests: bool = False, + **creation_parameters: Any, +) -> xr.Dataset: + """ + Download bias-corrected ERA5 data for a given feature. + + Parameters + ---------- + cutout : Cutout + Cutout object defining the spatial and temporal extent. + feature : str + Feature to retrieve. Only 'wind_corrected' is supported. + tmpdir : PathLike + Directory for temporary download files. + lock : SerializableLock or None, optional + Lock for thread-safe file creation. + monthly_requests : bool, optional + If True, split API requests by month. Otherwise, request one year at a + time. Default False. + concurrent_requests : bool, optional + If True, use dask.delayed for parallel downloads. Default False. + **creation_parameters : Any + Additional cutout parameters, ignored. + + Notes + ----- + Timeframes longer than 12 months are rejected because the CDS API request + would exceed its size limit. Each wind variable is downloaded separately. + + Returns + ------- + xr.Dataset + Data for the requested feature, aligned to cutout coordinates. + """ + if feature not in features: + raise NotImplementedError(f"Feature {feature!r} not supported by era5-corrected") + if (cutout.dx, cutout.dy) != (0.25, 0.25): + raise NotImplementedError( + "era5-corrected only supports the native 0.25°×0.25° grid." + ) + + coords = cutout.coords + retrieval_params: dict[str, Any] = { + "area": _area(coords), + "chunks": cutout.chunks, + "tmpdir": tmpdir, + "lock": lock, + } + + logger.info("Requesting data for feature %s...", feature) + + # Keep requests to at most 12 months to stay within the CDS request limit. + time = coords["time"].to_index() + if len(time) and time[-1] > time[0] + pd.DateOffset(months=12): + logger.warning( + "The era5-corrected API supports cutouts of at most 12 months; " + "the requested timeframe is longer." + ) + raise ValueError("era5-corrected timeframe cannot exceed 12 months.") + + # The CDS delivers whole months regardless of 'day', so request a single + # day per month. Group those months by year unless monthly requests were + # explicitly requested. + if monthly_requests: + time_chunks = [ + { + "year": [str(month.year)], + "month": [month.strftime("%m")], + "day": ["01"], + } + for month in time.to_period("M").unique() + ] + else: + time_chunks = [ + { + "year": [str(year)], + "month": [ + month.strftime("%m") + for month in time[time.year == year].to_period("M").unique() + ], + "day": ["01"], + } + for year in time.year.unique() + ] + + def retrieve_once(t: dict[str, Any]) -> xr.Dataset: + datasets = [ + retrieve_data(variable=[variable], **retrieval_params, **t) + for variable in WIND_VARIABLES + ] + return xr.merge(datasets) + + if concurrent_requests: + datasets = compute(*[delayed(retrieve_once)(c) for c in time_chunks]) + else: + datasets = map(retrieve_once, time_chunks) + + result = xr.concat(datasets, dim="time").sel(time=coords["time"]) + assert isinstance(result, xr.Dataset) + return result diff --git a/atlite/wind.py b/atlite/wind.py index a7b28645..f0b4bef5 100644 --- a/atlite/wind.py +++ b/atlite/wind.py @@ -80,7 +80,9 @@ def extrapolate_wind_speed( if from_height is None: heights: NDArray = np.asarray([ - int(str(s)[3:-1]) for s in ds if re.match(r"wnd\d+m", str(s)) + int(m.group(1)) + for s in ds + if (m := re.fullmatch(r"wnd(\d+)m", str(s))) ]) if len(heights) == 0: From f6a4455b03ce1133bee637ae85e4d2e08aa40f65 Mon Sep 17 00:00:00 2001 From: jome1 Date: Tue, 13 Oct 2026 13:53:01 +0200 Subject: [PATCH 2/4] add download for wind power law coefficients (shear exponent) --- RELEASE_NOTES.rst | 10 +++-- atlite/datasets/era5_corrected.py | 72 ++++++++++++++++++++++++++++--- 2 files changed, 71 insertions(+), 11 deletions(-) diff --git a/RELEASE_NOTES.rst b/RELEASE_NOTES.rst index dd9c94d9..7c8f1b1d 100644 --- a/RELEASE_NOTES.rst +++ b/RELEASE_NOTES.rst @@ -20,10 +20,12 @@ Upcoming Release **Features** * Add support for Python 3.14. -* Add the ``era5-corrected`` module providing the bias-corrected C3S Energy ERA5 - 100 m wind speed ([``sis-energy-global-reanalysis``](https://cds.climate.copernicus.eu/datasets/sis-energy-global-reanalysis?tab=overview)) as ``wnd100m_corrected`` - and ``wnd100m_corrected_extrapolated``. Combine it with ``era5`` to get the shear exponent in order to extrapolate the corrected wind speed. Assign - the preferred variable to ``cutout.data["wnd100m"]`` to use it in ``cutout.wind``. +* Add the ``era5-corrected`` module providing the bias-corrected ([``sis-energy-global-reanalysis``](https://cds.climate.copernicus.eu/datasets/sis-energy-global-reanalysis?tab=overview)) + ERA5 100 m wind speed as ``wnd100m_corrected`` and ``wnd100m_corrected_extrapolated``, plus the power law coefficients as + ``wnd_shear_exp_corrected`` (a month-by-hour climatology, fetched in a separate request). + Assign the preferred variables to ``cutout.data["wnd100m"]`` and + ``cutout.data["wnd_shear_exp"]`` to use them in ``cutout.wind``. Combine with ``era5`` for ``roughness`` + (logarithmic extrapolation). **Bug fixes** diff --git a/atlite/datasets/era5_corrected.py b/atlite/datasets/era5_corrected.py index 0680848d..eab45f34 100644 --- a/atlite/datasets/era5_corrected.py +++ b/atlite/datasets/era5_corrected.py @@ -9,17 +9,23 @@ https://cds.climate.copernicus.eu/datasets/sis-energy-global-reanalysis Only the 100 m wind speed is available, in two variants which are both stored -in the cutout (feature ``wind_corrected``): +in the cutout together with the power law coefficients (feature +``wind_corrected``): * ``wnd100m_corrected``: bias-corrected 100 m wind speed (CDS variable ``100m_wind_speed``) * ``wnd100m_corrected_extrapolated``: wind speed extrapolated to 100 m with a power law from the bias-corrected lower level wind speed (CDS variable ``100m_wind_speed_extrapolated``) +* ``wnd_shear_exp_corrected``: power law coefficient (alpha) matching the + corrected wind speed (CDS variable ``wind_power_law_coefficients``). It is a + month-by-hour climatology, repeated for every year of the cutout, and is + requested separately from the wind speeds. -Combine the module with ``era5`` to get the roughness needed for hub heights -other than 100 m. Afterwards choose which 100 m wind speed is used by -:py:meth:`atlite.Cutout.wind` by assigning it to ``wnd100m``. +Hub heights other than 100 m need either the power law coefficient (use +``method="power"`` and assign ``wnd_shear_exp_corrected`` to ``wnd_shear_exp``) +or the roughness from ``era5`` (logarithmic method). Choose which 100 m wind +speed is used by :py:meth:`atlite.Cutout.wind` by assigning it to ``wnd100m``. Example ------- @@ -61,7 +67,13 @@ # Model and CRS Settings crs = 4326 -features = {"wind_corrected": ["wnd100m_corrected", "wnd100m_corrected_extrapolated"]} +features = { + "wind_corrected": [ + "wnd100m_corrected", + "wnd100m_corrected_extrapolated", + "wnd_shear_exp_corrected", + ] +} DATASET = "sis-energy-global-reanalysis" # CDS variable -> (tag in the delivered file name, name in the cutout) @@ -69,6 +81,15 @@ "100m_wind_speed": ("_mbc_", "wnd100m_corrected"), "100m_wind_speed_extrapolated": ("_pwl_", "wnd100m_corrected_extrapolated"), } +# Power law coefficients (alpha) are a month-by-hour climatology, not a time +# series, and are requested separately to keep the request size small. +POWER_LAW_VARIABLE = "wind_power_law_coefficients" +POWER_LAW_NAME = "wnd_shear_exp_corrected" +# CDS variable -> (file name tag, name in the delivered file, name in the cutout) +_RAW_VARIABLES = { + **{v: (tag, "ws100", name) for v, (tag, name) in WIND_VARIABLES.items()}, + POWER_LAW_VARIABLE: ("ALPH-coeff", "alpha", POWER_LAW_NAME), +} def _rename_and_clean_coords(ds: xr.Dataset, add_lon_lat: bool = True) -> xr.Dataset: @@ -184,7 +205,7 @@ def retrieve_data( with zipfile.ZipFile(archive) as zf: names = list(dict.fromkeys(n for n in zf.namelist() if n.endswith(".nc"))) for variable in request["variable"]: - tag, new_name = WIND_VARIABLES[variable] + tag, raw_name, new_name = _RAW_VARIABLES[variable] members = sorted(n for n in names if tag in n) assert members, ( f"Expected at least one NetCDF file for {variable}, got none" @@ -196,7 +217,7 @@ def retrieve_data( with os.fdopen(fd, "wb") as out: out.write(zf.read(member)) with xr.open_dataset(target) as raw: - var_datasets.append(raw.load().rename({"ws100": new_name})) + var_datasets.append(raw.load().rename({raw_name: new_name})) finally: os.remove(target) if len(var_datasets) == 1: @@ -319,4 +340,41 @@ def retrieve_once(t: dict[str, Any]) -> xr.Dataset: result = xr.concat(datasets, dim="time").sel(time=coords["time"]) assert isinstance(result, xr.Dataset) + + # The delivered coefficients do not depend on the requested year or month, + # so a single request for the first month is sufficient. + alpha = retrieve_data( + variable=[POWER_LAW_VARIABLE], + year=[str(time[0].year)], + month=[time[0].strftime("%m")], + day=["01"], + **retrieval_params, + ) + result[POWER_LAW_NAME] = _climatology_to_time(alpha[POWER_LAW_NAME], time) return result + + +def _climatology_to_time(clim: xr.DataArray, time: pd.DatetimeIndex) -> xr.DataArray: + """ + Map a month-by-hour climatology onto the timestamps of a cutout. + + Parameters + ---------- + clim : xr.DataArray + Climatology with a 'time' dimension whose timestamps carry the month + and hour of day they apply to; the year and day are irrelevant. + time : pd.DatetimeIndex + Target timestamps. + + Returns + ------- + xr.DataArray + Climatology values for every target timestamp. + """ + clim_time = clim.indexes["time"] + positions = pd.Index(clim_time.month * 24 + clim_time.hour).get_indexer( + time.month * 24 + time.hour + ) + if (positions < 0).any(): + raise ValueError("Power law coefficients do not cover all months and hours.") + return clim.isel(time=positions).assign_coords(time=time) From 146bf5509918ceaaf311b97b6756bf6645855c02 Mon Sep 17 00:00:00 2001 From: jome1 Date: Tue, 13 Oct 2026 16:15:52 +0200 Subject: [PATCH 3/4] singel features: roughness, wind_corrected, wind_corrected_extrapolated --- atlite/datasets/era5.py | 25 +++++++++++ atlite/datasets/era5_corrected.py | 75 +++++++++++++++++-------------- atlite/datasets/era5_edh.py | 25 +++++++++++ 3 files changed, 91 insertions(+), 34 deletions(-) diff --git a/atlite/datasets/era5.py b/atlite/datasets/era5.py index 8c76f73d..5c6a875a 100644 --- a/atlite/datasets/era5.py +++ b/atlite/datasets/era5.py @@ -60,6 +60,7 @@ def nullcontext(): # noqa: D103 features = { "height": ["height"], "wind": ["wnd100m", "wnd_shear_exp", "wnd_azimuth", "roughness"], + "roughness": ["roughness"], "influx": [ "influx_toa", "influx_direct", @@ -199,6 +200,30 @@ def get_data_wind(retrieval_params: dict[str, Any]) -> xr.Dataset: return _process_wind(ds) +def get_data_roughness(retrieval_params: dict[str, Any]) -> xr.Dataset: + """ + Retrieve the forecast surface roughness from ERA5. + + Parameters + ---------- + retrieval_params : dict[str, Any] + CDS API retrieval parameters including area, time, and format. + + Returns + ------- + xr.Dataset + Dataset with variable: roughness. + """ + ds = retrieve_data(variable=["forecast_surface_roughness"], **retrieval_params) + ds = _rename_and_clean_coords(ds) + return ds.rename({"fsr": "roughness"}) + + +def sanitize_roughness(ds: xr.Dataset) -> xr.Dataset: + """Clip negative roughness values to a minimum of 2e-4.""" + return sanitize_wind(ds) + + def sanitize_wind(ds: xr.Dataset) -> xr.Dataset: """ Clip negative roughness values to a minimum of 2e-4. diff --git a/atlite/datasets/era5_corrected.py b/atlite/datasets/era5_corrected.py index eab45f34..843d9d57 100644 --- a/atlite/datasets/era5_corrected.py +++ b/atlite/datasets/era5_corrected.py @@ -8,19 +8,20 @@ ``sis-energy-global-reanalysis``, see https://cds.climate.copernicus.eu/datasets/sis-energy-global-reanalysis -Only the 100 m wind speed is available, in two variants which are both stored -in the cutout together with the power law coefficients (feature -``wind_corrected``): +Only the 100 m wind speed is available, in two variants, each selectable as a +separate feature: -* ``wnd100m_corrected``: bias-corrected 100 m wind speed +* ``wind_corrected``: ``wnd100m_corrected``, bias-corrected 100 m wind speed (CDS variable ``100m_wind_speed``) -* ``wnd100m_corrected_extrapolated``: wind speed extrapolated to 100 m with a - power law from the bias-corrected lower level wind speed - (CDS variable ``100m_wind_speed_extrapolated``) -* ``wnd_shear_exp_corrected``: power law coefficient (alpha) matching the - corrected wind speed (CDS variable ``wind_power_law_coefficients``). It is a - month-by-hour climatology, repeated for every year of the cutout, and is - requested separately from the wind speeds. +* ``wind_corrected_extrapolated``: ``wnd100m_corrected_extrapolated``, wind + speed extrapolated to 100 m with a power law from the bias-corrected lower + level wind speed (CDS variable ``100m_wind_speed_extrapolated``) + +The feature ``wind_shear_exp_corrected`` provides ``wnd_shear_exp_corrected``, +the power law coefficient (alpha) matching the corrected wind speed (CDS +variable ``wind_power_law_coefficients``). It is a month-by-hour climatology, +repeated for every year of the cutout, and is requested separately from the +wind speeds. Hub heights other than 100 m need either the power law coefficient (use ``method="power"`` and assign ``wnd_shear_exp_corrected`` to ``wnd_shear_exp``) @@ -68,11 +69,9 @@ crs = 4326 features = { - "wind_corrected": [ - "wnd100m_corrected", - "wnd100m_corrected_extrapolated", - "wnd_shear_exp_corrected", - ] + "wind_corrected": ["wnd100m_corrected"], + "wind_corrected_extrapolated": ["wnd100m_corrected_extrapolated"], + "wind_shear_exp_corrected": ["wnd_shear_exp_corrected"], } DATASET = "sis-energy-global-reanalysis" @@ -81,6 +80,11 @@ "100m_wind_speed": ("_mbc_", "wnd100m_corrected"), "100m_wind_speed_extrapolated": ("_pwl_", "wnd100m_corrected_extrapolated"), } +# Feature -> CDS wind variable requested for it +_FEATURE_WIND_VARIABLE = { + "wind_corrected": "100m_wind_speed", + "wind_corrected_extrapolated": "100m_wind_speed_extrapolated", +} # Power law coefficients (alpha) are a month-by-hour climatology, not a time # series, and are requested separately to keep the request size small. POWER_LAW_VARIABLE = "wind_power_law_coefficients" @@ -252,7 +256,8 @@ def get_data( cutout : Cutout Cutout object defining the spatial and temporal extent. feature : str - Feature to retrieve. Only 'wind_corrected' is supported. + Feature to retrieve: 'wind_corrected', 'wind_corrected_extrapolated' + or 'wind_shear_exp_corrected'. tmpdir : PathLike Directory for temporary download files. lock : SerializableLock or None, optional @@ -268,7 +273,8 @@ def get_data( Notes ----- Timeframes longer than 12 months are rejected because the CDS API request - would exceed its size limit. Each wind variable is downloaded separately. + would exceed its size limit. The wind speed and the power law coefficients + are downloaded in separate requests. Returns ------- @@ -326,12 +332,24 @@ def get_data( for year in time.year.unique() ] + if feature == "wind_shear_exp_corrected": + # The delivered coefficients do not depend on the requested year or + # month, so a single request for the first month is sufficient. + alpha = retrieve_data( + variable=[POWER_LAW_VARIABLE], + year=[str(time[0].year)], + month=[time[0].strftime("%m")], + day=["01"], + **retrieval_params, + ) + return xr.Dataset( + {POWER_LAW_NAME: _climatology_to_time(alpha[POWER_LAW_NAME], time)} + ) + def retrieve_once(t: dict[str, Any]) -> xr.Dataset: - datasets = [ - retrieve_data(variable=[variable], **retrieval_params, **t) - for variable in WIND_VARIABLES - ] - return xr.merge(datasets) + return retrieve_data( + variable=[_FEATURE_WIND_VARIABLE[feature]], **retrieval_params, **t + ) if concurrent_requests: datasets = compute(*[delayed(retrieve_once)(c) for c in time_chunks]) @@ -340,17 +358,6 @@ def retrieve_once(t: dict[str, Any]) -> xr.Dataset: result = xr.concat(datasets, dim="time").sel(time=coords["time"]) assert isinstance(result, xr.Dataset) - - # The delivered coefficients do not depend on the requested year or month, - # so a single request for the first month is sufficient. - alpha = retrieve_data( - variable=[POWER_LAW_VARIABLE], - year=[str(time[0].year)], - month=[time[0].strftime("%m")], - day=["01"], - **retrieval_params, - ) - result[POWER_LAW_NAME] = _climatology_to_time(alpha[POWER_LAW_NAME], time) return result diff --git a/atlite/datasets/era5_edh.py b/atlite/datasets/era5_edh.py index 848cbb16..76e2c9dd 100644 --- a/atlite/datasets/era5_edh.py +++ b/atlite/datasets/era5_edh.py @@ -57,6 +57,7 @@ features: dict[str, list[str]] = { "height": ["height"], "wind": ["wnd100m", "wnd_shear_exp", "wnd_azimuth", "roughness"], + "roughness": ["roughness"], "influx": [ "influx_toa", "influx_direct", @@ -74,6 +75,7 @@ _FEATURE_VARS: dict[str, list[str]] = { "wind": ["u10", "v10", "u100", "v100", "fsr"], + "roughness": ["fsr"], "influx": ["ssrd", "ssr", "fdir", "tisr"], "temperature": ["t2m", "d2m", "stl4"], "runoff": ["ro"], @@ -350,6 +352,27 @@ def get_data_wind(cutout: Cutout) -> xr.Dataset: return ds +def get_data_roughness(cutout: Cutout) -> xr.Dataset: + """ + Retrieve the forecast surface roughness only. + + Parameters + ---------- + cutout : atlite.Cutout + Cutout defining the target grid and time range. + + Returns + ------- + xarray.Dataset + Dataset containing ``roughness``. + """ + ds = _load_feature(cutout, "roughness").rename({"fsr": "roughness"}) + ds["roughness"] = ds["roughness"].assign_attrs( + units="m", long_name="Forecast surface roughness" + ) + return ds + + def get_data_influx(cutout: Cutout) -> xr.Dataset: """ Retrieve and prepare solar influx variables. @@ -434,6 +457,7 @@ def get_data_height(cutout: Cutout) -> xr.Dataset: _HANDLERS: dict[str, Handler] = { "wind": get_data_wind, + "roughness": get_data_roughness, "influx": get_data_influx, "temperature": get_data_temperature, "runoff": get_data_runoff, @@ -442,6 +466,7 @@ def get_data_height(cutout: Cutout) -> xr.Dataset: _SANITIZERS: dict[str, Sanitizer] = { "wind": sanitize_wind, + "roughness": sanitize_wind, "influx": sanitize_influx, "runoff": sanitize_runoff, } From 5b12b3738e49ac1db1dac55934b8067858eb24d4 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 13 Oct 2026 14:58:26 +0000 Subject: [PATCH 4/4] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- atlite/datasets/era5_corrected.py | 15 +++++++++------ atlite/wind.py | 4 +--- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/atlite/datasets/era5_corrected.py b/atlite/datasets/era5_corrected.py index 843d9d57..95dad7f8 100644 --- a/atlite/datasets/era5_corrected.py +++ b/atlite/datasets/era5_corrected.py @@ -45,6 +45,7 @@ import logging import os +import pathlib import tempfile import zipfile from contextlib import nullcontext @@ -223,7 +224,7 @@ def retrieve_data( with xr.open_dataset(target) as raw: var_datasets.append(raw.load().rename({raw_name: new_name})) finally: - os.remove(target) + pathlib.Path(target).unlink() if len(var_datasets) == 1: datasets.append(var_datasets[0]) else: @@ -231,7 +232,7 @@ def retrieve_data( combined = combined.sortby("time") datasets.append(combined) finally: - os.remove(archive) + pathlib.Path(archive).unlink() ds = xr.merge(datasets) ds = _rename_and_clean_coords(ds) @@ -282,7 +283,9 @@ def get_data( Data for the requested feature, aligned to cutout coordinates. """ if feature not in features: - raise NotImplementedError(f"Feature {feature!r} not supported by era5-corrected") + raise NotImplementedError( + f"Feature {feature!r} not supported by era5-corrected" + ) if (cutout.dx, cutout.dy) != (0.25, 0.25): raise NotImplementedError( "era5-corrected only supports the native 0.25°×0.25° grid." @@ -342,9 +345,9 @@ def get_data( day=["01"], **retrieval_params, ) - return xr.Dataset( - {POWER_LAW_NAME: _climatology_to_time(alpha[POWER_LAW_NAME], time)} - ) + return xr.Dataset({ + POWER_LAW_NAME: _climatology_to_time(alpha[POWER_LAW_NAME], time) + }) def retrieve_once(t: dict[str, Any]) -> xr.Dataset: return retrieve_data( diff --git a/atlite/wind.py b/atlite/wind.py index f0b4bef5..f40aa058 100644 --- a/atlite/wind.py +++ b/atlite/wind.py @@ -80,9 +80,7 @@ def extrapolate_wind_speed( if from_height is None: heights: NDArray = np.asarray([ - int(m.group(1)) - for s in ds - if (m := re.fullmatch(r"wnd(\d+)m", str(s))) + int(m.group(1)) for s in ds if (m := re.fullmatch(r"wnd(\d+)m", str(s))) ]) if len(heights) == 0: