.. py:currentmodule:: pywasp.wasp

.. _calculate_aep:

=============
Calculate AEP
=============

PyWAsP can estimate the Annual Energy Production (AEP) of wind turbines.

For now, this includes estimating the gross and potential AEP (as per the workflow below),
with only limited support, for now, for estimating other losses and P50 and P90.

.. figure:: ../_static/yield_workflow.png
    :align: center

    The yield assessment workflow, from :cite:`mortensen2015offshore`

We will start by importing needed libraries and some test data from the
Serra Santa Luzia site.

.. ipython:: python
    :okwarning:

    import numpy as np
    import pandas as pd
    import xarray as xr
    import matplotlib.pyplot as plt

    import windkit as wk
    import pywasp as pw

.. ipython:: python
    :okwarning:

    bwc = wk.read_bwc("./source/tutorials/data/SerraSantaLuzia.omwc", "EPSG:4326")
    bwc = wk.spatial.reproject(bwc, "EPSG:32629")

    elev_map = wk.read_vector_map(
        "./source/tutorials/data/SerraSantaLuzia.map",
        crs="EPSG:32629",
        map_type="elevation"
    )

    lc_map, lc_tbl = wk.read_vector_map(
        "./source/tutorials/data/SerraSantaLuzia.map",
        crs="EPSG:32629",
        map_type="roughness"
    )

    topo_map = pw.wasp.TopographyMap(elev_map, lc_map, lc_tbl)

We have used the site in the :ref:`wasp_flow_model` page, but here we will
go further and include some turbines. Their locations are stored in a ``.csv``,
so we will create a `xarray.Dataset` to hold their spatial information: x, y, and hub height.

We can plot their locations together with the elevation map to get an idea of where they are located:

.. ipython:: python
    :okwarning:

    wtg_locs = pd.read_csv('./source/tutorials/data/turbine_positions.csv')

    output_locs =  wk.spatial.create_dataset(
        wtg_locs.Easting.values,
        wtg_locs.Northing.values,
        wtg_locs['Hub height'].values,
        crs="EPSG:32629"
    )

    fig, ax = plt.subplots()
    ax.set(xlim=(510000, 521500), ylim=(4616000, 4627500))
    elev_map.plot("elev", ax=ax);
    @savefig turbine_positions.png align=center width=6in
    ax.scatter(
        output_locs.west_east.values,
        output_locs.south_north.values,
        c="k",
        marker="x",
        s=15,
        zorder=2,
    );

To estimate the AEP, we will fist estimate wind climate at the hub height of each turbine:

.. ipython:: python
    :okwarning:

    wwc = pw.wasp.predict_wwc(bwc, topo_map, output_locs)
    print(wwc)


Gross AEP
=========

Before we estiamte the AEP we have to get a Wind Turbine Generator (WTG) corresponding to the turbine models used.
We use :py:func:`windkit.read_wtg` to read the ``.wtg`` file.

.. ipython:: python
    :okwarning:

    wtg = wk.read_wtg("./source/tutorials/data/Bonus_1_MW.wtg")
    print(wtg)

The wtg holds power and thrust curves for a number of wind speeds and operational modes (in this case just 1).
Now we can estimate the gross AEP using :py:func:`gross_aep`:

.. ipython:: python
    :okwarning:

    aep = pw.wasp.gross_aep(wwc, wtg)

    print(aep)

Passing simply a ``wwc`` and a ``WTG`` works well when your climate is predicted for exactly the hub heights and locations of the turbines
you want to know the AEP for, but sometimes you have different turbines associated with different WTGs.
To handle this, you can pass a dictionary of WTG's and a dataset of turbine locations, heights, turbine id's, group id's, and wtg keys that
map to the WTG's in the dictionary. Let's try it here:

.. ipython:: python
    :okwarning:

    wtg_dict = {"Bonus_1_MW": wtg}
    wind_turbines = wk.create_wind_turbines_from_arrays(
        west_east=output_locs.west_east.values,
        south_north=output_locs.south_north.values,
        height=output_locs.height.values,
        wtg_keys=["Bonus_1_MW"]*15,
        group_ids=np.zeros(15, dtype=int),
        turbine_ids=np.arange(15, dtype=int),
    )

    aep = pw.wasp.gross_aep(wwc, wtg_dict, wind_turbines)

    print(aep)

