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

atlite/convert.py

@@ -478,19 +478,28 @@
 
 
 # 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 @@
     return da
 
 
-def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params):
+def wind(
+    cutout,
+    turbine,
+    smooth=False,
+    add_cutout_windspeed=False,
+    real_long_run_average_windspeed=None,
+    real_long_run_average_height=100,
+    **params,
+):
     """
     Generate wind generation time-series.
 
@@ -529,6 +546,10 @@
         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.
+    real_long_run_average_speed : Path or rasterio dataset or None
+        Raster dataset with wind speeds to bias correct long run average wind speeds
+    real_long_run_average_height : int = 100
+        Height in meters of provided real_long_average_speed dataset
 
     Note
     ----
@@ -549,8 +570,24 @@
     if smooth:
         turbine = windturbine_smooth(turbine, params=smooth)
 
+    windspeed_bias_correction = None
+    from_height = None
+
+    if real_long_run_average_windspeed is not None:
+        windspeed_bias_correction = windm.calculate_windspeed_bias_correction(
+            cutout.data,
+            real_long_run_average_windspeed,
+            real_long_run_average_height,
+            crs=cutout.crs,
+        )
+        from_height = real_long_run_average_height
+
     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=from_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
@@ -45,6 +46,7 @@
 features = {
     "height": ["height"],
     "wind": ["wnd100m", "wnd_azimuth", "roughness"],
+    "wind_lra": ["wnd100m_lra"],
     "influx": [
         "influx_toa",
         "influx_direct",
@@ -58,6 +60,7 @@
 }
 
 static_features = {"height"}
+longrunaverage_features = {"wind_lra"}
 
 
 def _add_height(ds):
@@ -101,6 +104,32 @@
     return ds
 
 
+def get_data_wind_lra(retrieval_params):
+    """
+    Get long run average wind data for given retrieval parameters.
+    """
+    ds = retrieve_data(
+        variable=[
+            "100m_u_component_of_wind",
+            "100m_v_component_of_wind",
+        ],
+        **retrieval_params,
+    )
+    ds = _rename_and_clean_coords(ds)
+
+    ds["wnd100m_lra"] = (
+        sqrt(ds["u100"] ** 2 + ds["v100"] ** 2)
+        .mean("date")
+        .assign_attrs(
+            units=ds["u100"].attrs["units"],
+            long_name="100 metre wind speed as long run average",
+        )
+    )
+    ds = ds.drop_vars(["u100", "v100"])
+
+    return ds
+
+
 def get_data_wind(retrieval_params):
     """
     Get wind data for given retrieval parameters.
@@ -255,7 +284,7 @@
     return [y1, x0, y0, x1]
 
 
-def retrieval_times(coords, static=False, monthly_requests=False):
+def retrieval_times(coords, static=False, monthly_requests=False, longrunaverage=False):
     """
     Get list of retrieval cdsapi arguments for time dimension in coordinates.
 
@@ -287,6 +316,12 @@
             "day": str(time[0].day),
             "time": time[0].strftime("%H:00"),
         }
+    elif longrunaverage:
+        return {
+            "year": [str(y) for y in range(1980, datetime.date.today().year)],
+            "month": [f"{m:02}" for m in range(1, 12 + 1)],
+            "time": ["00:00"],
+        }
 
     # Prepare request for all months and years
     times = []
@@ -323,7 +358,7 @@
         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 +375,7 @@
     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()
@@ -418,7 +453,8 @@
     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,14 +467,18 @@
 
     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)
         return ds
 
     if feature in static_features:
         return retrieve_once(retrieval_times(coords, static=True)).squeeze()
+    elif feature in longrunaverage_features:
+        retrieval_params["dataset"] += "-monthly-means"
+        retrieval_params["product_type"] = "monthly_averaged_reanalysis"
+        return retrieve_once(retrieval_times(coords, longrunaverage=True))
 
     time_chunks = retrieval_times(coords, monthly_requests=monthly_requests)
     if concurrent_requests:

atlite/wind.py

@@ -7,8 +7,13 @@
 
 import logging
 import re
+from pathlib import Path
 
 import numpy as np
+import rasterio as rio
+import xarray as xr
+
+from .gis import _as_transform
 
 logger = logging.getLogger(__name__)
 
@@ -51,14 +56,11 @@
     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 +69,46 @@
 
     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(ds, real_lra, lra_height, crs):
+    data_lra = ds[f"wnd{lra_height}m_lra"]
+
+    if isinstance(real_lra, (str, Path)):
+        real_lra = rio.open(real_lra)
+
+    if isinstance(real_lra, rio.DatasetReader):
+        real_lra = rio.band(real_lra, 1)
+
+    if isinstance(real_lra, rio.Band):
+        real_lra, transform = rio.warp.reproject(
+            real_lra,
+            np.empty(data_lra.shape),
+            dst_crs=crs,
+            dst_transform=_as_transform(
+                x=data_lra.indexes["x"], y=data_lra.indexes["y"]
+            ),
+            resampling=rio.enums.Resampling.average,
+        )
+
+        real_lra = xr.DataArray(
+            real_lra, [data_lra.indexes["y"], data_lra.indexes["x"]]
+        )
 
-    return wnd_spd.rename(to_name)
+    return real_lra / data_lra