Repository: PyPSA/atlite

Issue:
Wind Conversion: potential bug when power curve does not end with zero after cutout speed
<!---
SPDX-FileCopyrightText: 2021 The Atlite Authors

SPDX-License-Identifier: CC0-1.0
--->


<!-- Provide a general summary of the issue -->

## Description
<!-- Provide a more detailed introduction to the issue itself, and why you consider it to be a bug -->
<!-- If you can, add a minimal example which reproduces the bug -->
Hi there,

I am not sure whether this is a bug or if I did get something wrong. In the `convert_wind()` function, you use [`np.interp()`](https://numpy.org/doc/stable/reference/generated/numpy.interp.html) in order to fit the wind speeds from `wnd_hub` to the specific power output:

https://github.com/PyPSA/atlite/blob/549d0fdf22925e06ac3f726e801f215d993671dc/atlite/convert.py#L466-L467

I think this assumes that the last value in the power curve (i.e. the value corresponding to the highest wind speed) to be a zero in order to account for cutout at high wind speeds of a wind turbine correctly. By default, np.interp assumes fp[-1] for x > xp[-1] (documented in the "right" argument in np.interp). This will lead to overestimation of specific generation for high wind speed if the power curve does not end with zero. See the following example:

```python
In [1]: import numpy as np
   ...: # example curve from https://github.com/PyPSA/atlite/blob/master/atlite/resources/windturbine/NREL_ReferenceTurbine_2019ORCost_12MW_offshore.yaml
   ...: V = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]
   ...: POW = [0.0, 0.0, 0.0, 0.4, 1.141, 2.189, 3.581, 5.323, 7.579, 10.397, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0]
   ...: P = np.max(POW)
   ...: print("specific generation at 30 m/s:", np.interp(30, V, POW / P))
   ...:
specific generation at 30 m/s: 1.0
```

I noticed that most of the OEDB power curves do not end with a zero in the power curve values (only for 5 turbines in the OEDB table this is true). See below a list of all turbines in OEDB, where the power curve does not end with a zero:
```python
In [2]: import requests
   ...: import pandas as pd
   ...: import ast
   ...: oedb_json = requests.get("https://openenergy-platform.org/api/v0/schema/supply/tables/wind_turbine_library/rows")
   ...: oedb = pd.DataFrame.from_dict(oedb_json.json())
   ...: 
   ...: for row in oedb[oedb["has_power_curve"]].itertuples():
   ...:     if list(ast.literal_eval(row.power_curve_values))[-1] != 0:
   ...:         print(f"id: {row.id}, manufacturer: {row.manufacturer}, name: {row.name}")
   ...: 
id: 0, manufacturer: Enercon, name: E-141/4200 EP 4
id: 1, manufacturer: Enercon, name: E-126/4200 EP4
id: 2, manufacturer: Enercon, name: E-101/3500 E2
id: 3, manufacturer: Enercon, name: E-115/3200
id: 5, manufacturer: Enercon, name: E-115/3000
id: 6, manufacturer: Enercon, name: E-82/3000
id: 7, manufacturer: Enercon, name: E-92/2350
id: 8, manufacturer: Enercon, name: E-82/2350
id: 9, manufacturer: Enercon, name: E-82/2300 E2
id: 10, manufacturer: Enercon, name: E-70/2300 E4
id: 11, manufacturer: Enercon, name: E-82/2000 E2
id: 12, manufacturer: Enercon, name: E-70/2000
id: 13, manufacturer: Enercon, name: E-53/800
id: 15, manufacturer: Nordex, name: None
id: 16, manufacturer: Nordex, name: None
id: 17, manufacturer: Nordex, name: None
id: 20, manufacturer: Nordex, name: None
id: 63, manufacturer: Siemens, name: SWT-3.3-130
id: 22, manufacturer: Nordex, name: None
id: 23, manufacturer: Nordex, name: None
id: 33, manufacturer: Vestas, name: V164-8.0 MW
id: 97, manufacturer: Senvion/REpower, name: MM92
id: 42, manufacturer: Vestas, name: V126-3.45 MW
id: 43, manufacturer: Vestas, name: V126-3.3 MW
id: 44, manufacturer: Vestas, name: V126-3.0 MW
id: 47, manufacturer: Vestas, name: V117-3.45 MW
id: 48, manufacturer: Vestas, name: V117-3.3 MW
id: 49, manufacturer: Vestas, name: V112-3.45 MW
id: 50, manufacturer: Vestas, name: V112-3.3 MW
id: 51, manufacturer: Vestas, name: V112-3.075 MW
id: 52, manufacturer: Vestas, name: V112-3.0 MW
id: 53, manufacturer: Vestas, name: V90-3.0 MW
id: 54, manufacturer: Vestas, name: V90-2.0 MW
id: 55, manufacturer: Vestas, name: V80-2.0 MW
id: 60, manufacturer: Siemens, name: SWT-3.6-130
id: 61, manufacturer: Siemens, name: SWT-3.6-120
id: 64, manufacturer: Siemens, name: SWT-3.2-113
id: 67, manufacturer: Siemens, name: SWT-3.15-142
id: 73, manufacturer: Siemens, name: SWT-2.3-113
id: 84, manufacturer: Senvion/REpower, name: 6.3M152
id: 86, manufacturer: Senvion/REpower, name: 6.2M126
id: 89, manufacturer: Senvion/REpower, name: 3.4M114
id: 90, manufacturer: Senvion/REpower, name: 3.4M104
id: 93, manufacturer: Senvion/REpower, name: 3.2M114 NES
id: 94, manufacturer: Senvion/REpower, name: 3.2M122
id: 95, manufacturer: Senvion/REpower, name: 3.0M122
id: 96, manufacturer: Senvion/REpower, name: MM100
id: 114, manufacturer: GE Wind, name: GE 3.2-130
id: 124, manufacturer: Eno, name: ENO 126 3.5
id: 116, manufacturer: GE Wind, name: GE 2.75-120
id: 117, manufacturer: GE Wind, name: GE 2.75-103
id: 119, manufacturer: GE Wind, name: GE 2.5-120
id: 127, manufacturer: Eno, name: ENO 114 3.5
id: 128, manufacturer: Eno, name: ENO 100
id: 131, manufacturer: Vestas, name: None
id: 132, manufacturer: Enercon, name: E-48 800
id: 133, manufacturer: Vestas, name: V90/2000 GS
id: 134, manufacturer: Vestas, name: V100/1800
id: 135, manufacturer: Vestas, name: V100/1800 GS
id: 136, manufacturer: Vestas, name: V117-3.6
id: 137, manufacturer: aerodyn, name: SCD 168 8000
id: 138, manufacturer: Enercon, name: E-126/7500
id: 140, manufacturer: IEA, name: IEA 15 MW offshore reference turbine
```

The same is true for some of the atlite-builtin power curves:

```python
In [3]: import atlite
   ...: for turbine in atlite.windturbines.keys():
   ...:     cfg = atlite.resource.get_windturbineconfig(turbine)
   ...:     if cfg["POW"][-1] != 0:
   ...:         print(turbine)
   ...: 
NREL_ReferenceTurbine_2016CACost_10MW_offshore
NREL_ReferenceTurbine_2016CACost_6MW_offshore 
NREL_ReferenceTurbine_2016CACost_8MW_offshore 
NREL_ReferenceTurbine_2019ORCost_12MW_offshore
NREL_ReferenceTurbine_2019ORCost_15MW_offshore
NREL_ReferenceTurbine_2020ATB_12MW_offshore   
NREL_ReferenceTurbine_2020ATB_15MW_offshore   
NREL_ReferenceTurbine_2020ATB_18MW_offshore
```
## Expected Behavior
<!-- Tell us what should happen -->
I think atlite should do one of the two following options:
   1) check if a high wind speed cutoff in the power curve is present and warn the user if this is not the case
   2) by default assume a specific generation of 0 for wind speeds higher than the last value in the power curve

