Fix power curves missing cutout wind speed (#316)

* warn for missing power curve cutoff, add option to modify power curve

* move power curve and turbine validation to new function

* add check for increasing V in the turbine config

* rename `cutoff` to `cutout_windspeed`

* specify release where change applies in FutureWarning

* add test and release note for add_cutout_windspeed

RELEASE_NOTES.rst

@@ -12,6 +12,11 @@ Release Notes
 .. Upcoming Release
 .. ================
 
+.. * Fix: the wind turbine power curve is checked for a missing cut-out wind speed and an option to add a
+..   cut-out wind speed at the end of the power curve is introduced. From the next release v0.2.13, adding
+..   a cut-out wind speed will be the default behavior (`GH #316 <https://github.com/PyPSA/atlite/pull/316>`_)
+
+
 Version 0.2.11
 ==============
 

atlite/convert.py

@@ -480,7 +480,7 @@ def _interpolate(da):
     return da
 
 
-def wind(cutout, turbine, smooth=False, **params):
+def wind(cutout, turbine, smooth=False, add_cutout_windspeed=False, **params):
     """
     Generate wind generation time-series.
 
@@ -501,6 +501,11 @@ def wind(cutout, turbine, smooth=False, **params):
         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.
 
     Note
     ----
@@ -512,8 +517,10 @@ def wind(cutout, turbine, smooth=False, **params):
     [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)):
-        turbine = get_windturbineconfig(turbine)
+    if isinstance(turbine, (str, Path, dict)):
+        turbine = get_windturbineconfig(
+            turbine, add_cutout_windspeed=add_cutout_windspeed
+        )
 
     if smooth:
         turbine = windturbine_smooth(turbine, params=smooth)

atlite/resource.py

@@ -11,6 +11,7 @@
 import json
 import logging
 import re
+import warnings
 from operator import itemgetter
 from pathlib import Path
 
@@ -32,13 +33,13 @@
 CSPINSTALLATION_DIRECTORY = RESOURCE_DIRECTORY / "cspinstallation"
 
 
-def get_windturbineconfig(turbine):
+def get_windturbineconfig(turbine, add_cutout_windspeed=False):
     """
     Load the wind 'turbine' configuration.
 
     Parameters
     ----------
-    turbine : str or pathlib.Path
+    turbine : str or pathlib.Path or dict
         if str:
             The name of a preshipped turbine from alite.resources.windturbine .
             Alternatively, if a str starting with 'oedb:<name>' is passed the Open
@@ -47,32 +48,52 @@ def get_windturbineconfig(turbine):
             `atlite.resource.get_oedb_windturbineconfig(...)`
         if `pathlib.Path` is provided the configuration is read from this local
             path instead
+        if dict:
+            a user provided config dict. Needs to have the keys "POW", "V", "P", and
+            "hub_height". Values for "POW" and "V" need to be list or np.ndarray with
+            equal length.
+    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.
 
     Returns
     ----------
     config : dict
         Config with details on the turbine
     """
-    assert isinstance(turbine, (str, Path))
+    assert isinstance(turbine, (str, Path, dict))
 
-    if isinstance(turbine, str) and turbine.startswith("oedb:"):
-        return get_oedb_windturbineconfig(turbine[len("oedb:") :])
-
-    elif isinstance(turbine, str):
-        turbine_path = windturbines[turbine.replace(".yaml", "")]
+    if add_cutout_windspeed is False:
+        msg = (
+            "'add_cutout_windspeed' for wind turbine\npower curves will default to "
+            "True in atlite relase v0.2.13."
+        )
+        warnings.warn(msg, FutureWarning)
 
-    elif isinstance(turbine, Path):
-        turbine_path = turbine
+    if isinstance(turbine, str) and turbine.startswith("oedb:"):
+        conf = get_oedb_windturbineconfig(turbine[len("oedb:") :])
+
+    elif isinstance(turbine, (str, Path)):
+        if isinstance(turbine, str):
+            turbine_path = windturbines[turbine.replace(".yaml", "")]
+
+        elif isinstance(turbine, Path):
+            turbine_path = turbine
+
+        with open(turbine_path, "r") as f:
+            conf = yaml.safe_load(f)
+            conf = dict(
+                V=np.array(conf["V"]),
+                POW=np.array(conf["POW"]),
+                hub_height=conf["HUB_HEIGHT"],
+                P=np.max(conf["POW"]),
+            )
 
-    with open(turbine_path, "r") as f:
-        conf = yaml.safe_load(f)
+    elif isinstance(turbine, dict):
+        conf = turbine
 
-    return dict(
-        V=np.array(conf["V"]),
-        POW=np.array(conf["POW"]),
-        hub_height=conf["HUB_HEIGHT"],
-        P=np.max(conf["POW"]),
-    )
+    return _validate_turbine_config_dict(conf, add_cutout_windspeed)
 
 
 def get_solarpanelconfig(panel):
@@ -260,6 +281,82 @@ def smooth(velocities, power):
     return turbine
 
 
+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):
+    """
+    Checks the turbine config dict format and power curve.
+
+    Parameters
+    ----------
+    turbine : dict
+        turbine configuration dict. Needs the keys "POW", "V", "P", and "hub_height".
+        Values for "V" and "POW" need to be list or np.ndarray.
+    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.
+
+    Returns
+    -------
+    dict
+        validated and potentially modified turbine config dict
+    """
+    if not all(key in turbine for key in ("POW", "V", "P", "hub_height")):
+        err_msg = (
+            "turbine config dict needs at least the following keys: ['POW', 'V', 'P', "
+            f"'hub_height']\nbut are currently: {list(turbine.keys())}"
+        )
+        raise ValueError(err_msg)
+
+    if not all(isinstance(turbine[p], (np.ndarray, list)) for p in ("POW", "V")):
+        err_msg = "turbine entries 'POW' and 'V' must be np.ndarray or list"
+        raise ValueError(err_msg)
+
+    # convert lists from user provided turbine dicts to numpy arrays
+    if any(isinstance(turbine[p], list) for p in ("POW", "V")):
+        turbine["V"] = np.array(turbine["V"])
+        turbine["POW"] = np.array(turbine["POW"])
+
+    if len(turbine["POW"]) != len(turbine["V"]):
+        err_msg = "turbine wind speed and power arrays do not have equal length."
+        raise ValueError(err_msg)
+
+    if not np.all(np.diff(turbine["V"]) >= 0):
+        # This check is not strict as it uses `>=` instead of `>` and thus allows equal
+        # wind speeds in the array. However, many power curves have two entries for the
+        # same wind speed at the cut-in and cut-out speeds which would make them fail if
+        # using `>` only.
+        err_msg = (
+            "wind speed 'V' in the turbine config dict is expected to be increasing, "
+            f"but is currently not in ascending order:\n{turbine['V']}"
+        )
+        raise ValueError(err_msg)
+
+    if add_cutout_windspeed is True and not _max_v_is_zero_pow(turbine):
+        turbine["V"] = np.pad(turbine["V"], (0, 1), "maximum")
+        turbine["POW"] = np.pad(turbine["POW"], (0, 1), "constant", constant_values=0)
+        logger.info(
+            (
+                "adding a cut-out wind speed to the turbine power curve at "
+                f"V={turbine['V'][-1]} m/s."
+            )
+        )
+
+    if not _max_v_is_zero_pow(turbine):
+        logger.warning(
+            (
+                "The power curve does not have a cut-out wind speed, i.e. the power"
+                " output corresponding to the\nhighest wind speed is not zero. You can"
+                " either change the power curve manually or set\n"
+                "'add_cutout_windspeed=True' in the Cutout.wind conversion method."
+            )
+        )
+    return turbine
+
+
 def get_oedb_windturbineconfig(search=None, **search_params):
     """
     Download a windturbine configuration from the OEDB database.

test/test_resource.py

@@ -11,7 +11,7 @@
 """
 import pytest
 
-from atlite.resource import get_oedb_windturbineconfig
+from atlite.resource import get_oedb_windturbineconfig, get_windturbineconfig
 
 
 def test_oedb_windturbineconfig():
@@ -23,3 +23,12 @@ def test_oedb_windturbineconfig():
 
     # test string search with param
     assert get_oedb_windturbineconfig("E-101/3500 E2", hub_height=99)
+
+
+@pytest.mark.parametrize("add_cutout, last_pow", [(True, 0.0), (False, 1.0)])
+def test_windturbineconfig_add_cutout(add_cutout, last_pow):
+    t = get_windturbineconfig(
+        {"V": [0, 25], "POW": [0.0, 1.0], "hub_height": 1.0, "P": 1.0},
+        add_cutout_windspeed=add_cutout,
+    )
+    assert t["POW"][-1] == last_pow