Here, a single type of turbine is used, but it is possible to add more turbines to the collection by passing
several WTG's in the dictionary and adding more turbine locations to the dataset. The turbine id's and group id's
are used to group turbines together for calculating the AEP. The turbine id's are used to identify individual turbines,


Potential AEP
=============================

To add wind farm effects (such as wakes and blockage) to the AEP estimation and obtain the
"potential AEP", PyWAsP integrates with :doc:`PyWake <py_wake:index>`. PyWAsP defines a number
of named wind farm models (see :py:func:`potential_aep`), including "PARK2_onshore" and "PARK2_offshore",
that are indentical to the wake model in the WAsP GUI. Here we will use the onshore version:

.. ipython:: python
    :okwarning:

    aep_wakes = pw.wasp.potential_aep(wwc, wtg_dict, wind_turbines, wind_farm_model="PARK2_onshore", turbulence_intensity=0.0)

    print(aep_wakes)


Wind farm flow map
==================

To map an area around a wind farm PyWAsP has :py:func:`wind_farm_flow_map`, which
takes a resource grid as input and enriches it with wind farm variables including
effects of nearby turbines:

.. ipython:: python
    :okwarning:

    output_grid = wk.spatial.create_dataset(
        np.linspace(511700.0, 517700.0, 41),
        np.linspace(4619500.0, 4625500.0, 41),
        [50.0],
        crs="EPSG:32629",
        struct="cuboid",
    )

    wwc_grid = pw.wasp.downscale(gwc, topo_map, output_grid, interp_method="nearest")

    flow_map = pw.wasp.wind_farm_flow_map(
        wwc_grid,
        wtg_dict,
        wind_turbines,
        output_grid,
        wind_farm_model="PARK2_onshore",
        turbulence_intensity=0.0
    )

    print(flow_map)


We can plot the wind speed deficit for one direction to show the wake effects:

.. ipython:: python
    :okwarning:

    @savefig wspd_deficit.png align=center width=6in
    flow_map["wspd_deficit_sector"].sel(sector=210).plot(vmax=0.1)


PyWake integration in PyWAsP
============================

PyWAsP can use :doc:`PyWake <py_wake:index>` wind farm models to estimate the potential AEP. This is done by passing the name of the PyWake wind farm model to the :py:func:`potential_aep` and :py:func:`wind_farm_flow_map` functions.
To ensure stability and consistency of PyWAsP, the PyWake versions for each PyWAsP release is pinned to a specific version of PyWake. When installing PyWAsP via mamba, the correct version of PyWake is installed automatically.
We periodically update the pinned version of PyWake to the latest version that is compatible with PyWAsP.


Using custom PyWake wind farm models
------------------------------------

Instead of using a named wind farm model, it is possible to construct a
PyWake wind farm model manually, using :py:func:`functools.partial`, and use it for adding wind farm effects in PyWAsP:

.. ipython:: python
    :okwarning:

    from functools import partial
    import py_wake

    wind_farm_model = partial(
        py_wake.wind_farm_models.engineering_models.All2AllIterative,
        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(),
        ),
        superpositionModel=py_wake.superposition_models.LinearSum(),
        blockage_deficitModel=py_wake.deficit_models.Rathmann(),
        deflectionModel=None,
        turbulenceModel=None,
    )

    aep_wakes = pw.wasp.potential_aep(wwc, wtg_dict, wind_turbines, wind_farm_model=wind_farm_model, turbulence_intensity=0.0)

    print(aep_wakes)

Calculating AEP including wake effects
---------------------------------------

In PyWAsP, the "park" method is used in :py:meth:`py_wake.Site.XRSsite.from_pwc` for calculating the speed-up when passing the wind climate to a py_wake site object.

The steps used to calculate the AEP including wake effects are as follows:

1. Flow cases are defined to ensure all possible/relevant directions and speeds are covered. The wind speed in each flow case is relative to the maximum mean wind speed in the sector.
2. The wind farm simulation results from py_wake contains the variable ``WS``, which is inhomogeneous background wind speed used to obtain the probability
   density from the indivudal sector-wise weibull distributions. The variables ``WS_eff`` contains the effective wind speed including wind farm effects.
3. The AEP is determined by finding the effective power curves relative to the background wind speed ``P_eff(WS)`` and performing an integral over wind speed.