feat(wind): Add power law interpolation method for wind conversion (#402)

* feat(wind): Add power law interpolation method for wind conversion

Refer to #231 for context.

* fix wind docstring

* Update atlite/convert.py

Co-authored-by: Lukas Trippe <[email protected]>

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Add typing for wind conversion

* fix(wind): Raise RuntimeError from extrapolate_wind_speed with an explanation

If auxiliary data is missing from the cutout.

---------

Co-authored-by: Jonas Hoersch <[email protected]>
Co-authored-by: Lukas Trippe <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: 3 people

atlite/convert.py

@@ -5,11 +5,14 @@
 All functions for converting weather data into energy system model data.
 """
 
+from __future__ import annotations
+
 import datetime as dt
 import logging
 from collections import namedtuple
 from operator import itemgetter
 from pathlib import Path
+from typing import TYPE_CHECKING
 
 import geopandas as gpd
 import numpy as np
@@ -39,6 +42,11 @@
 
 logger = logging.getLogger(__name__)
 
+if TYPE_CHECKING:
+    from typing import Literal
+
+    from atlite.resource import TurbineConfig
+
 
 def convert_and_aggregate(
     cutout,
@@ -478,19 +486,25 @@ def solar_thermal(
 
 
 # wind
-def convert_wind(ds, turbine):
+def convert_wind(
+    ds: xr.Dataset,
+    turbine: TurbineConfig,
+    interpolation_method: Literal["logarithmic", "power"],
+) -> xr.DataArray:
     """
     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)
+    wnd_hub = windm.extrapolate_wind_speed(
+        ds, to_height=hub_height, method=interpolation_method
+    )
 
-    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,12 +517,19 @@ def _interpolate(da):
     return da
 
 
-def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params):
+def wind(
+    cutout,
+    turbine: str | Path | dict,
+    smooth: bool | dict = False,
+    add_cutout_windspeed: bool = False,
+    interpolation_method: Literal["logarithmic", "power"] = "logarithmic",
+    **params,
+) -> xr.DataArray:
     """
     Generate wind generation time-series.
 
-    Extrapolates 10m wind speed with monthly surface roughness to hub
-    height and evaluates the power curve.
+    Extrapolates wind speed to hub height (using logarithmic or power law) and
+    evaluates the power curve.
 
     Parameters
     ----------
@@ -529,6 +550,9 @@ def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params):
         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.
+    interpolation_method : {"logarithmic", "power"}
+        Law to interpolate wind speed to turbine hub height. Refer to
+        :py:func:`atlite.wind.extrapolate_wind_speed`.
 
     Note
     ----
@@ -537,20 +561,20 @@ def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params):
 
     References
     ----------
-    [1] Andresen G B, Søndergaard A A and Greiner M 2015 Energy 93, Part 1
-    1074 – 1088. doi:10.1016/j.energy.2015.09.071
+    .. [1] Andresen G B, Søndergaard A A and Greiner M 2015 Energy 93, Part 1
+       1074 – 1088. doi:10.1016/j.energy.2015.09.071
 
     """
-    if isinstance(turbine, (str, Path, dict)):
-        turbine = get_windturbineconfig(
-            turbine, add_cutout_windspeed=add_cutout_windspeed
-        )
+    turbine = get_windturbineconfig(turbine, add_cutout_windspeed=add_cutout_windspeed)
 
     if smooth:
         turbine = windturbine_smooth(turbine, params=smooth)
 
     return cutout.convert_and_aggregate(
-        convert_func=convert_wind, turbine=turbine, **params
+        convert_func=convert_wind,
+        turbine=turbine,
+        interpolation_method=interpolation_method,
+        **params,
     )
 
 

atlite/datasets/era5.py

