Add wind power potential and wind profile indicators

CHANGES.rst

@@ -8,6 +8,7 @@ Contributors to this version: Éric Dupuis (:user:`coxipi`), Trevor James Smith
 
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* Add ``wind_power_potential`` to estimate the potential for wind power production given wind speed at the turbine hub height and turbine specifications, along with  ``wind_profile`` to estimate the wind speed at different heights based on wind speed at a reference height. (:issue:`1458`, :pull:`1471`)
 * `xclim` now has a dedicated console command for prefetching testing data from `xclim-testdata` with branch options (e.g.: `$ xclim prefetch_testing_data --branch some_development_branch`). This command can be used to download the testing data to a local cache, which can then be used to run the testing suite without internet access or in "offline" mode. For more information, see the contributing documentation section for `Updating Testing Data`. (:issue:`1468`, :pull:`1473`).
 * The testing suite now offers a means of running tests in "offline" mode (using `pytest-socket <https://github.com/miketheman/pytest-socket>`_ to block external connections). This requires a local copy of `xclim-testdata` to be present in the user's home cache directory and for certain `pytest` options and markers to be set when invoked. For more information, see the contributing documentation section for `Running Tests in Offline Mode`. (:issue:`1468`, :pull:`1473`).
 * The `SKIP_NOTEBOOKS` flag to speed up docs builds is now documented. See the contributing documentation section `Get Started!` for details. (:issue:`1470`, :pull:`1476`).
@@ -35,6 +36,7 @@ Internal changes
 * GitHub deployment workflows now employs use of deployment environments for workflow security and uses the `Trusted Publisher <https://docs.pypi.org/trusted-publishers/using-a-publisher/>`_ feature to sign and publish the `xclim` wheel and source distributions. (:pull:`1469`).
 * Mastodon publishing now uses `chuhlomin/render-template <https://github.com/chuhlomin/render-template>`_ and a standard formatting markdown document to format Mastodon toots. (:pull:`1469`).
 