## Potential Fixes

In case of option 2), one could either modify the `get_windturbineconfig()` function to append a zero to the end of the power curve or change the `_interpolate()` function to:

```python
def _interpolate(da):
    return np.interp(da, V, POW / P, right=0)
```

Discussion:
euronion: Thanks @joAschauer !

A few suggestions:
* make the allowed values for `add_cutoff = True|False|"warn"` with `warn` being the default (-> backwards compatibility)
* Move the warning from `convert_wind(...)` and implementation from `wind(..)` both into `get_windturbineconfig(turbine, add_cutoff)` right before the `return` statement at the end.
https://github.com/PyPSA/atlite/blob/549d0fdf22925e06ac3f726e801f215d993671dc/atlite/resource.py#L35-L75

(Since that method is called internally by `.wind(...)` it will yield the same result, but people can also pre-load the turbine and modify the power curve manually before passing it to `.wind(...)`

* Only print the warning if `add_cutoff=warn`

* We could change the behaviour in an upcoming version to `add_cutoff = True` by default. What do you think? In that case we'd need a warning that the implementation will change.
joAschauer: Hi @euronion,

thanks for the feedback, I can get back to this soon. 

I put the warning and modification code in `convert_wind()` and `.wind()` respectively, since the user can also pass a dict as turbine argument to `.wind()`. In this case, a missing cutoff will be overlooked if the code is in `get_windturbineconfig()`. If this is acceptable to you, I can change it as you suggested above.

I can also try to add a deprecation warning, I think it's reasonable to default to `add_cutoff=True` in future versions. 
euronion: Good thinking with the custom power curves!

I generally lean towards making everyone responsible for the correctness of their own input (custom power curves), but I think keeping the warning is also fair.

How about putting the warning + the modification into a single function and adding that function to `get_windturbineconfig()` (warning + modifying) and `.wind()` (warning only) then? Then we have both cases covered with the same functionality.
joAschauer: Hi @euronion, I made an update and moved the turbine config validation to a new function. I did not implement the `"warn"` option for `add_cutoff` as I think it is clearer to leave the two options True/False. Keeping the default at False for now will not introduce any backwards incompatible changes. 

I only call the validation function in `get_windturbineconfig()` cause otherwise we would end up with two concurrently thrown equal warnings.
euronion: Very nice!

A few final points:

* The correct term is "cut-out", not "cut-off", could you adjust all occurrences? (Sorry for using the incorrect term until now). It might also be better to use `add_cutout_windspeed` to avoid confusion with `atlite.Cutout` (your call!)
* FutureWarning: Announce change for `v0.2.13` instead of "future version of atlite"
* Do you think we should add a check to `_validate_turbine_config_dict(...)` to ensure `V` is sorted ascending order?
* Can you add brief release notes into RELEASE_NOTES.rst?
* A test for the CI would be nice, this should suffice (add to `test/test_resource.py`:
```
def test_windturbineconfig_add_cutout():
    assert atlite.resource.get_windturbineconfig({'V': [ 0,  25], 'POW': [0.,1.], 'hub_height': 1., 'P': 1.}, add_cutout=True)["POW"][-1] == 0
```

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index fa07480..b333723 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -45,7 +45,7 @@ repos:
 
 # Formatting with "black" coding style
 - repo: https://github.com/psf/black
-  rev: 23.7.0
+  rev: 23.9.1
   hooks:
   # Format Python files
   - id: black
diff --git a/RELEASE_NOTES.rst b/RELEASE_NOTES.rst
index ffef895..839c144 100644
--- a/RELEASE_NOTES.rst
+++ b/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
 ==============
 
diff --git a/atlite/convert.py b/atlite/convert.py
index 26d443a..1baaa07 100644
--- a/atlite/convert.py
+++ b/atlite/convert.py
@@ -480,7 +480,7 @@ def convert_wind(ds, turbine):
     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)
diff --git a/atlite/resource.py b/atlite/resource.py
index ce1ecdd..fb165d9 100644
--- a/atlite/resource.py
+++ b/atlite/resource.py
@@ -11,6 +11,7 @@ panel configurations.
 import json
 import logging
 import re
+import warnings
 from operator import itemgetter
 from pathlib import Path
 
@@ -32,13 +33,13 @@ SOLARPANEL_DIRECTORY = RESOURCE_DIRECTORY / "solarpanel"
 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 windturbine_smooth(turbine, params=None):
     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 patch:
diff --git a/test/test_resource.py b/test/test_resource.py
index a2977c4..0b4be09 100644
--- a/test/test_resource.py
+++ b/test/test_resource.py
@@ -11,7 +11,7 @@ Created on Tue Jun 22 10:46:27 2021.
 """
 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