.. 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

    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

    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

    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

    pwc = pw.wasp.generalize_and_downscale(output_locs, bwc, topo_map)
    print(pwc)


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

    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

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

    print(aep)

Instead of a WTG, a more flexible option is to use the :py:class:`windkit.WindTurbines` object instead. It holds
a collection of WTGs and associated wind turbine locations. This is useful when estimating the AEP of a wind farm
with multiple turbine types:

.. ipython:: python

    wind_turbines = wk.WindTurbines([(wtg, output_locs)])

    aep = pw.wasp.gross_aep(pwc, 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 tuples of WTGs and locations to the :py:class:`windkit.WindTurbines` constructor.


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

    aep_wakes = pw.wasp.potential_aep(pwc, wind_turbines, wind_farm_model="PARK2_onshore")

    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

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

    pwc_grid = pw.wasp.downscale(gwc, topo_map, output_grid, genwc_interp="nearest")

    flow_map = pw.wasp.wind_farm_flow_map(
        pwc_grid,
        wind_turbines,
        output_grid,
        wind_farm_model="PARK2_onshore",
    )

    print(flow_map)


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

.. ipython:: python

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



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

    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(pwc, wtg, wind_farm_model=wind_farm_model)

    print(aep_wakes)