@@ -44,7 +44,7 @@ def nullcontext():
 
 features = {
     "height": ["height"],
-    "wind": ["wnd100m", "wnd_azimuth", "roughness"],
+    "wind": ["wnd100m", "wnd_shear_exp", "wnd_azimuth", "roughness"],
     "influx": [
         "influx_toa",
         "influx_direct",
@@ -107,6 +107,8 @@ def get_data_wind(retrieval_params):
     """
     ds = retrieve_data(
         variable=[
+            "10m_u_component_of_wind",
+            "10m_v_component_of_wind",
             "100m_u_component_of_wind",
             "100m_v_component_of_wind",
             "forecast_surface_roughness",
@@ -115,13 +117,19 @@ def get_data_wind(retrieval_params):
     )
     ds = _rename_and_clean_coords(ds)
 
-    ds["wnd100m"] = sqrt(ds["u100"] ** 2 + ds["v100"] ** 2).assign_attrs(
-        units=ds["u100"].attrs["units"], long_name="100 metre wind speed"
-    )
+    for h in [10, 100]:
+        ds[f"wnd{h}m"] = sqrt(ds[f"u{h}"] ** 2 + ds[f"v{h}"] ** 2).assign_attrs(
+            units=ds[f"u{h}"].attrs["units"], long_name=f"{h} metre wind speed"
+        )
+    ds["wnd_shear_exp"] = (
+        np.log(ds["wnd10m"] / ds["wnd100m"]) / np.log(10 / 100)
+    ).assign_attrs(units="", long_name="wind shear exponent")
+
     # span the whole circle: 0 is north, π/2 is east, -π is south, 3π/2 is west
     azimuth = arctan2(ds["u100"], ds["v100"])
     ds["wnd_azimuth"] = azimuth.where(azimuth >= 0, azimuth + 2 * np.pi)
-    ds = ds.drop_vars(["u100", "v100"])
+
+    ds = ds.drop_vars(["u100", "v100", "u10", "v10", "wnd10m"])
     ds = ds.rename({"fsr": "roughness"})
 
     return ds

atlite/resource.py

@@ -6,12 +6,15 @@
 panel configurations.
 """
 
+from __future__ import annotations
+
 import json
 import logging
 import re
 import warnings
 from operator import itemgetter
 from pathlib import Path
+from typing import TYPE_CHECKING
 
 import numpy as np
 import pandas as pd
@@ -30,8 +33,24 @@
 SOLARPANEL_DIRECTORY = RESOURCE_DIRECTORY / "solarpanel"
 CSPINSTALLATION_DIRECTORY = RESOURCE_DIRECTORY / "cspinstallation"
 
+if TYPE_CHECKING:
+    from typing import TypedDict
+
+    from typing_extensions import NotRequired
+
+    class TurbineConfig(TypedDict):
+        V: np.ndarray
+        POW: np.ndarray
+        P: float
+        hub_height: float | int
+        name: NotRequired[str]
+        manufacturer: NotRequired[str]
+        source: NotRequired[str]
+
 
-def get_windturbineconfig(turbine, add_cutout_windspeed=False):
+def get_windturbineconfig(
+    turbine: str | Path | dict, add_cutout_windspeed: bool = False
+) -> TurbineConfig:
     """
     Load the wind 'turbine' configuration.
 
@@ -61,7 +80,7 @@ def get_windturbineconfig(turbine, add_cutout_windspeed=False):
         Config with details on the turbine
 
     """
-    assert isinstance(turbine, (str, Path, dict))
+    assert isinstance(turbine, str | Path | dict)
 
     if add_cutout_windspeed is False:
         msg = (
@@ -73,7 +92,7 @@ def get_windturbineconfig(turbine, add_cutout_windspeed=False):
     if isinstance(turbine, str) and turbine.startswith("oedb:"):
         conf = get_oedb_windturbineconfig(turbine[len("oedb:") :])
 
-    elif isinstance(turbine, (str, Path)):
+    elif isinstance(turbine, str | Path):
         if isinstance(turbine, str):
             turbine_path = windturbines[turbine.replace(".yaml", "")]
 
@@ -287,7 +306,9 @@ def _max_v_is_zero_pow(turbine):
     return np.any(turbine["POW"][turbine["V"] == turbine["V"].max()] == 0)
 
 
-def _validate_turbine_config_dict(turbine: dict, add_cutout_windspeed: bool):
+def _validate_turbine_config_dict(
+    turbine: dict, add_cutout_windspeed: bool
+) -> TurbineConfig:
     """
     Checks the turbine config dict format and power curve.
 
@@ -356,7 +377,9 @@ def _validate_turbine_config_dict(turbine: dict, add_cutout_windspeed: bool):
     return turbine
 
 
-def get_oedb_windturbineconfig(search=None, **search_params):
+def get_oedb_windturbineconfig(
+    search: int | str | None = None, **search_params
+) -> TurbineConfig:
     """
     Download a windturbine configuration from the OEDB database.
 

atlite/wind.py

@@ -5,51 +5,72 @@
 Functions for use in conjunction with wind data generation.
 """
 
+from __future__ import annotations
+
 import logging
 import re
+from typing import TYPE_CHECKING
 
 import numpy as np
+import xarray as xr
 
 logger = logging.getLogger(__name__)
 
 
-def extrapolate_wind_speed(ds, to_height, from_height=None):
+if TYPE_CHECKING:
+    from typing import Literal
+
+
+def extrapolate_wind_speed(
+    ds: xr.Dataset,
+    to_height: int | float,
+    from_height: int | None = None,
+    method: Literal["logarithmic", "power"] = "logarithmic",
+) -> xr.DataArray:
     """
     Extrapolate the wind speed from a given height above ground to another.
 
     If ds already contains a key refering to wind speeds at the desired to_height,
     no conversion is done and the wind speeds are directly returned.
 
-    Extrapolation of the wind speed follows the logarithmic law as desribed in [1].
+    Extrapolation of the wind speed can either use the "logarithmic" law as
+    described in [1]_ or the "power" law as described in [2]_. See also discussion
+    in GH issue: https://github.com/PyPSA/atlite/issues/231 .
 
     Parameters
     ----------
     ds : xarray.Dataset
         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)
-        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
         Height (m) to which the wind speeds are extrapolated to.
+    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.
+    method : {"logarithmic", "power"}
+        Method to use for extra/interpolating wind speeds
 
     Returns
     -------
     da : xarray.DataArray
         DataArray containing the extrapolated wind speeds. Name of the DataArray
         is 'wnd{to_height:d}'.
 
+    Raises
+    ------
+    RuntimeError
+        If the cutout is missing the data for the chosen `method`
+
     References
     ----------
-    [1] Equation (2) in Andresen, G. et al (2015): 'Validation of Danish wind
-    time series from a new global renewable energy atlas for energy system
-    analysis'.
-
-    [2] https://en.wikipedia.org/w/index.php?title=Roughness_length&oldid=862127433,
-    Retrieved 2019-02-15.
+    .. [1] Equation (2) in Andresen, G. et al (2015): 'Validation of Danish
+       wind time series from a new global renewable energy atlas for energy
+       system analysis'.
 
+    .. [2] Gualtieri, G. (2021): 'Reliability of ERA5 Reanalysis Data for
+       Wind Resource Assessment: A Comparison against Tall Towers'
+       https://doi.org/10.3390/en14144169 .
     """
     # Fast lane
     to_name = f"wnd{int(to_height):0d}m"
@@ -67,15 +88,40 @@ def extrapolate_wind_speed(ds, to_height, from_height=None):
 
     from_name = f"wnd{int(from_height):0d}m"
 
-    # Wind speed extrapolation
-    wnd_spd = ds[from_name] * (
-        np.log(to_height / ds["roughness"]) / np.log(from_height / ds["roughness"])
-    )
+    if method == "logarithmic":
+        try:
+            roughness = ds["roughness"]
+        except KeyError:
+            raise RuntimeError(
+                "The logarithmic interpolation method requires surface roughness (roughness);\n"
+                "make sure you choose a compatible dataset like ERA5"
+            )
+        wnd_spd = ds[from_name] * (
+            np.log(to_height / roughness) / np.log(from_height / roughness)
+        )
+        method_desc = "logarithmic method with roughness"
+    elif method == "power":
+        try:
+            wnd_shear_exp = ds["wnd_shear_exp"]
+        except KeyError:
+            raise RuntimeError(
+                "The power law interpolation method requires a wind shear exponent (wnd_shear_exp);\n"
+                "make sure you choose a compatible dataset like ERA5 and update your cutout"
+            )
+        wnd_spd = ds[from_name] * (to_height / from_height) ** wnd_shear_exp
+        method_desc = "power method with wind shear exponent"
+    else:
+        raise ValueError(
+            f"Interpolation method must be 'logarithmic' or 'power', "
+            f" but is: {method}"
+        )
 
     wnd_spd.attrs.update(
         {
-            "long name": f"extrapolated {to_height} m wind speed using logarithmic "
-            f"method with roughness and {from_height} m wind speed",
+            "long name": (
+                f"extrapolated {to_height} m wind speed using {method_desc} "
+                f" and {from_height} m wind speed"
+            ),
             "units": "m s**-1",
         }
     )