+
 v0.45.0 (2023-09-05)
 --------------------
 Contributors to this version: David Huard (:user:`huard`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`), Gabriel Rondeau-Genesse (:user:`RondeauG`), Marco Braun (:user:`vindelico`), Éric Dupuis (:user:`coxipi`).

docs/references.bib

@@ -2046,3 +2046,30 @@ @misc{usda_2012
 	author = {{USDA Agricultural Research Service}},
 	year = {2012},
 }
+
+@article{chen_2020,
+title = {Impacts of climate change on wind resources over North America based on NA-CORDEX},
+journal = {Renewable Energy},
+volume = {153},
+pages = {1428-1438},
+year = {2020},
+issn = {0960-1481},
+doi = {https://doi.org/10.1016/j.renene.2020.02.090},
+url = {https://www.sciencedirect.com/science/article/pii/S0960148120302858},
+author = {Liang Chen},
+keywords = {Climate change, Wind energy, NA-CORDEX, North America}
+}
+
+@article{tobin_2018,
+doi = {10.1088/1748-9326/aab211},
+url = {https://dx.doi.org/10.1088/1748-9326/aab211},
+year = {2018},
+month = {apr},
+publisher = {IOP Publishing},
+volume = {13},
+number = {4},
+pages = {044024},
+author = {I Tobin and W Greuell and S Jerez and F Ludwig and R Vautard and M T H van Vliet and F-M Bréon},
+title = {Vulnerabilities and resilience of European power generation to 1.5°C, 2°C and 3°C warming},
+journal = {Environmental Research Letters}
+}

tests/test_atmos.py

@@ -253,6 +253,41 @@ def test_wind_chill_index(atmosds):
     )
 
 
+def test_wind_profile(atmosds):
+    out = atmos.wind_profile(
+        wind_speed=atmosds.sfcWind, h_r="10 m", h="100 m", alpha=1 / 7
+    )
+    assert out.attrs["units"] == "m s-1"
+    assert (out > atmosds.sfcWind).all()
+
+
+def test_wind_power_potential(atmosds):
+    out = atmos.wind_power_potential(wind_speed=atmosds.sfcWind)
+    assert out.attrs["units"] == ""
+    assert (out >= 0).all()
+    assert (out <= 1).all()
+
+
+def test_wind_power_potential_from_3h_series():
+    """Test a typical computation workflow from 3-hourly time series to daily power production in MWh."""
+    from xclim.core.units import convert_units_to
+    from xclim.indices.generic import select_resample_op
+    from xclim.testing.helpers import test_timeseries
+
+    w = test_timeseries(
+        np.ones(96) * 15, variable="sfcWind", start="7/1/2000", units="m s-1", freq="3H"
+    )
+    out = atmos.wind_power_potential(wind_speed=w)
+
+    # Multiply with nominal capacity
+    power = out * 100
+    power.attrs["units"] = "MW"
+    annual_power = convert_units_to(
+        select_resample_op(power, op="sum", freq="D"), "MWh"
+    )
+    assert (annual_power == 100 * 24).all()
+
+
 class TestDrynessIndex:
     def test_simple(self, atmosds):
         ds = atmosds.isel(location=3)

tests/test_indices.py

@@ -3626,3 +3626,29 @@ def test_late_frost_days_lat(self, tasmin_series):
         tasmin = tasmin_series(np.array([-1, 1, 2, -4, 0]) + K2C, start="30/3/2023")
         lfd = xci.frost_days(tasmin, date_bounds=("04-01", "06-30"))
         np.testing.assert_allclose(lfd, 1)
+
+
+class TestWindProfile:
+    def test_simple(self, sfcWind_series):
+        a = np.linspace(0, 100)
+        v = xci.wind_profile(sfcWind_series(a), h="100 m", h_r="10 m")
+        np.testing.assert_allclose(v, a * 10 ** (1 / 7))
+
+
+class TestWindPowerPotential:
+    def test_simple(self, sfcWind_series):
+        v = [2, 6, 20, 30]
+        p = xci.wind_power_potential(
+            sfcWind_series(v, units="m/s"), cut_in="4 m/s", rated="8 m/s"
+        )
+        np.testing.assert_allclose(p, [0, (6**3 - 4**3) / (8**3 - 4**3), 1, 0])
+
+        # Test discontinuities at the default thresholds
+        v = np.array([3.5, 15])
+        a = sfcWind_series(v - 1e-7, units="m/s")
+        b = sfcWind_series(v + 1e-7, units="m/s")
+
+        pa = xci.wind_power_potential(a)
+        pb = xci.wind_power_potential(b)
+
+        np.testing.assert_array_almost_equal(pa, pb, decimal=6)

xclim/core/units.py

@@ -579,15 +579,15 @@ def to_agg_units(
             units="", is_dayofyear=np.int32(1), calendar=get_calendar(orig)
         )
 
-    elif op in ["count", "integral"]:
+    elif op in ["count", "integral", "sum"]:
         m, freq_u_raw = infer_sampling_units(orig[dim])
         orig_u = str2pint(orig.units)
         freq_u = str2pint(freq_u_raw)
         out = out * m
 
         if op == "count":
             out.attrs["units"] = freq_u_raw
-        elif op == "integral":
+        elif op in ["integral", "sum"]:
             out.attrs["units"] = pint2cfunits(orig_u * freq_u)
     else:
         raise ValueError(

xclim/data/fr.json

@@ -792,6 +792,18 @@
     "long_name": "Indice forêt météo",
     "description": "Évaluation numérique de l'intensité du feu."
   },
+  "WIND_POWER_POTENTIAL": {
+    "long_name": "Potentiel de production éolienne",
+    "description": "Potentiel de production éolienne estimé par une courbe de puissance semi-idéalisée d'une turbine, paramétrée par une vitesse minimale de démarrage {cut_in}, la vitesse nominale {rated}, et la vitesse d'arrêt {cut_out}.",
+    "title": "Potentiel de production éolienne",
+    "abstract": "Estimation à partir d'une loi de puissance semi-idéalisée de la production d'énergie éolienne selon la vitesse du vent."
+  },
+  "WIND_PROFILE": {
+    "long_name": "Vitesse du vent à la hauteur {h}",
+    "description": "Vitesse du vent à une hauteur de {h} calculée de la vitesse à {h_r} par un profil de loi de puissance.",
+    "title": "Profil du vent",
+    "abstract": "Estimation de la vitesse du vent à partir du vent à une hauteur de référence."
+  },
   "WIND_SPEED_FROM_VECTOR": {
     "title": "Vitesse et direction du vent à partir du vecteur vent",
     "abstract": "Calcul de la magnitude et la direction de la vitesse du vent à partir des deux composantes ouest-est et sud-nord."

xclim/data/variables.yml

@@ -1,4 +1,10 @@
 variables:
+  air_density:
+    canonical_units: kg m-3
+    cell_methods: "time: mean"
+    description: Air density.
+    dimensions: "[density]"
+    standard_name: air_density
   areacello:
     canonical_units: m2
     cell_methods: "area: sum"
@@ -316,6 +322,12 @@ variables:
     description: Northward component of the wind velocity (at the surface).
     dimensions: "[speed]"
     standard_name: northward_wind
+  wind_speed:
+    canonical_units: m s-1
+    cell_methods: "time: mean"
+    description: Wind speed.
+    dimensions: "[speed]"
+    standard_name: wind_speed
   wsgsmax:
     cmip6: False
     canonical_units: m s-1

xclim/indicators/atmos/_conversion.py

@@ -28,6 +28,8 @@
     "water_budget",
     "water_budget_from_tas",
     "wind_chill_index",
+    "wind_power_potential",
+    "wind_profile",
     "wind_speed_from_vector",
     "wind_vector_from_speed",
 ]
@@ -123,6 +125,30 @@ def cfcheck(self, **das):
     compute=indices.sfcwind_2_uas_vas,
 )
 
+wind_power_potential = Converter(
+    title="Wind power potential",
+    identifier="wind_power_potential",
+    units="",
+    long_name="Wind power potential",
+    description="Wind power potential using a semi-idealized turbine power curve using a cut_in speed of {cut_in}, "
+    "a rated speed of {rated}, and a cut_out speed of {cut_out}.",
+    cell_methods="",
+    abstract="Calculation of the wind power potential using a semi-idealized turbine power curve.",
+    compute=indices.wind_power_potential,
+)
+
+wind_profile = Converter(
+    title="Wind profile",
+    identifier="wind_profile",
+    var_name="wind_speed",
+    units="m s-1",
+    standard_name="wind_speed",
+    long_name="Wind speed at height {h}",
+    description="Wind speed at a height of {h} computed from the wind speed at {h_r} using a power law profile.",
+    cell_methods="",
+    abstract="Calculation of the wind speed at a given height from the wind speed at a reference height.",
+    compute=indices.wind_profile,
+)
 
 saturation_vapor_pressure = Converter(
     title="Saturation vapour pressure (e_sat)",

xclim/indices/_conversion.py

@@ -53,6 +53,8 @@
     "uas_vas_2_sfcwind",
     "universal_thermal_climate_index",
     "wind_chill_index",
+    "wind_power_potential",
+    "wind_profile",
 ]
 
 
@@ -2017,3 +2019,161 @@ def mean_radiant_temperature(
         0.25,
     )
     return mrt.assign_attrs({"units": "K"})
+
+
+@declare_units(wind_speed="[speed]", h="[length]", h_r="[length]")
+def wind_profile(
+    wind_speed: xr.DataArray,
+    h: Quantified,
+    h_r: Quantified,
+    method: str = "power_law",
+    **kwds,
+):
+    r"""Wind speed at a given height estimated from the wind speed at a reference height.
+
+    Estimate the wind speed based on a power law profile relating wind speed to height above the surface.
+
+    Parameters
+    ----------
+    wind_speed : xarray.DataArray
+        Wind speed at the reference height.
+    h : Quantified
+        Height at which to compute the wind speed.
+    h_r : Quantified
+        Reference height.
+    method : {"power_law"}
+        Method to use. Currently only "power_law" is implemented.
+    kwds : dict
+        Additional keyword arguments to pass to the method. For power_law, this is alpha, which takes a default value
+        of 1/7, but is highly variable based on topography, surface cover and atmospheric stability.
+
+    Notes
+    -----
+    The power law profile is given by
+
+    .. math::
+
+        v = v_r \left( \frac{h}{h_r} \right)^{\alpha},
+
+    where :math:`v_r` is the wind speed at the reference height, :math:`h` is the height at which the wind speed is
+    desired, and :math:`h_r` is the reference height.
+
+    """
+    # Convert units to meters
+    h = convert_units_to(h, "m")
+    h_r = convert_units_to(h_r, "m")
+
+    if method == "power_law":
+        alpha = kwds.pop("alpha", 1 / 7)
+        out = wind_speed * (h / h_r) ** alpha
+        out.attrs["units"] = wind_speed.attrs["units"]
+        return out
+    else:
+        raise NotImplementedError(f"Method {method} not implemented.")
+
+
+@declare_units(
+    wind_speed="[speed]",
+    air_density="[air_density]",
+    cut_in="[speed]",
+    rated="[speed]",
+    cut_out="[speed]",
+)
+def wind_power_potential(
+    wind_speed: xr.DataArray,
+    air_density: xr.DataArray = None,
+    cut_in: Quantified = "3.5 m/s",
+    rated: Quantified = "13 m/s",
+    cut_out: Quantified = "25 m/s",
+) -> xr.DataArray:
+    r"""Wind power potential estimated from an idealized wind power production factor.
+
+    The actual power production of a wind farm can be estimated by multiplying its nominal (nameplate) capacity by the
+    wind power potential, which depends on wind speed at the hub height, the turbine specifications and air density.
+
+    Parameters
+    ----------
+    wind_speed : xarray.DataArray
+        Wind speed at the hub height. Use the `wind_profile` function to estimate from the surface wind speed.
+    air_density: xarray.DataArray
+        Air density at the hub height. Defaults to 1.225 kg/m³. This is worth changing if applying in cold or
+        mountainous regions with non-standard air density.
+    cut_in : Quantified
+        Cut-in wind speed. Default is 3.5 m/s.
+    rated : Quantified
+        Rated wind speed. Default is 13 m/s.
+    cut_out : Quantified
+        Cut-out wind speed. Default is 25 m/s.
+
+    Returns
+    -------
+    xr.DataArray
+      The power production factor. Multiply by the nominal capacity to get the actual power production.
+
+    See Also
+    --------
+    wind_profile : Estimate wind speed at the hub height from the surface wind speed.
+
+    Notes
+    -----
+    This estimate of wind power production is based on an idealized power curve with four wind regimes specified
+    by the cut-in wind speed (:math:`u_i`), the rated speed (:math:`u_r`) and the cut-out speed (:math:`u_o`).
+    Power production is zero for wind speeds below the cut-in speed, increases cubically between the cut-in
+    and rated speed, is constant between the rated and cut-out speed, and is zero for wind speeds above the cut-out
+    speed to avoid damage to the turbine :cite:p:`tobin_2018`:
+
+    .. math::
+
+        \begin{cases}
+        0,  &  v < u_i \\
+        (v^3 - u_i^3) / (u_r^3 - u_i^3),  & u_i ≤ v < u_r \\
+        1, & u_r ≤ v < u_o \\
+        0, & v ≥ u_o
+        \end{cases}
+
+    For non-standard air density (:math:`\rho`), the wind speed is scaled using
+    :math:`v_n = v \left( \frac{\rho}{\rho_0} \right)^{1/3}`.
+
+    The temporal resolution of wind time series has a significant influence on the results: mean daily wind
+    speeds yield lower values than hourly wind speeds. Note however that percent changes in the wind power potential
+    climate projections are similar across resolutions :cite:p:`chen_2020`.
+
+    To compute the power production, multiply the power production factor by the nominal
+    turbine capacity (e.g. 100), set the units attribute (e.g. "MW"), resample and sum with
+    `xclim.indices.generic.select_resample_op(power, op="sum", freq="D")`, then convert to
+    the desired units (e.g. "MWh") using `xclim.core.units.convert_units_to`.
+
+    References
+    ----------
+    :cite:cts:`chen_2020,tobin_2018`.
+
+    """
+    # Convert units
+    cut_in = convert_units_to(cut_in, wind_speed)
+    rated = convert_units_to(rated, wind_speed)
+    cut_out = convert_units_to(cut_out, wind_speed)
+
+    # Correct wind speed for air density
+    if air_density is not None:
+        default_air_density = convert_units_to("1.225 kg/m^3", air_density)
+        f = (air_density / default_air_density) ** (1 / 3)
+    else:
+        f = 1
+
+    v = wind_speed * f
+
+    out = xr.apply_ufunc(_wind_power_factor, v, cut_in, rated, cut_out)
+    out.attrs["units"] = ""
+    return out
+
+
+@vectorize
+def _wind_power_factor(v, cut_in, rated, cut_out):
+    """Wind power factor function"""
+    if v < cut_in:
+        return 0.0
+    if v < rated:
+        return (v**3 - cut_in**3) / (rated**3 - cut_in**3)
+    if v < cut_out:
+        return 1.0
+    return 0.0