feat(wind): Bias correct wind speeds based on scaling to a known average

Adds a new function :py:func:`atlite.datasets.era5.retrieve_average_windspeed`
to retrieve average windspeeds. This dataset function is called by
:py:func:`atlite.wind.calculate_windspeed_bias_correction` to derive a bias
correction which can be passed to the default wind generation.

Example usage:

```python
windspeed_bias_correction = atlite.wind.calculate_windspeed_bias_correction(
    cutout,
    real_average="gwa3_250_windspeed_100m.tif",
    height=100,
)
cutout.wind(
    atlite.windturbines.Enercon_E101_3000kW,
    windspeed_bias_correction=windspeed_bias_correction,
    windspeed_height=100,
)
```

atlite/convert.py

@@ -478,19 +478,28 @@ def solar_thermal(
 
 
 # wind
-def convert_wind(ds, turbine):
+def convert_wind(
+    ds: xr.Dataset, turbine, windspeed_bias_correction=None, from_height=None
+):
     """
     Convert wind speeds for turbine to wind energy generation.
     """
     V, POW, hub_height, P = itemgetter("V", "POW", "hub_height", "P")(turbine)
 
-    wnd_hub = windm.extrapolate_wind_speed(ds, to_height=hub_height)
+    if windspeed_bias_correction is not None:
+        ds = ds.assign(
+            {f"wnd{from_height}m": ds[f"wnd{from_height}m"] * windspeed_bias_correction}
+        )
+
+    wnd_hub = windm.extrapolate_wind_speed(
+        ds, to_height=hub_height, from_height=from_height
+    )
 
-    def _interpolate(da):
+    def apply_power_curve(da):
         return np.interp(da, V, POW / P)
 
     da = xr.apply_ufunc(
-        _interpolate,
+        apply_power_curve,
         wnd_hub,
         input_core_dims=[[]],
         output_core_dims=[[]],
@@ -503,7 +512,15 @@ def _interpolate(da):
     return da
 
 
-def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params):
+def wind(
+    cutout,
+    turbine: str | dict,
+    smooth: bool | dict = False,
+    add_cutout_windspeed: bool = False,
+    windspeed_bias_correction: xr.DataArray | None = None,
+    windspeed_height: int | None = None,
+    **params,
+):
     """
     Generate wind generation time-series.
 
@@ -513,22 +530,32 @@ def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params):
     Parameters
     ----------
     turbine : str or dict
-        A turbineconfig dictionary with the keys 'hub_height' for the
-        hub height and 'V', 'POW' defining the power curve.
-        Alternatively a str refering to a local or remote turbine configuration
-        as accepted by atlite.resource.get_windturbineconfig(). Locally stored turbine
-        configurations can also be modified with this function. E.g. to setup a different hub
-        height from the one used in the yaml file,one would write
-                "turbine=get_windturbineconfig(“NREL_ReferenceTurbine_5MW_offshore”)|{“hub_height”:120}"
+        A turbineconfig dictionary with the keys 'hub_height' for the hub height
+        and 'V', 'POW' defining the power curve. Alternatively a str refering to
+        a local or remote turbine configuration as accepted by
+        atlite.resource.get_windturbineconfig(). Locally stored turbine
+        configurations can also be modified with this function. E.g. to setup a
+        different hub height from the one used in the yaml file, one would write
+        >>> turbine = (
+        >>>     get_windturbineconfig("NREL_ReferenceTurbine_5MW_offshore")
+        >>>     | {"hub_height":120}
+        >>> )
     smooth : bool or dict
-        If True smooth power curve with a gaussian kernel as
-        determined for the Danish wind fleet to Delta_v = 1.27 and
-        sigma = 2.29. A dict allows to tune these values.
+        If True smooth power curve with a gaussian kernel as determined for the
+        Danish wind fleet to Delta_v = 1.27 and sigma = 2.29. A dict allows to
+        tune these values.
     add_cutout_windspeed : bool
-        If True and in case the power curve does not end with a zero, will add zero power
-        output at the highest wind speed in the power curve. If False, a warning will be
-        raised if the power curve does not have a cut-out wind speed. The default is
-        False.
+        If True and in case the power curve does not end with a zero, will add
+        zero power output at the highest wind speed in the power curve. If
+        False, a warning will be raised if the power curve does not have a
+        cut-out wind speed. The default is False.
+    windspeed_bias_correction : DataArray, optional
+        Correction factor that is applied to the windspeed at
+        `windspeed_height`. Such a correction factor can be calculated using
+        :py:func:`atlite.wind.calculate_windspeed_bias_correction` with a raster
+        dataset of mean wind speeds.
+    windspeed_height : int, optional
+        Height in meters of windspeed data from which to extrapolate
 
     Note
     ----
@@ -541,7 +568,7 @@ def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params):
     1074 – 1088. doi:10.1016/j.energy.2015.09.071
 
     """
-    if isinstance(turbine, (str, Path, dict)):
+    if isinstance(turbine, str | Path | dict):
         turbine = get_windturbineconfig(
             turbine, add_cutout_windspeed=add_cutout_windspeed
         )
@@ -550,7 +577,11 @@ def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params):
         turbine = windturbine_smooth(turbine, params=smooth)
 
     return cutout.convert_and_aggregate(
-        convert_func=convert_wind, turbine=turbine, **params
+        convert_func=convert_wind,
+        turbine=turbine,
+        windspeed_bias_correction=windspeed_bias_correction,
+        from_height=windspeed_height,
+        **params,
     )
 
 

