"""Python interface to WAsP's AEP calculation
The Annual Energy Production (AEP) of a wind turbine is usually
expressed in terms of MWh or GWh and depends on the wind climate
(described by sectorwise Weibull A's and k's and frequencies),
the air density and the wind turbine power curve. Usually one might
have only a few power curves from the manufacturer, which are valid
for specific air densities. Since WAsP 12 a inter- or extrapolation
routine is used to estimate the power curves for different air
densities based on one or several tables with power curves.
"""
from functools import partial, wraps
from functools import wraps
import warnings
import numpy as np
import geopandas as gpd
import xarray as xr
try:
import py_wake
except ModuleNotFoundError:
py_wake = None # type: ignore
try:
from numba import njit
except ModuleNotFoundError:
# Set njit to identity decorator (do nothing to function)
def njit(func):
@wraps(func)
def wrapper_njit(*args, **kwargs):
return func(*args, **kwargs)
return wrapper_njit
import windkit as wk
from windkit.metadata import (
update_var_attrs,
update_history,
_AEP_ATTRS,
_WF_FLOW_MAP_ATTRS,
)
from ..core import Rvea0326
_FORT_AEP = Rvea0326.calc_points.aep_points
_FORT_POWER_INTEGRAL = Rvea0326.aep_externals_mod.aepwp
AIR_DENSITY_DEFAULT = 1.225
def _prepare_wtg_for_fort(wtg):
"""Convert a PyWAsP Wind Turbine Generator dataset
to numpy arrays for ingestion to waspcore fortran AEP calculation
Parameters
----------
wtg : xr.Dataset
PyWAsP Wind Turbine Generator Dataset.
Returns
-------
int :
Number of turbine modes
int :
Number of wind speed data points
np.ndarray :
1D array of floats containing Wind speed values
np.ndarray :
1D array of floats containing Power output values
"""
if "mode" in wtg.coords:
n_mode = wtg.coords["mode"].size
else:
n_mode = 1
n_wspd = wtg.coords["wind_speed"].size
wspds = np.repeat(wtg["wind_speed"].values[None, :], n_mode, axis=0)
power = wtg["power_output"].values
if power.ndim == 1:
power = np.repeat(power[None, :], n_mode, axis=0)
table_points = np.repeat(n_wspd, n_mode)
table_dens = np.array(wtg.air_density)
wspd_arr = np.column_stack((table_points, wspds)).ravel()
power_arr = np.column_stack((table_dens, power)).ravel()
control_system = int(wtg.regulation_type.data)
return (n_mode, n_wspd, wspd_arr, power_arr, control_system)
def _gross_aep(wwc, wtg, use_sectors=True, interpolate=False, air_density=None):
"""Calculate annual energy production (AEP), using the WAsP core
Fortran implementation, from Weibull wind climate(s)
and Wind Turbine Generator(s) or WindTurbines object.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
wwc : xr.Dataset
Weibull Wind Climate (WWC) Dataset containing
sectorwise Weibull A's and k's and frequencies.
wtg : xr.Dataset
One Wind Turbine Generator as a xr.Dataset
use_sectors : bool
Whether to use sectors ("A" and "k") or total ("A_combined" and "k_combined")
interpolate : bool
Whether to use power curve interpolation.
air_density : float
Air density to use for the AEP calculation. If not passed, the air density
from the WWC will be used. If air_density is passed, it will be used,
and overwrite the air density in the WWC, if present.
Returns
-------
xr.Dataset
Gross Annual Energy Production in GWh for each turbine/location.
Raises
------
TypeError
If wtg is not a xr.Dataset or wk.WindTurbines object.
ValueError
If required variables are not in the WWC.
ValueError
If air_density is not in the WWC and no explicit value was passed.
ValueError
If any license error occurs during the AEP calculation.
"""
wwc = wwc.copy()
# We always convert to wk.WindTurbines
# If wtg is a xr.Dataset, we assume the WWC locations are the same as the WTG locations
# and we create a wk.WindTurbines object with the WWC locations and the WTG power curves
if not isinstance(wtg, xr.Dataset):
raise TypeError(f"wtg must be a xr.Dataset, not {type(wtg)}")
if air_density is not None:
if "air_density" not in wwc.data_vars:
wwc["air_density"] = (("point",), np.ones([wwc.point.size]) * air_density)
else:
warnings.warn(
"air_density is already in the WWC but also passed explicit. It will be overwritten."
)
elif "air_density" not in wwc.data_vars:
raise ValueError(
"air_density is not in the WWC and no explicit value was passed."
)
if use_sectors:
required_vars = ["A", "k", "wdfreq", "air_density"]
tvars = ("sector", "point")
else:
required_vars = ["A_combined", "k_combined", "air_density"]
tvars = ("point",)
if not all(k in wwc.data_vars for k in required_vars):
raise ValueError(
f"""Need to have {required_vars}
variables to calculate the AEP."""
)
wwc = wwc.transpose(..., *tvars)
if not use_sectors:
output_locs = (
wwc[required_vars].expand_dims("sector").transpose("sector", "point")
)
output_locs["wdfreq"] = (
("sector", "point"),
np.ones_like(output_locs.A_combined),
)
output_locs = output_locs.rename({"A_combined": "A", "k_combined": "k"})
output_locs["air_density"] = output_locs["air_density"].isel(sector=0)
else:
output_locs = wwc[required_vars]
n_tables, n_points, wspd_arr, power_arr, control_system = _prepare_wtg_for_fort(wtg)
aep, ierror = xr.apply_ufunc(
_FORT_AEP,
output_locs.A,
output_locs.k,
output_locs.wdfreq,
output_locs.air_density,
n_tables,
n_points,
wspd_arr,
power_arr,
control_system,
1.0,
int(interpolate),
input_core_dims=[
["sector", "point"],
["sector", "point"],
["sector", "point"],
["point"],
[],
[],
[],
[],
[],
[],
[],
],
output_core_dims=[["point"], []],
)
if ierror != 0:
raise ValueError(f"License Error {ierror}")
# Convert ot xr.DataArray
output_locs["gross_AEP"] = (
("point",), # TODO: switch to sectorwise results
aep.data * 1e-9,
)
output_locs = update_var_attrs(output_locs, _AEP_ATTRS)
output_locs = update_history(output_locs)
return output_locs[["gross_AEP"]]
def _gross_aep_wind_turbines(
wwc, wtg_dict, wind_turbines, use_sectors=True, interpolate=False, air_density=None
):
"""
Calculate the gross AEP for a set of wind turbines. This functions is more flexible than
the :py:func:`gross_aep` function, as it allows for multiple wind turbines generators and wind turbine groups
to be passed in. This is useful when you have multiple wind turbine types or groups in your wind farm.
Parameters
----------
wwc : xr.Dataset
Weibull Wind Climate (WWC) Dataset containing
sectorwise Weibull A's and k's and frequencies. Must have "point" dimension and be
ordered in the same way as ``wind_turbines``.
wtg_dict : dict
Dictionary mapping wtg_keys to Wind Turbine Generator datasets.
wind_turbines : xr.Dataset
Wind Turbine locations, hub heights, group_id's, and wtg_keys matching the keys in ``wtg_dict``.
use_sectors : bool
Whether to use sectors ("A" and "k") or total ("A_combined" and "k_combined")
interpolate : bool
Whether to use power curve interpolation.
air_density : float
Air density to use for the AEP calculation. If not passed, the air density
from the WWC will be used. If air_density is passed, it will be used,
and overwrite the air density in the WWC, if present.
Returns
-------
xr.Dataset
Gross Annual Energy Production in GWh for each turbine/location.
Raises
------
TypeError
If wwc is not a xr.Dataset, wind_turbines is not a xr.Dataset, or wtg_dict is not a dict.
ValueError
If any of the wind turbine keys in the wind_turbines object are not in the wtg_dict.
"""
if not isinstance(wwc, xr.Dataset):
raise TypeError(f"wwc must be a xr.Dataset, not {type(wwc)}")
if not isinstance(wind_turbines, xr.Dataset):
raise TypeError(
f"wind_turbines must be a xr.Dataset, not {type(wind_turbines)}"
)
if not isinstance(wtg_dict, dict):
raise TypeError(f"wtg_dict must be a dict, not {type(wtg_dict)}")
wwc = wwc.copy()
wind_turbines = wind_turbines.copy()
# Ensure all wtg_key values in the dataframe are in the supplied dictionary.
if not all(k in wtg_dict.keys() for k in wind_turbines.wtg_key.data):
raise ValueError(
"Wind_turbines object includes wtg_keys that aren't in the wtg_dict."
)
group_ids = np.unique(wind_turbines.group_id.data)
order = np.arange(wind_turbines.point.size)
reorder = []
results = []
for group_id in group_ids:
indx = np.where(wind_turbines.group_id.data == group_id)[0]
reorder.append(order[indx])
wwc_grp = wwc.isel(point=indx)
wtg_grp = wtg_dict[str(wind_turbines.isel(point=indx[0]).wtg_key.data)]
aep_grp = _gross_aep(
wwc_grp,
wtg_grp,
use_sectors=use_sectors,
interpolate=interpolate,
air_density=air_density,
)
results.append(aep_grp)
aep = xr.concat(results, dim="point")
aep = aep.isel(point=np.argsort(np.concatenate(reorder)))
return aep
[docs]
@wk.spatial.stack_then_unstack
def gross_aep(
wwc,
wtg,
/,
wind_turbines=None,
*,
use_sectors=True,
interpolate=False,
air_density=None,
):
"""Calculate annual energy production (AEP), using the WAsP core
Fortran implementation, from Weibull wind climate(s)
and Wind Turbine Generator(s) or WindTurbines object.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
wwc : xr.Dataset
Weibull Wind Climate (WWC) Dataset containing
sectorwise Weibull A's and k's and frequencies.
wtg : xr.Dataset, dict, wk.WindTurbines
- If a xr.Dataset is passed it should be a WTG formatted dataset.
PyWAsP will assume WWC points match the turbines hub height and that the
WTG is the same for all turbines
- If a dict of wtg_keys and WTG xr.Datasets is passed, a ``turbines``
dataset must be passed as well, and the key's must match between the two.
- A wk.WindTurbines object can also be passed, but is deprecated and
will be removed in a future version.
wind_turbines : xr.Dataset, optional
Wind Turbine locations, hub heights, group_id's, and wtg_keys.
use_sectors : bool
Whether to use sectors ("A" and "k") or total ("A_combined" and "k_combined")
interpolate : bool
Whether to use power curve interpolation.
air_density : float
Air density to use for the AEP calculation. If not passed, the air density
from the WWC will be used. If air_density is passed, it will be used,
and overwrite the air density in the WWC, if present.
Returns
-------
xr.Dataset
Gross Annual Energy Production in GWh for each turbine/location.
Raises
------
TypeError
If wtg is not a xr.Dataset or wk.WindTurbines object.
ValueError
If required variables are not in the WWC.
ValueError
If air_density is not in the WWC and no explicit value was passed.
ValueError
If any license error occurs during the AEP calculation.
"""
if wind_turbines is not None:
if not isinstance(wtg, dict):
raise TypeError(
f"wtg must be a dict of wtg_keys and xr.Dataset values when passing 'wind_turbines', not {type(wtg)}"
)
return _gross_aep_wind_turbines(
wwc,
wtg,
wind_turbines,
use_sectors=use_sectors,
interpolate=interpolate,
air_density=air_density,
)
else:
# TODO: When wk.WindTurbines are removed, this block can also be removed.
if isinstance(wtg, wk.WindTurbines):
warnings.warn(
"The WindTurbins class is deprecated and will be removed in a future version. "
+ "Use a wtg dict and xr.Dataset of wind turbines instead. ",
DeprecationWarning,
stacklevel=2,
)
wtg_dict, wind_turbines = wtg.to_wtg_dict_and_turbines_ds()
return _gross_aep_wind_turbines(
wwc,
wtg_dict,
wind_turbines,
use_sectors=use_sectors,
interpolate=interpolate,
air_density=air_density,
)
elif not isinstance(wtg, xr.Dataset):
raise TypeError(
f"When 'wind_turbines' are not passed, 'wtg' must be a xr.Dataset, not {type(wtg)}"
)
return _gross_aep(
wwc,
wtg,
use_sectors=use_sectors,
interpolate=interpolate,
air_density=air_density,
)
def _wind_farm_model_from_string(string):
"""Get py_wake wind farm model from string/name.
Parameters
----------
string : str
Name of predefined wind farm model. Choices includes:
- "PARK1"
- "PARK2_onshore"
- "PARK2_offshore"
Returns
-------
py_wake.wind_farm_model
Wind farm model
Raises
------
NameError
If name is not in the valid list of names.
"""
valid_strings = [
"PARK1",
"PARK2_onshore",
"PARK2_offshore",
]
if string not in valid_strings:
raise NameError(f"No wind farm model found for: '{string}'")
if string == "PARK1":
class _PARK1Deficit(py_wake.deficit_models.noj.NOJDeficit):
def calc_deficit(self, WS_ilk, WS_eff_ilk, ct_ilk, **kwargs):
if not self.deficit_initalized:
self._calc_layout_terms(
WS_ilk=WS_ilk, WS_eff_ilk=WS_eff_ilk, ct_ilk=ct_ilk, **kwargs
)
ct_ilk = np.minimum(ct_ilk, 1) # treat ct_ilk for np.sqrt()
# NOJ: term_numerator_ilk = (1 - np.sqrt(1 - ct_ilk))
with np.errstate(divide="ignore", invalid="ignore"):
term_numerator_ilk = 1 - WS_eff_ilk / WS_ilk * np.sqrt(1 - ct_ilk)
return term_numerator_ilk[:, None] * self.layout_factor_ijlk
wind_farm_model = partial(
py_wake.wind_farm_models.engineering_models.PropagateDownwind,
wake_deficitModel=_PARK1Deficit(
k=0.1,
ct2a=py_wake.deficit_models.utils.ct2a_mom1d,
groundModel=py_wake.ground_models.Mirror(),
rotorAvgModel=py_wake.rotor_avg_models.AreaOverlapAvgModel(),
),
superpositionModel=py_wake.superposition_models.SquaredSum(),
deflectionModel=None,
turbulenceModel=None,
)
if string == "PARK2_onshore":
wind_farm_model = partial(
py_wake.wind_farm_models.engineering_models.PropagateDownwind,
wake_deficitModel=py_wake.deficit_models.noj.NOJLocalDeficit(
a=[0, 0.09],
ct2a=py_wake.deficit_models.utils.ct2a_mom1d,
use_effective_ws=True,
use_effective_ti=False,
rotorAvgModel=py_wake.rotor_avg_models.AreaOverlapAvgModel(),
groundModel=py_wake.ground_models.NoGround(),
),
superpositionModel=py_wake.superposition_models.LinearSum(),
deflectionModel=None,
turbulenceModel=None,
)
if string == "PARK2_offshore":
wind_farm_model = partial(
py_wake.wind_farm_models.engineering_models.PropagateDownwind,
wake_deficitModel=py_wake.deficit_models.noj.NOJLocalDeficit(
a=[0, 0.06],
ct2a=py_wake.deficit_models.utils.ct2a_mom1d,
use_effective_ws=True,
use_effective_ti=False,
rotorAvgModel=py_wake.rotor_avg_models.AreaOverlapAvgModel(),
groundModel=py_wake.ground_models.NoGround(),
),
superpositionModel=py_wake.superposition_models.LinearSum(),
deflectionModel=None,
turbulenceModel=None,
)
return wind_farm_model
@njit
def argrelmax(data):
"""Find the relative maxima of a 1D array.
Parameters
----------
data : np.ndarray
1D array of floats
Returns
-------
np.ndarray
Indices of the relative maxima
"""
result = np.zeros_like(data, dtype=np.bool_)
result[1:-1] = (data[1:-1] > data[:-2]) & (data[1:-1] > data[2:])
return np.nonzero(result)
@njit
def argrelmin(data):
"""Find the relative minima of a 1D array.
Parameters
----------
data : np.ndarray
1D array of floats
Returns
-------
np.ndarray
Indices of the relative minima
"""
result = np.zeros_like(data, dtype=np.bool_)
result[1:-1] = (data[1:-1] < data[:-2]) & (data[1:-1] < data[2:])
return np.nonzero(result)
@njit
def _merge_power_curve_piece(pc_table, ws, ws_eff, power):
"""Function to find P(U_free) from P(U_eff) and U_free(U_eff)
Author: Original Author Tobias Ahsbahs (ttah)
Parameters
----------
pc_table : np.ndarray
2D array of floats of size (N, 2) power curve table
N is the number of entries in the power curve and
wind speed and power are the first and second "columns"
in the array.
ws : np.ndarray
1D array of floats. Free wind speeds.
ws_eff : np.ndarray
1D array of floats. Effective wind speeds.
power : np.ndarray
1D array of floats. Effective power output.
Returns
-------
np.ndarray
2D array of floats of size (M, 2) containing
[M, [U_free, P]] and M is the number of entries.
"""
# Interpolating free stream wind speed for PC entries
# Need ws_eff to be increasing
if ws_eff[-1] - ws_eff[0] < 0:
ws_eff = ws_eff[::-1]
ws = ws[::-1]
power = power[::-1]
pc_mask = (pc_table[:, 0] >= ws_eff.min()) & (pc_table[:, 0] <= ws_eff.max())
ws_pc, power_pc = pc_table[pc_mask, 0], pc_table[pc_mask, 1]
ws_free_pc = np.interp(ws_pc, ws_eff, ws)
# # Merging the entries
pc_table_merged = np.stack(
(
np.concatenate((ws, ws_free_pc)),
np.concatenate((power, power_pc)),
),
axis=1,
)
# # Sorting by ws
pc_table_merged = pc_table_merged[pc_table_merged[:, 0].argsort()]
# # Removing the first element to avoid double entries in the merging
return pc_table_merged[1:, :]
@njit
def _piece_wise_limits(ind_max, power):
"""Find piece wise limits of a power curve with local maxima and minima
Parameters
----------
ind_max : np.ndarray
1D array of floats. Indices of the local maxima.
power : np.ndarray
1D array of floats. Power output.
Returns
-------
np.ndarray:
2D array of integers containing pairs of indicies (the piece wise limits).
"""
# Find where power is zero and a previous value is higher. These are the points
# where the power is zero within the production range
bool_zero = power == 0.0
bool_step_down = ((power[1:] - power[:-1]) < 0) & bool_zero[1:]
# Finding the positions. Add 1 due to shifting above
pos_step_down = np.argwhere(bool_step_down).flatten() + 1
# Position of local minima
local_minima = argrelmin(power)[0].flatten()
# Both zero positions and local minima for piecewise definition
ind_min = np.unique(np.sort(np.concatenate((pos_step_down, local_minima))))
# Removing the last value after cut_out
ind_min = ind_min[: len(ind_max)]
# Adding start 0 and end length of array. Assuming that start and end
# double points are taken care of in the inputs.
ind_min = np.concatenate((np.array([0]), ind_min))
ind_max = np.concatenate((ind_max, np.array([len(power)])))
# Adding together entries of local maxima and minima for piecewise function
piece_wise_limits = np.array(
list(zip(ind_min, ind_max)) + list(zip(ind_max[:-1], ind_min[1:]))
)
piece_wise_limits = piece_wise_limits[piece_wise_limits[:, 0].argsort()]
return piece_wise_limits
@njit
def _power_curve_drop_consec_duplicates(pc_table):
"""Drop consecutive duplicates in a power curve table.
Parameters
----------
pc_table : np.ndarray
2D array of floats of size (N, 2) power curve table
Returns
-------
np.ndarray:
2D array of floats of size (N-N_duplicates, 2) power curve table.
"""
# Remove any consecutive duplicates from the merged PC
mask_no_duplicates = np.ones(pc_table.shape[0], dtype=np.bool_)
mask_no_duplicates[1:] = pc_table[1:, 0] != pc_table[:-1, 0]
return pc_table[mask_no_duplicates, :]
@njit
def _merge_power_curve_piecewise(piece_wise_limits, pc_table, ws, ws_eff, power):
"""Merge the flow results and power curve piecewise to obtain P(U_free).
This is done for each piece corresponding to monotonically increasing or
decreasing power between local maxima. The sorted results are concatenated.
Parameters
----------
piece_wise_limits : np.ndarray
2D array of integers containing pairs of indicies (the piece wise limits).
pc_table : np.ndarray
2D array of floats of size (N, 2) power curve table
ws : np.ndarray
1D array of floats. Free wind speeds.
ws_eff : np.ndarray
1D array of floats. Effective wind speeds.
power : np.ndarray
1D array of floats. Effective power output.
Returns
-------
np.ndarray:
2D array of floats of size (N, 2). Merged P(U_free) power curve table.
"""
# Create initial array of shape (1, 2) with zero-entries to start
# concatenating from. This allows numba to infer the type of the array.
pc_table_merged = np.array([[0.0, 0.0]], dtype=np.float64)
for i0, i1 in piece_wise_limits:
pc_table_piece = _merge_power_curve_piece(
pc_table, ws[i0 : i1 + 1], ws_eff[i0 : i1 + 1], power[i0 : i1 + 1]
)
pc_table_merged = np.concatenate((pc_table_merged, pc_table_piece), axis=0)
return pc_table_merged
@njit
def _merge_power_curve(pc_table, ws, ws_eff, power):
"""
Original author (of Python code): Tobias Ahsbahs (ttah) from Ole Rathmanns Delphi code
Merges the entries of a power curve with reduced wind speeds and powers.
This is a faster replacement of the one to one code transfer from pascal.
Parameters
----------
pc_table : np.ndarray
2D array of shape (N, M) where M=2 (wind speed, power)
and N is number of P(U) entries.
ws : np.ndarray
1D array with the free wind speed at the turbine location without wake effect
ws_eff : np.ndarray
1D array with the wake-effected wind speed (effective wind speed) at the turbine location
power : np.ndarray
1D array with the wake-effected (reduced) power corresponding to ws_eff
Returns
-------
np.ndarray
2D array of shape (N, M), where M=3 containing [ws, ws_eff, power] columns
and N is the total number of entries in merged power_curve table
Also see:
-------
windkit.WindTurbine.char_data_tables[0][:,[0,2]] as example for list_PC
py_wake.WakeModel.calc_wake for calculating WS, WS_eff, power
Notes
-----
For interpolating the power curve free wind speeds on the wake results the
function WS_eff(WS) needs to invertible. In general this is not true as the
wake effect from upstrem turbine can reduce the effective wind speed even
as the reference wind speed is rising. Therefore, it is necessary to check
if this actually occurs. If so, results from the wake calculation are split
into piecewise invertable functions.
"""
# Local maxima for piecewise definition
ind_max = argrelmax(power)[0]
# Check if there are local maxima in the power curve. If this is true then
# the interpolation needs to be defined piecewise.
if ind_max.size != 0:
piece_wise_limits = _piece_wise_limits(ind_max, power)
pc_table_merged = _merge_power_curve_piecewise(
piece_wise_limits, pc_table, ws, ws_eff, power
)
# If no piecewise definition is needed
else:
# Calculate with all entries. The first element of the merged PC is cut
# off which is the leading 0 introduced above the if statement.
pc_table_merged = _merge_power_curve_piece(pc_table, ws, ws_eff, power)
# Remove any consecutive duplicates from the merged PC
pc_table_merged = _power_curve_drop_consec_duplicates(pc_table_merged)
return pc_table_merged
def _calc_aep_single_turbine_one_wd(
ws, ws_eff, power, A, k, pc_table_index, pc_table_list, *, hours_per_year=8766
):
"""Original author (of Python code): Tobias Ahsbahs (ttah)
Calcualtes the AEP of a single turbine using the power curve and the
reduced wind speed from the wake calculation.
Parameters
----------
ws : np.ndarray
1D array of type float containing free wind speeds
ws_eff : np.ndarray
1D array of type float containing effective wind speeds
power : np.ndarray
1D array of type float containing power output of
the turbine corresponding to ws_eff
A : float
Weibull scale parameter
k : float
Weibull shape parameter
pc_table_index : int
Index of power curve table to use.
pc_table_list : list
List of Array-based power curves tables for each turbine type.
Order in list should correspond to the index in pc_table_index.
hours_per_year : int, optional
Hours in a full year used in AEP calculation(s), by default 8766
Returns
-------
float
AEP in units of Wh
"""
if np.isnan(A) or np.isnan(k):
return np.nan
pc_table = pc_table_list[pc_table_index]
# Merging the power curve with the wake calculations
merged_pc_table = _merge_power_curve(pc_table, ws, ws_eff, power)
# Calculate AEP
aep_single_turb_one_wd, _ = _FORT_POWER_INTEGRAL(
A, k, merged_pc_table[:, 0], merged_pc_table[:, 1]
)
return aep_single_turb_one_wd
def _preprocess_inputs_potential_aep(
wwc, wtg, /, wind_turbines=None, *, shear=None, turbulence_intensity=None
):
"""
Preprocess inputs for the potential_aep and wind_farm_flow_map functions.
The function checks the inputs and converts them to the expected format
for further processing:
- wwc including turbulence_intensity and shear if not present
- wtg's are converted to a dict with keys matching the keys in wind_turbines
- wind_turbines are converted to a xr.Dataset if not already
Parameters
----------
wwc : xr.Dataset
PyWAsP Weibull Wind Climate Dataset containing
wind climates for all the turbine locations.
wtg : xr.Dataset, dict, wk.WindTurbines
Wind Turbine Generator as a xr.Dataset, a dict of wtg_keys and xr.Datasets, or
a wk.WindTurbines object containing groups of WTGs and locations.
If a xr.Dataset is passed, the WWC locations are assumed to be the same
as the turbine locations. wk.WindTurbines is deprecated and will be removed in a future version.
wind_turbines : xr.Dataset, optional
Wind Turbine locations, hub heights, group_id's, and wtg_keys.
shear : float, optional
Shear exponent to use in the AEP calculation.
turbulence_intensity : float, optional
Turbulence intensity to use in the AEP calculation.
Returns
-------
xr.Dataset, dict, xr.Dataset
Updated wwc, wtg_dict, wind_turbines
"""
# We always convert to wtg_dict and wind turbines xr.Dataset
# If wtg is a xr.Dataset, we assume convert it to a dict and add the key to
# the wind turbines xr.Dataset
if isinstance(wtg, xr.Dataset):
wtg_dict = {str(wtg.name.data): wtg}
if wind_turbines is None:
wind_turbines = wk.create_wind_turbines_from_arrays(
west_east=wwc.west_east.data,
south_north=wwc.south_north.data,
height=wwc.height.data,
crs=wk.spatial.get_crs(wwc),
group_ids=np.zeros(wwc.point.size, dtype=int),
turbine_ids=np.arange(wwc.point.size),
wtg_keys=[str(wtg.name.data)] * wwc.point.size,
)
else:
wind_turbines = wind_turbines.copy()
wind_turbines["wtg_key"] = (("point",), [wtg.name.data])
elif isinstance(wtg, wk.WindTurbines):
warnings.warn(
"The WindTurbins class is deprecated and will be removed in a future version. "
+ "Use a wtg dict and xr.Dataset of wind turbines instead. ",
DeprecationWarning,
stacklevel=2,
)
if wind_turbines is not None:
raise ValueError(
"When passing a wk.WindTurbines object as 'wtg', 'wind_turbines' must not be passed."
)
wtg_dict, wind_turbines = wtg.to_wtg_dict_and_turbines_ds()
elif isinstance(wtg, dict):
if wind_turbines is None:
raise ValueError(
"When passing a dict as 'wtg', 'wind_turbines' must also be passed."
)
wtg_dict = wtg
else:
raise TypeError(
f"wtg must be a xr.Dataset, wk.WindTurbines object, or dict, not {type(wtg)}"
)
if turbulence_intensity is not None:
if "turbulence_intensity" not in wwc.data_vars:
wwc["turbulence_intensity"] = xr.full_like(wwc.A, turbulence_intensity)
else:
warnings.warn(
"turbulence_intensity is already in the WWC but also passed explicit. It will be overwritten."
)
elif "turbulence_intensity" not in wwc.data_vars:
raise ValueError(
"turbulence_intensity is not in the WWC and no explicit value was passed."
)
if shear is not None:
if "shear" not in wwc.data_vars:
wwc["shear"] = xr.full_like(wwc.A, shear)
else:
warnings.warn(
"shear is already in the WWC but also passed explicit. It will be overwritten."
)
return wwc, wtg_dict, wind_turbines
[docs]
@wk.spatial.stack_then_unstack
def potential_aep(
wwc,
wtg,
/,
wind_turbines=None,
*,
wind_farm_model="PARK2_onshore",
ws_stepsize=1.0,
ws_upper_limit=None,
ws_lower_limit=0.0,
n_subsector=5,
site_interp_method="nearest",
site_interp_bounds="ignore",
turbulence_intensity=None,
shear=None,
n_cpu_pywake=1,
):
"""Calculate Potential Annual Energy Production
using wind farm effects from PyWake.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
wwc : xr.Dataset
PyWAsP Weibull Wind Climate Dataset containing
wind climates for all the turbine locations.
wtg : xr.Dataset, dict, wk.WindTurbines
- If a xr.Dataset is passed it should be a WTG formatted dataset.
PyWAsP will assume WWC points match the turbines hub height and that the
WTG is the same for all turbines
- If a dict of wtg_keys and WTG xr.Datasets is passed, a ``turbines``
dataset must be passed as well, and the key's must match between the two.
- A wk.WindTurbines object can also be passed, but is deprecated and
will be removed in a future version.
wind_turbines : xr.Dataset, optional
Wind Turbine locations, hub heights, group_id's, and wtg_keys.
wind_farm_model : str, function
Wind farm model to use for deficit calculations. Can either be a name of
a predefined wind farm model: "PARK1", "PARK2_onshore", "PARK2_offshore", or a
predefined py_wake wind_farm_model object (wrapped in a 'functools.partial' function)
ws_stepsize : float, optional
Wind speed step sizes between wind speed cases, by default 1.0 m/s
n_subsector : int, optional
Number of subsectors to use per sector, by default None
Which means 5 subsectors will be used.
ws_upper_limit : float, optional
Upper limit of wind speed range, by default None
Which means the upper limit of the wind speed range will be the maximum
ws_lower_limit : float, optional
Lower limit of wind speed range, by default 0.0
n_subsector : int, optional
Number of subsectors to use per sector, by default None
Which means 5 subsectors will be used.
site_interp_method : str, optional
The interpolation method to use in py_wake XRSite, by default "nearest"
site_interp_bounds : str, optional
How to handle values outside the site boundary, by default "ignore"
Which means values outside the site bondary will not be checked for.
turbulence_intensity : float, optional
Turbulence intensity to use for the AEP calculation. If not passed, the turbulence intensity
from the WWC will be used. If turbulence intensity is passed, it will be used,
and overwrite the turbulence intensity in the WWC, if present.
shear : float, optional
shear to use for the AEP calculation. If not passed, the shear
from the WWC will be used. If shear is passed, it will be used,
and overwrite the shear in the WWC, if present.
n_cpu_pywake : int, optional
Number of CPUs to use for the PyWake calculations, by default 1
If n_cpu_pywake is None, the number of CPUs will be set to the number of
CPUs available on the machine.
Returns
-------
xr.Dataset
Potential Annual Energy Production in GWh for each turbine/location and sector.
Raises
------
ModuleNotFoundError
If PyWake is not installed.
TypeError
If wtg is not a xr.Dataset or wk.WindTurbines object.
ValueError
If required variables are not in the WWC.
ValueError
If turbulence_intensity is not in the WWC and no explicit value was passed.
"""
if py_wake is None:
raise ModuleNotFoundError(
"PyWake not installed. It is required for AEP calculations with wakes!"
)
wwc = wwc.copy()
wwc, wtg_dict, wind_turbines = _preprocess_inputs_potential_aep(
wwc, wtg, wind_turbines, turbulence_intensity=turbulence_intensity, shear=shear
)
flow = _run_pywake_flowcases(
wwc,
wtg_dict,
wind_turbines,
output_locs=None,
wind_farm_model=wind_farm_model,
ws_stepsize=ws_stepsize,
ws_upper_limit=ws_upper_limit,
ws_lower_limit=ws_lower_limit,
n_subsector=n_subsector,
site_interp_method=site_interp_method,
site_interp_bounds=site_interp_bounds,
n_cpu_pywake=n_cpu_pywake,
)
flow = _postproc_pywake_flow_results(flow, wtg_dict, wind_turbines)
# Calculate AEP and aggregate subsectors
flow["potential_AEP_sector"] = (
_calculate_aep_wasp_method(flow, wtg_dict) * flow["wdfreq"]
)
# Calculate the AEP deficit by sector
flow["AEP_deficit_sector"] = -(
(flow["potential_AEP_sector"] / flow["gross_AEP_sector"]) - 1.0
)
# Aggregate subsectors -> sectors
flow = _agg_subsectors(
flow[
[
"potential_AEP_sector",
"gross_AEP_sector",
"AEP_deficit_sector",
"wspd_sector",
"wspd_eff_sector",
"wspd_deficit_sector",
"turbulence_intensity_eff_sector",
"wdfreq",
]
]
)
# Calculate the average wind speed, AEP, and deficit values over sectors
flow["potential_AEP"] = flow["potential_AEP_sector"].sum(dim="sector")
flow["gross_AEP"] = flow["gross_AEP_sector"].sum(dim="sector")
flow["AEP_deficit"] = (flow["wdfreq"] * flow["AEP_deficit_sector"]).sum(
dim="sector"
)
flow["wspd"] = (flow["wdfreq"] * flow["wspd_sector"]).sum(dim="sector")
flow["wspd_eff"] = (flow["wdfreq"] * flow["wspd_eff_sector"]).sum(dim="sector")
flow["wspd_deficit"] = (flow["wdfreq"] * flow["wspd_deficit_sector"]).sum(
dim="sector"
)
flow["turbulence_intensity_eff"] = (
flow["wdfreq"] * flow["turbulence_intensity_eff_sector"]
).sum(dim="sector")
flow = update_var_attrs(flow, _AEP_ATTRS)
flow = update_history(flow)
return flow
def _split_wwc_into_subsectors(
wwc,
n_subsector=5,
dim_sector_and_subsector="sector_and_subsector",
dim_subsector="subsector",
dim_sector="sector",
):
"""Split the sectors of a weibull wind climate into
a number of subsectors. The sectors and subsectors are stacked
onto a new dimension that combines both.
Parameters
----------
wwc : xr.Dataset
PyWAsP Weibull Wind Climate Dataset containing
wind climates for all the turbine locations.
n_subsector : int, optional
number of subsectors to split into, by default 5
dim_sector_and_subsector : str, optional
name of combined sector+subsector dimension, by default "sector_and_subsector"
dim_subsector : str, optional
name of subsector, by default "subsector"
dim_sector : str, optional
name of sector, by default "sector"
Returns
-------
xr.Dataset
Weibull wind climate with sectors split into subsectors.
"""
wwc_split = wwc.copy()
sector_width = 360.0 / wwc_split[dim_sector].size
subsector_width = sector_width / n_subsector
subsector_relative_center = np.linspace(
-subsector_width * (n_subsector // 2),
subsector_width * (n_subsector // 2),
n_subsector,
)
wwc_split = wwc_split.expand_dims(
**{dim_subsector: subsector_relative_center}
).stack(**{dim_sector_and_subsector: (dim_sector, dim_subsector)})
return wwc_split
def _get_ws_flowcases(
wwc,
wtg_dict,
ws_stepsize=1.0,
ws_upper_limit=None,
ws_lower_limit=None,
):
"""Get the needed flowcases for the wwc and wtg combination.
Parameters
----------
wwc : xr.Dataset
Weibull Wind Climate Dataset containing
wind climates for all the turbine locations.
wtg_dict : dict
Dictionary of Wind Turbine Generator xr.Datasets
ws_stepsize : float, optional
Wind speed step sizes between wind speed cases, by default 1.0
ws_upper_limit: float, optional
Maximum value of the wind speed cases, by default no max is set.
ws_lower_limit: float, optional
Minimum value (starting point) of the wind speed cases, be default
no minimum is set and the cut-in wind speed is used. However, to
obtain accurate mean wind speed values during a later integration, the
bins should start at 0.0.
Returns
-------
np.ndarray
1D array of (free) wind speed cases to simulate.
"""
# Find wind speed values to simulate
ws_cutin = min(wtg.wind_speed_cutin.values[0] for wtg in wtg_dict.values())
ws_cutout = max(wtg.wind_speed_cutout.values[0] for wtg in wtg_dict.values())
ws_mean = wk.weibull_moment(wwc["A"], wwc["k"], n=1.0)
spatial_dims = wk.spatial._struct._get_xy_dims(wwc)
ws_sectorwise_max = ws_mean.max(dim=spatial_dims)
ws_sectorwise_min = ws_mean.min(dim=spatial_dims)
ws_speedup_max = (ws_sectorwise_max / ws_sectorwise_min).max()
ws_min = ws_cutin
ws_max = ws_speedup_max * ws_cutout
if ws_lower_limit is None:
ws = np.arange(ws_min, ws_max + 2 * ws_stepsize, ws_stepsize)
else:
ws = np.arange(0.0, ws_max + 2 * ws_stepsize, ws_stepsize)
if ws_upper_limit is not None:
ws = ws[ws <= ws_upper_limit]
return ws
[docs]
def wind_farm_flow_map(
wwc,
wtg_dict,
wind_turbines,
output_locs,
wind_farm_model="PARK2_onshore",
turbulence_intensity=None,
shear=None,
ws_stepsize=1.0,
ws_upper_limit=None,
ws_lower_limit=0.0,
n_subsector=5,
wtg_id=0,
site_interp_method="nearest",
site_interp_bounds="ignore",
n_cpu_pywake=1,
):
"""Generate a flow map around a wind farm using py_wake for wake and
blockage effects. The map is calculated at the same points as the
weibull wind climate(s) are provided.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
wwc : xr.Dataset
Weibull Wind Climate Cuboid Dataset containing
wind climates for all the turbine locations.
wtg_dict : dict
A dict of wtg_keys and WTG xr.Datasets should be passed, a ``turbines``
dataset must be passed as well, and the key's must match between the two.
wind_turbines : xr.Dataset
Wind Turbine locations, hub heights, group_id's, and wtg_keys.
output_locs : xr.Dataset
Locations to calculate the flow map at. Must be a "cuboid" xr.Dataset
with 3 spatial dimensions: 'west_east', 'south_north', 'height'.
The coordinates of the output locations must be covered by the
cuboid spanned by the weibull wind climate(s).
wind_farm_model : str, function, optional
Wind farm model to use for deficit calculations. Can
Either be a name of a predefined wind farm model: "PARK1", "PARK2_onshore",
"PARK2_offshore", or a predefined py_wake wind_farm_model
object. By default "PARK2_onshore" is used.
turbulence_intensity : float, optional
Turbulence intensity to use for the AEP calculation. If not passed, the turbulence intensity
from the WWC will be used. If turbulence intensity is passed, it will be used,
and overwrite the turbulence intensity in the WWC, if present.
shear : float, optional
shear to use for the AEP calculation. If not passed, the shear
from the WWC will be used. If shear is passed, it will be used,
and overwrite the shear in the WWC, if present.
ws_stepsize : float, optional
Wind speed bins stepsize, by default 1.0
ws_upper_limit : float, optional
Upper limit of wind speed range, by default None
Which means the upper limit of the wind speed range will be the maximum
ws_lower_limit : float, optional
Lower limit of wind speed range, by default 0.0
n_subsector : int, optional
Number of subsectors to use per sector, by default 5
wtg_id : int, optional
The wind turbine type to use for calculating the AEP flow map, by default 0, the first type
The WTG is taken from the corresponding entry in the wind_turbines object.
site_interp_method : str, optional
The interpolation method to use in py_wake XRSite, by default "nearest"
site_interp_bounds : str, optional
The extrapolation method to use in py_wake XRSite, by default "ignore"
n_cpu_pywake : int, optional
Number of CPUs to use in py_wake, by default 1
Returns
-------
xr.Dataset
Wind farm flow map for the output_locs containing the variables:
- potential_AEP_sector:
Potential AEP in GWh for each sector
- gross_AEP_sector:
Gross AEP in GWh for each sector
- AEP_deficit_sector:
AEP deficit in units of fraction for each sector
- wspd_sector:
Average wind speed in m/s for each sector
- wspd_eff_sector:
Average effective wind speed in m/s for each sector
- wspd_deficit_sector:
Average wind speed deficit in units of fraction for each sector
- turbulence_intensity_eff_sector:
Average effective turbulence intensity for each sector
- wdfreq:
Wind direction frequency in units of fraction for each sector
- potential_AEP:
Potential AEP in GWh for each location
- gross_AEP:
Gross AEP in GWh for each location
- AEP_deficit:
AEP deficit in units of fraction for each location
- wspd:
Average wind speed in m/s for each location
- wspd_eff:
Average effective wind speed in m/s for each location
- wspd_deficit:
Average wind speed deficit in units of fraction for each location
- turbulence_intensity_eff:
Average effective turbulence intensity for each location
Raises
------
TypeError
If output_locs is not a xr.Dataset
ValueError
If output_locs is not a "cuboid" xr.Dataset
ValueError
If wwc is not a "cuboid" xr.Dataset
ValueError
If output_locs is not covered by the cuboid spanned by the weibull wind climate(s)
ValueError
If the wind turbine locations are not covered by the cuboid spanned by the weibull wind climate(s)
"""
if not isinstance(output_locs, xr.Dataset):
raise TypeError(f"output_locs must be a xr.Dataset, not {type(output_locs)}")
if not wk.spatial.is_cuboid(output_locs):
raise ValueError(
"output_locs must be a 'cuboid' (i.e. have 3 spatial dimensions: 'west_east', 'south_north', 'height')"
)
if not wk.spatial.is_cuboid(wwc):
raise ValueError(
"wwc must be a 'cuboid' (i.e. have 3 spatial dimensions: 'west_east', 'south_north', 'height')"
)
if not wk.spatial.covers(wwc, output_locs):
raise ValueError(
"The output_locs must be covered by the cuboid spanned by the weibull wind climate(s)"
)
if not wk.spatial.covers(wwc, wind_turbines):
raise ValueError(
"The wind turbine locations must be covered by the cuboid spanned by the weibull wind climate(s)"
)
wwc, wtg_dict, wind_turbines = _preprocess_inputs_potential_aep(
wwc,
wtg_dict,
wind_turbines,
turbulence_intensity=turbulence_intensity,
shear=shear,
)
flow = _run_pywake_flowcases(
wwc,
wtg_dict,
wind_turbines,
output_locs=output_locs,
wind_farm_model=wind_farm_model,
ws_stepsize=ws_stepsize,
ws_upper_limit=ws_upper_limit,
ws_lower_limit=ws_lower_limit,
n_subsector=n_subsector,
wtg_id=wtg_id,
site_interp_method=site_interp_method,
site_interp_bounds=site_interp_bounds,
n_cpu_pywake=n_cpu_pywake,
)
flow = _postproc_pywake_flow_results(flow, wtg_dict, wind_turbines)
flow["potential_AEP_sector"] = (
_calculate_aep_wasp_method(flow, wtg_dict) * flow["wdfreq"]
)
# Calculate the AEP deficit by sector
flow["AEP_deficit_sector"] = -(
(flow["potential_AEP_sector"] / flow["gross_AEP_sector"]) - 1.0
)
# Aggregate subsectors -> sectors
flow = _agg_subsectors(
flow[
[
"potential_AEP_sector",
"gross_AEP_sector",
"AEP_deficit_sector",
"wspd_sector",
"wspd_eff_sector",
"wspd_deficit_sector",
"turbulence_intensity_eff_sector",
"wdfreq",
]
]
)
# Calculate the average wind speed, AEP, and deficit values over sectors
flow["potential_AEP"] = flow["potential_AEP_sector"].sum(dim="sector")
flow["gross_AEP"] = flow["gross_AEP_sector"].sum(dim="sector")
flow["AEP_deficit"] = (flow["wdfreq"] * flow["AEP_deficit_sector"]).sum(
dim="sector"
)
flow["wspd"] = (flow["wdfreq"] * flow["wspd_sector"]).sum(dim="sector")
flow["wspd_eff"] = (flow["wdfreq"] * flow["wspd_eff_sector"]).sum(dim="sector")
flow["wspd_deficit"] = (flow["wdfreq"] * flow["wspd_deficit_sector"]).sum(
dim="sector"
)
flow["turbulence_intensity_eff"] = (
flow["wdfreq"] * flow["turbulence_intensity_eff_sector"]
).sum(dim="sector")
flow = flow.transpose(*wwc.dims, ...)
flow = update_var_attrs(flow, _WF_FLOW_MAP_ATTRS)
flow = update_history(flow)
return flow
def _drop_sector_vars(ds):
"""
Drop sector variables from a dataset.
Parameters
----------
ds : xr.Dataset
Dataset to drop sector variables from.
Returns
-------
xr.Dataset
Dataset with sector variables dropped.
"""
sector_vars = (
"sector",
"subsector",
"sector_and_subsector",
"sector_ceil",
"sector_floor",
)
return ds.drop_vars([c for c in sector_vars if c in ds.coords])
def _run_pywake_flowcases(
wwc,
wtg_dict,
wind_turbines,
output_locs=None,
wind_farm_model="PARK2_onshore",
turbulence_intensity=None,
shear=None,
ws_stepsize=1.0,
ws_lower_limit=None,
ws_upper_limit=None,
n_subsector=None,
wtg_id=0,
site_interp_method="nearest",
site_interp_bounds="ignore",
n_cpu_pywake=1,
):
"""Run wind farm flow cases using PyWake.
Parameters
----------
wwc : xr.Dataset
Weibull Wind Climate Dataset containing
wind climates for all the turbine locations.
wtg_dict : dict
Dictionary of Wind Turbine Generator xr.Datasets
wind_turbines : xr.Dataset
WindTurbines object containing turbine locations groups and wtg keys.
output_locs : xr.Dataset, optional
Locations to calculate the flow map at. Must be a "cuboid" xr.Dataset
wind_farm_model : str, function, optional
Wind farm model to use for deficit calculations. Can
Either be a name of a predefined wind farm model:
"PARK", "PARK2_onshore", "PARK2_offshore", or a predefined py_wake
wind_farm_model object be default "PARK2_onshore" is used.
turbulence_intensity : float, optional
Turbulence intensity to use for the AEP calculation. If not passed, the turbulence intensity
from the WWC will be used. If turbulence intensity is passed, it will be used,
and overwrite the turbulence intensity in the WWC, if present.
shear : float, optional
shear to use for the AEP calculation. If not passed, the shear
from the WWC will be used. If shear is passed, it will be used,
and overwrite the shear in the WWC, if present.
ws_stepsize : float, optional
Wind speed bins stepsize, by default 1.0
ws_upper_limit : float, optional
Upper limit of wind speed range, by default None
Which means the upper limit of the wind speed range will be the maximum
ws_lower_limit : float, optional
Lower limit of wind speed range, by default 0.0
n_subsector : int, optional
Number of subsectors to use per sector, by default None
Which means 5 subsectors will be used.
wtg_id : int, optional
The wind turbine type to use for calculating the AEP flow map, by default 0, the first type
The WTG is taken from the corresponding entry in the wind_turbines object.
site_interp_method : str, optional
The interpolation method to use in py_wake XRSite, by default "nearest"
site_interp_boudns : str, optional
The boundary condition to use in py_wake XRSite, by default "ignore"
n_cpu_pywake : int, optional
Number of CPUs to use in py_wake, by default 1
If n_cpu_pywake is None, the number of CPUs will be set to the number of
CPUs available on the machine.
Returns
-------
xr.Dataset
Wind farm flow simulation results merged with the weibull wind climate.
Raises
------
ModuleNotFoundError
Raised if 'py_wake' is not installed in environment.
ValueError
Raised if 'wind_farm_model' is not string or partial function
"""
if py_wake is None:
raise ModuleNotFoundError(
"PyWake not installed. It is required for AEP calculations with wakes!"
)
wind_turbines_pywake = py_wake.wind_turbines.WindTurbines.from_WindTurbine_lst(
[wtg_to_pywake(wtg) for wtg in wtg_dict.values()]
)
type_indx = [
list(wtg_dict.keys()).index(wtg_key) for wtg_key in wind_turbines.wtg_key.values
]
# Translate WWC to pywake XRSite format
site_pywake = py_wake.site.xrsite.XRSite.from_pwc(
wwc, interp_method=site_interp_method, bounds=site_interp_bounds
)
# Fetch/check wind farm model
if isinstance(wind_farm_model, str):
wind_farm_model = _wind_farm_model_from_string(wind_farm_model)
elif not isinstance(wind_farm_model, partial):
raise ValueError("wind_farm_model must be a string or py_wake wind farm model")
# Add site and wtg data to wind farm model
wind_farm = wind_farm_model(site_pywake, wind_turbines_pywake)
# Make sure no "mutable" effects leaves scope
wwc = wwc.copy()
if n_subsector is not None:
wwc = _split_wwc_into_subsectors(wwc, n_subsector=n_subsector)
# Find wind direction values to simulate
wd_cases = np.mod((wwc["sector"] + wwc["subsector"]).values, 360.0)
dim_sector = "sector_and_subsector"
else:
wd_cases = wwc["sector"].values
dim_sector = "sector"
ws_cases = _get_ws_flowcases(
wwc,
wtg_dict=wtg_dict,
ws_stepsize=ws_stepsize,
ws_lower_limit=ws_lower_limit,
ws_upper_limit=ws_upper_limit,
)
# Run the py_wake flow cases
sim_res = wind_farm(
x=wind_turbines.coords["west_east"].values,
y=wind_turbines.coords["south_north"].values,
h=wind_turbines.coords["height"].values,
type=type_indx,
wd=wd_cases,
ws=ws_cases,
n_cpu=n_cpu_pywake,
)
spatial_struct = wk.spatial.get_spatial_struct(wwc)
if spatial_struct == "cuboid":
horizontal_grids = []
for output_height in output_locs["height"].values:
horizontal_grid = py_wake.flow_map.HorizontalGrid(
x=output_locs["west_east"].values,
y=output_locs["south_north"].values,
h=output_height,
)
# Calculate the flow around the farm
flow = sim_res.flow_map(horizontal_grid)
flow["Power"] = flow.power_xylk(type=wtg_id)
# For flow maps the first wtg entry is used to calculate the AEP map
flow["wtg_type"] = (flow["WS"].isel(ws=0, wd=0) * wtg_id).astype(int)
flow = xr.Dataset(flow)
horizontal_grids.append(flow)
flow = xr.concat(horizontal_grids, dim="h")
wwc = wwc.reindex_like(output_locs, method="nearest")
else:
flow = sim_res.sel(wd=wd_cases)
flow = flow[["WS", "WS_eff", "Power", "TI_eff"]]
flow["wtg_type"] = (("i",), type_indx)
for var in [v for v in flow.data_vars if "wt" in flow[v].dims]:
flow[var] = flow[var].rename({"wt": "i"})
flow = xr.Dataset(flow)
# Rename to pywasp names
new_names = {
"x": "west_east",
"y": "south_north",
"h": "height",
"wd": dim_sector,
"ws": "wind_speed",
"i": "point",
}
flow = flow.rename({k: v for k, v in new_names.items() if k in flow})
flow = flow.transpose(*wwc.dims, ...)
flow = _drop_sector_vars(flow)
flow = wwc.merge(flow, compat="override")
return flow
def _calculate_aep_wasp_method(flow_results, wtg_dict):
"""
Caculate AEP using the WAsP methodology.
1. Caculate WS_eff and power results for a matrix
of free stream wind speeds and directions. Turbine
with the highest mean wind speed by direction is used as the
referene turbine and all other turbines are scaled according
to the relative speed-up.
2. Merge simulation results and the power curve to obtain
the power curve in terms of free stream wind speed
as opposed to the local wind speed.
3. Integrate power over wind speed for all the turbines and directions.
Parameters
----------
flow_results : xr.Dataset
Flow results from the PyWake simulations results obtained from
the "_run_pywake_flowcases" function .
wtg_dict : dict
Dictionary of Wind Turbine Generator xr.Datasets
Returns
-------
xr.DataArray
AEP (potential) in GWh for each sector
"""
pc_table_list = [
np.stack([wtg.wind_speed.values, wtg["power_output"].values[0, :]], axis=1)
for wtg in wtg_dict.values()
]
# # Calculate AEP by subsector and position
flow_results["potential_AEP_sector"] = xr.apply_ufunc(
_calc_aep_single_turbine_one_wd,
flow_results["WS"],
flow_results["WS_eff"],
flow_results["Power"],
flow_results["A"],
flow_results["k"],
flow_results["wtg_type"],
input_core_dims=[["wind_speed"], ["wind_speed"], ["wind_speed"], [], [], []],
output_core_dims=[[]],
kwargs={"pc_table_list": pc_table_list},
vectorize=True,
dask="parallelized",
output_dtypes=["float64"],
)
return flow_results["potential_AEP_sector"] * 1e-9
def _agg_subsectors(
ds, dim_sector_and_subsector="sector_and_subsector", dim_subsector="subsector"
):
"""Unstack subsectors from sectors and compute the average across subsectors.
Parameters
----------
ds : xr.Dataset
Dataset
dim_sector_and_subsector : str, optional
stacked sectors and subsectors dimension, by default "sector_and_subsector"
dim_subsector : str, optional
subsector dimension, by default "subsector"
Returns
-------
xr.Dataset
Dataset with unstacked and averaged subsectors
"""
if dim_sector_and_subsector in ds.dims:
# Need this to keep attrs
with xr.set_options(keep_attrs=True):
# Unstack sectors and subsectors, and average subsector results
ds = ds.unstack(dim_sector_and_subsector)
ds = ds.mean(dim=dim_subsector)
return ds
# @timeit
def _postproc_pywake_flow_results(
flow_results,
wtg_dict,
wind_turbines,
dim_wind_speed="wind_speed",
hours_per_year=8766,
wtg_id=0,
):
"""Calculate mean wind speed and mean effective wind speed from
PyWake flow results obtained from the "_run_pywake_flowcases" function above
it should cover the total range of wind speed (bins) in the
possible probability space (including speed-ups).
Parameters
----------
flow_results : xr.Dataset
Flow results from the PyWake simulations results obtained from
the "_run_pywake_flowcases" function above.
wtg_dict : dict
Dictionary of Wind Turbine Generator xr.Datasets
wind_turbines : xr.Dataset
WindTurbines object containing turbine locations, id's, group id's, and wtg keys mapped to
wtg_dict.
dim_wind_speed : str, optional
The dimension of wind speed bins in the simulation results,
by default "wind_speed"
hours_per_year : int, optional
Number of hours per year, by default 8766
wtg_id : int, optional
The wind turbine type to use for calculating the AEP flow map, by default 0, the first type
Returns
-------
xr.Dataset:
Flow simulation results with mean wind speed
and mean effective wind speed addded.
"""
def _calc_prob(A, k, ws_lower, ws_upper, dim="wind_speed"):
"""Calculate probability mass based on Weibull distributions
and wind speed bins.
#TODO: integrate/move to windkit.weibull
Parameters
----------
A : xr.DataArray
Weibull A parameter
k : xr.DataArray
Weibull k parameter
ws_lower : xr.DataArray
Lower edge value of the wind speed bins
ws_upper : xr.DataArray
Upper edge value of the wind speed bins
dim : str, optional
Dimension to sum over, by default "wind_speed"
Returns
-------
xr.DataArray
Bin probabilities
"""
def _weibull_pmf(A, k, ws_lower, ws_upper):
# Expand dims to match rank of ws_lower/ws_upper
A = A[..., None]
k = k[..., None]
cdf_l = 1.0 - np.exp(-((ws_lower / A) ** k))
cdf_u = 1.0 - np.exp(-((ws_upper / A) ** k))
pmf = cdf_u - cdf_l
return pmf
return xr.apply_ufunc(
_weibull_pmf,
A,
k,
ws_lower,
ws_upper,
input_core_dims=[[], [], [dim], [dim]],
output_core_dims=[[dim]],
dask="parallelized",
output_dtypes=["float64"],
)
def _calc_average(lower, upper, prob, dim="wind_speed"):
"""Caculate average by summing up values weigthed by
probability density over binned wind speed range.
The average value in each bin is used as the representative value.
Parameters
----------
lower : xr.DataArray
Lower edge value of the quantity
upper : xr.DataArray
Upper edge value of the quantity
prob : xr.DataArray
Bin probabilities
dim : str, optional
Dimension to sum over, by default "wind_speed"
Returns
-------
xr.DataArray
Average value over wind speed bins.
"""
return (0.5 * (lower + upper) * prob).sum(dim=dim)
n_wtgs = len(wtg_dict.keys())
wtgs = list(wtg_dict.values())
flow_results["gross_power"] = wk.wtg_power(wtgs[0], flow_results["WS"])
if n_wtgs > 1:
for iwtg in range(1, n_wtgs):
flow_results["gross_power"] = xr.where(
flow_results["wtg_type"] == iwtg,
wk.wtg_power(wtgs[iwtg], flow_results["WS"]),
flow_results["gross_power"],
)
fr_lower = flow_results.drop_vars(dim_wind_speed).isel(
**{dim_wind_speed: slice(None, -1)}
)
fr_upper = flow_results.drop_vars(dim_wind_speed).isel(
**{dim_wind_speed: slice(1, None)}
)
# calculate probabilty mass of each bin in the flow results and mean wind speed(s)
ws_lower = fr_lower["WS"]
ws_upper = fr_upper["WS"]
prob = _calc_prob(
flow_results["A"], flow_results["k"], ws_lower, ws_upper, dim=dim_wind_speed
)
# Calculate mean wind speed
flow_results["wspd_sector"] = _calc_average(
ws_lower, ws_upper, prob, dim=dim_wind_speed
)
# Caculate effective mean wind speed(s)
flow_results["wspd_eff_sector"] = _calc_average(
fr_lower["WS_eff"],
fr_upper["WS_eff"],
prob,
dim=dim_wind_speed,
)
# Calculate wind speed defifict (positive = more deficit)
flow_results["wspd_deficit_sector"] = -(
(flow_results["wspd_eff_sector"] / flow_results["wspd_sector"]) - 1.0
)
# Caculate gross AEP
flow_results["gross_AEP_sector"] = (
_calc_average(
fr_lower["gross_power"],
fr_upper["gross_power"],
prob,
dim=dim_wind_speed,
)
* 1e-9
* flow_results["wdfreq"]
* hours_per_year
)
# Caculate effective mean wind speed(s) (probability mass is the same as above)
flow_results["turbulence_intensity_eff_sector"] = _calc_average(
fr_lower["TI_eff"], fr_upper["TI_eff"], prob, dim=dim_wind_speed
)
return flow_results
[docs]
def wtg_to_pywake(wtg):
"""Converts a PyWAsP WTG to a PyWake WTG
Parameters
----------
wtg : xr.Dataset
PyWAsP Wind Turbine Generator (WTG) Dataset
Returns
-------
WindTurbine
PyWake WindTurbine object
Raises
------
ModuleNotFoundError
Raised if 'py_wake' is not installed in environment.
"""
if py_wake is None:
raise ModuleNotFoundError("This function requires 'py_wake' to run!")
wtg_renamed = (
wtg[["power_output", "thrust_coefficient", "air_density"]]
.rename(
{
"power_output": "power",
"thrust_coefficient": "ct",
"air_density": "rho",
"wind_speed": "ws",
}
)
.transpose("ws", "mode", missing_dims="warn")
)
if wtg.coords["mode"].size > 1:
wtg_renamed = wtg_renamed.swap_dims(mode="rho").sortby("rho").drop_vars("mode")
power_ct_curves = py_wake.wind_turbines.power_ct_functions.PowerCtXr(
wtg_renamed, "W"
)
else:
wspd = wtg_renamed["ws"].values
power = wtg_renamed["power"].values.squeeze()
ct = wtg_renamed["ct"].values.squeeze()
power_ct_curves = py_wake.wind_turbines.power_ct_functions.PowerCtTabular(
wspd,
power,
"W",
ct,
ws_cutin=wtg.wind_speed_cutin.values,
ws_cutout=wtg.wind_speed_cutout.values,
)
return py_wake.wind_turbines.WindTurbine(
name=wtg.name.values,
diameter=wtg.rotor_diameter.values,
hub_height=wtg.hub_height.values,
powerCtFunction=power_ct_curves,
ws_cutin=wtg.wind_speed_cutin.values,
ws_cutout=wtg.wind_speed_cutout.values,
)