Source code for pywasp.wasp.aep

"""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, )