atlite/datasets/era5.py

@@ -8,6 +8,7 @@
 https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation
 """
 
+import datetime
 import logging
 import os
 import warnings
@@ -323,7 +324,7 @@ def noisy_unlink(path):
         logger.error(f"Unable to delete file {path}, as it is still in use.")
 
 
-def retrieve_data(product, chunks=None, tmpdir=None, lock=None, **updates):
+def retrieve_data(dataset, chunks=None, tmpdir=None, lock=None, **updates):
     """
     Download data like ERA5 from the Climate Data Store (CDS).
 
@@ -340,7 +341,7 @@ def retrieve_data(product, chunks=None, tmpdir=None, lock=None, **updates):
     client = cdsapi.Client(
         info_callback=logger.debug, debug=logging.DEBUG >= logging.root.level
     )
-    result = client.retrieve(product, request)
+    result = client.retrieve(dataset, request)
 
     if lock is None:
         lock = nullcontext()
@@ -359,7 +360,7 @@ def retrieve_data(product, chunks=None, tmpdir=None, lock=None, **updates):
     ds = xr.open_dataset(target, chunks=chunks or {})
     if tmpdir is None:
         logger.debug(f"Adding finalizer for {target}")
-        weakref.finalize(ds._file_obj._manager, noisy_unlink, target)
+        weakref.finalize(ds._close.__self__.ds, noisy_unlink, target)
 
     # Remove default encoding we get from CDSAPI, which can lead to NaN values after loading with subsequent
     # saving due to how xarray handles netcdf compression (only float encoded as short int seem affected)
@@ -373,6 +374,55 @@ def retrieve_data(product, chunks=None, tmpdir=None, lock=None, **updates):
     return ds
 
 
+def retrieve_windspeed_average(cutout, height, first_year=1980, last_year=None):
+    """
+    Retrieve average windspeed from `first_year` to `last_year`
+
+    Parameters
+    ----------
+    cutout : atlite.Cutout
+        Cutout for which to retrieve windspeeds from CDS
+    height : int
+        Height of windspeeds (ERA5 typically knows about 10m, 100m, 150m?)
+    first_year : int
+        First year to take into account
+    last_year : int, optional
+        Last year to take into account (if omitted takes the previous year)
+
+    Returns
+    -------
+    DataArray
+        Mean windspeed at cutout dimension
+    """
+    if last_year is None:
+        last_year = datetime.date.today().year - 1
+
+    ds = retrieve_data(
+        "reanalysis-era5-single-levels-monthly-means",
+        chunks=cutout.chunks,
+        product_type="monthly_averaged_reanalysis",
+        variable=[
+            f"{height}m_u_component_of_wind",
+            f"{height}m_v_component_of_wind",
+        ],
+        area=_area(cutout.coords),
+        grid=[cutout.dx, cutout.dy],
+        year=[str(y) for y in range(first_year, last_year + 1)],
+        month=[f"{m:02}" for m in range(1, 12 + 1)],
+        time=["00:00"],
+    )
+    ds = _rename_and_clean_coords(ds)
+
+    return (
+        sqrt(ds[f"u{height}"] ** 2 + ds[f"v{height}"] ** 2)
+        .mean("date")
+        .assign_attrs(
+            units=ds[f"u{height}"].attrs["units"],
+            long_name=f"{height} metre wind speed as long run average",
+        )
+    )
+
+
 def get_data(
     cutout,
     feature,
@@ -418,7 +468,8 @@ def get_data(
     sanitize = creation_parameters.get("sanitize", True)
 
     retrieval_params = {
-        "product": "reanalysis-era5-single-levels",
+        "dataset": "reanalysis-era5-single-levels",
+        "product_type": "reanalysis",
         "area": _area(coords),
         "chunks": cutout.chunks,
         "grid": [cutout.dx, cutout.dy],
@@ -431,7 +482,7 @@ def get_data(
 
     logger.info(f"Requesting data for feature {feature}...")
 
-    def retrieve_once(time):
+    def retrieve_once(time, longrunaverage=False):
         ds = func({**retrieval_params, **time})
         if sanitize and sanitize_func is not None:
             ds = sanitize_func(ds)

atlite/wind.py

@@ -7,13 +7,20 @@
 
 import logging
 import re
+from pathlib import Path
 
 import numpy as np
+import rasterio as rio
+import xarray as xr
+
+from . import datasets
 
 logger = logging.getLogger(__name__)
 
 
-def extrapolate_wind_speed(ds, to_height, from_height=None):
+def extrapolate_wind_speed(
+    ds: xr.Dataset, to_height: int | float, from_height: int | None = None
+):
     """
     Extrapolate the wind speed from a given height above ground to another.
 
@@ -28,8 +35,7 @@ def extrapolate_wind_speed(ds, to_height, from_height=None):
         Dataset containing the wind speed time-series at 'from_height' with key
         'wnd{height:d}m' and the surface orography with key 'roughness' at the
         geographic locations of the wind speeds.
-    from_height : int
-        (Optional)
+    from_height : int, optional
         Height (m) from which the wind speeds are interpolated to 'to_height'.
         If not provided, the closest height to 'to_height' is selected.
     to_height : int|float
@@ -51,14 +57,11 @@ def extrapolate_wind_speed(ds, to_height, from_height=None):
     Retrieved 2019-02-15.
 
     """
-    # Fast lane
-    to_name = f"wnd{int(to_height):0d}m"
-    if to_name in ds:
-        return ds[to_name]
-
     if from_height is None:
         # Determine closest height to to_name
-        heights = np.asarray([int(s[3:-1]) for s in ds if re.match(r"wnd\d+m", s)])
+        heights = np.asarray(
+            [int(m.group(1)) for s in ds if (m := re.match(r"wnd(\d+)m", s))]
+        )
 
         if len(heights) == 0:
             raise AssertionError("Wind speed is not in dataset")
@@ -67,17 +70,81 @@ def extrapolate_wind_speed(ds, to_height, from_height=None):
 
     from_name = f"wnd{int(from_height):0d}m"
 
+    # Fast lane
+    if from_height == to_height:
+        return ds[from_name]
+
     # Wind speed extrapolation
     wnd_spd = ds[from_name] * (
         np.log(to_height / ds["roughness"]) / np.log(from_height / ds["roughness"])
     )
 
-    wnd_spd.attrs.update(
+    return wnd_spd.assign_attrs(
         {
             "long name": f"extrapolated {to_height} m wind speed using logarithmic "
             f"method with roughness and {from_height} m wind speed",
             "units": "m s**-1",
         }
-    )
+    ).rename(f"wnd{to_height}m")
+
+
+def calculate_windspeed_bias_correction(
+    cutout, real_average: str | rio.DatasetReader, height: int = 100
+):
+    """
+    Derive a bias correction factor for windspeed at lra_height
+
+    Regrids the raster dataset in real_average to cutout grid, retrieves the average
+    windspeed from the first dataset that offers
+    :py:func:`retrieve_longrunaverage_windspeed` (only ERA5, currently).
+
+    Parameters
+    ----------
+    cutout : Cutout
+        Atlite cutout
+    real_average : Path or rasterio.Dataset
+        Raster dataset with wind speeds to bias correct average wind speeds
+    height : int
+        Height in meters at which average windspeeds are provided
 
-    return wnd_spd.rename(to_name)
+    Returns
+    -------
+    DataArray
+        Ratio between windspeeds in `real_average` to those of average windspeeds in
+        dataset.
+    """
+    if isinstance(real_average, str | Path):
+        real_average = rio.open(real_average)
+
+    if isinstance(real_average, rio.DatasetReader):
+        real_average = rio.band(real_average, 1)
+
+    if isinstance(real_average, rio.Band):
+        real_average, transform = rio.warp.reproject(
+            real_average,
+            np.empty(cutout.shape),
+            dst_crs=cutout.crs,
+            dst_transform=cutout.transform,
+            dst_nodata=np.nan,
+            resampling=rio.enums.Resampling.average,
+        )
+
+        real_average = xr.DataArray(
+            real_average, [cutout.coords["y"], cutout.coords["x"]]
+        )
+
+    for module in np.atleast_1d(cutout.module):
+        retrieve_windspeed_average = getattr(
+            getattr(datasets, module), "retrieve_windspeed_average"
+        )
+        if retrieve_windspeed_average is not None:
+            break
+    else:
+        raise AssertionError(
+            "None of the datasets modules define retrieve_windspeed_average"
+        )
+
+    logger.info("Retrieving average windspeeds at %d from module %s", height, module)
+    data_average = retrieve_windspeed_average(cutout, height)
+
+    return real_average / data_average