PyWAsP tutorial 4 PyWAsP - PyWake integration

Introduction

The example works through a complete wind turbine siting operation, starting with some measured wind data and ending up with a prediction of the power yield by erecting wind turbies at a specific site. This example includes AEP calculation which consider wake losses.

pywasp note: This is a classical step-by-step example from the WAsP help but extended with PyWake calculations of AEP.

Next, you will set up a wind farm consisting of several wind turbines and predict the annual energy production from this farm. Finally, we’ll map the wind resource over an area.

pywasp note: pywasp relies on the wake models from PyWake and park2 to calculate energy losses caused by the rotor wakes.

pywake note: To install pywake make sure you are in a conda enviroment containing pywasp and execute following command in a terminal pip install py_wake

pywake note: To install pywake PARK plugin make sure you are in a conda enviroment containing pywasp and pywask and execute following command in a terminal pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/cuttingedge/pywake/pywake_park.git

The situation

The company Friends of Wind Energy, E-Corp Ltd. has asked you to provide a prediction of the power yield from locating a wind turbine in Serra Santa Luzia area in Portugal. They propose to erect a single 1-MW wind turbine at the summit of Serra Santa Luzia hill (they have modest energy requirements).

No wind measurements have been taken at the turbine site itself, but data have been collected from a meteorological station at nearby hill.

You are equipped with:

  • a contour map of the area

  • the wind data from the met station

  • a simple description of the land cover in the area

  • a description of the power-generating characteristics of the turbine

These data are located in the data subfolder as following files:

  • a digital map of elevations and roughness SerraSantaLuzia.map

  • a file containing wind data SerraSantaLuzia.omwc

  • a data file containing a power production curve for the turbine Bonus_1_MW.wtg

cf303d6bd0f4449897d28547298ee371

Working with WAsP to provide a prediction

From engineering data, you know how much power will be generated by the turbine at a given wind speed. If the plan was to erect the turbine at exactly the same place where the meteorological data had been collected, then it would be a really simple task to work out how much power to expect.

However, just from looking at the map it is obvious that the proposed turbine site is completely different from the meteorological station at the airport: the properties of the meteorological station itself will affect the wind data recorded there. In addition, the properties of the turbine site will have an effect on the way that the wind behaves near the turbine. It is also unlikely that the hub height of the turbine would be the same as the height of the anemometer.

What you need is a way to take the wind climate recorded at the meteorological station, and use it to predict the wind climate at the turbine site. That is what pywasp does.

Using pywasp, you can analyse the recorded wind data, correcting for the recording site effects to produce a site-independent characterization of the local wind climate. This site-independent characterization of the local wind climate is called a wind atlas data set or generalised wind climate. You can also use pywasp to apply site effects to generalised wind climate data to produce a site-specific interpretation of the local wind climate.

Providing a prediction in the Waspdale case will therefore be a two-stage process (depicted in the right image). First, the data from the meteorological station need to be analysed to produce a wind atlas (“going up” in WAsP lingo), and then the resulting generalised wind climate needs to be applied to the proposed turbine site to estimate the wind power (“going down” in WAsP lingo).

Import packages

Usually the first step when writing a python program is importing standard and external packages we will need. For this analysis we will import numpy, pandas, matplotlib, xarray, pywasp, pywake and PARK pywake plugin.

As you work your way through the notebook, make sure to run the python code in each cell in the order that they appear. You run the code by clicking on the cell (outlines around the cell should appear) and pressing <shift> + <enter> on the keyboard.

notebook note: if something looks wrong, or when errors occur, it can be helpful to restart the python kernel via the kernel tab in the top

[1]:
import warnings
warnings.filterwarnings("ignore")  # We will ignore warnings to avoid cluttering the notebook

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

%load_ext autoreload
%autoreload 2
import pywasp as pw
import windkit as wk
import py_wake

from pywasp_tutorial import plot_bwc;

It is common to make a short alias of the package when importing using the import - as - syntax. The functions and classes of imported packages must be accessed through explicitly typing the package name or alias, e.g. np.cos(3.14) will use the cosine function from numpy.

Observed wind climate

Now we will import our (observed) binned wind climate from the mast located at hill nearby Serra Santa Luzia hill where the turbine is going to be installed. As mentioned earlier, the binned wind climates is stored in SerraSantaLuzia.omwc in the data folder. It represents the observed wind climate for the period 1997-2002 at 25.3 m above ground level.

pywasp note pywasp includes functionality to both read and write wind climate files in many data formats, including ascii (.tab), xml (.owc and .omwc), and netCDF (.nc).

The geospatial coordinate in the file is in the World Geodetic System 1984, (EPSG:4326). We will tell pywasp this by explicitly passing a keyword argument srs=4326 to the open_bwc.

[2]:
bwc = wk.read_bwc('./data/SerraSantaLuzia.omwc', crs=4326)
print(bwc)
<xarray.Dataset>
Dimensions:       (point: 1, sector: 12, wsbin: 32)
Coordinates:
    height        (point) float64 25.3
    south_north   (point) float64 41.74
    west_east     (point) float64 -8.823
    crs           int8 0
  * wsbin         (wsbin) float64 0.5 1.5 2.5 3.5 4.5 ... 28.5 29.5 30.5 31.5
    wsceil        (wsbin) float64 1.0 2.0 3.0 4.0 5.0 ... 29.0 30.0 31.0 32.0
    wsfloor       (wsbin) float64 0.0 1.0 2.0 3.0 4.0 ... 28.0 29.0 30.0 31.0
  * sector        (sector) float64 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
    sector_ceil   (sector) float64 15.0 45.0 75.0 105.0 ... 285.0 315.0 345.0
    sector_floor  (sector) float64 345.0 15.0 45.0 75.0 ... 255.0 285.0 315.0
Dimensions without coordinates: point
Data variables:
    wdfreq        (sector, point) float64 0.05314 0.03321 ... 0.1148 0.0707
    wsfreq        (wsbin, sector, point) float64 0.02601 0.04219 ... 0.0 0.0
Attributes:
    Conventions:      CF-1.8
    history:          2024-06-05T12:32:31+00:00:\twindkit==0.7.1.dev52+gc85bf...
    wasp_header:      SerraSantaluzia
    Package name:     windkit
    Package version:  0.7.1.dev52+gc85bfbd
    Creation date:    2024-06-05T12:32:32+00:00
    Object type:      Binned Wind Climate
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind Energy

Notice that the bwc object is of type <xarray.Dataset> and that it contains four kinds of data:

  1. Dimensions: core named dimensions

  2. Coordinates: coordinate values along dimensions

  3. Data variables: named arrays with data along 0-N named dimensions

  4. Attributes: additional meta data attached to the dataset

xarray note the primitive datatype and dimensions of each variable are also shown, along with a small sample of the data. wsfreq is a four-dimensional float64 (double precision) variable along dimensions (wsbin, sector, height, point)

Xarray datasets wrap numpy arrays, annotating them with human-readable dimensions and coordinates, and allowing for easy subsetting, data manipulation, and plotting of the underlying data. An xarray.Dataset object is a collection of data variables, while each varible itself has type xarray.DataArray.

xarray note Use the .values object attribute to access the underlying numpy array

Beyond the wind speed and wind direction distributions, the wind climate contains information about the height of the measurements (height) and the geospatial location (west_east and south_north), which in this case hold the location information in the projected coordinates of the EPSG:4326 projection (i.e., south_north=latitude and west_east=longitude). However, as our terrain and rougness data are in the UTM Zone 29 (EPSG:32629) we will convert and update the mast coordinates. Once the coordiantes are converted we will store them in in variables loc_x and loc_y for later use.

[3]:
bwc = wk.spatial.reproject(bwc, "EPSG:32629")
loc_y = bwc.south_north
loc_x = bwc.west_east

The next step is to plot the wind rose and wind speed distributions in the binned wind climate. For convinience a plotting function plot_bwc has been implemented in pywasp_tutorial that will do this.

notebook note: you can view the documentation for a function in jupyter notebooks by placing a ? in front of the function, and you can get the entire function by using ??.

[4]:
wk.plot.histogram_lines(bwc)

As expected the prevailing wind direction at the site is from west, thus the Atlantic ocean.

Topography data

Now that we have loaded and inspected the wind climate, the next step is to the same with the terrain and roughness data. These data are used to calculate site effects used to generalize the observed wind climate (“going up” in WAsP lingo) and to downscale to nearby locations (“going down”).

The site effects are:

  • Orographic speed-up (by sector)

  • Internal boundar layer speed-up (by sector)

  • Effective upstream surface roughness (by sector)

pywasp note WAsP includes site effects due to obstacles, but we will not consider those here

Convienietly the terrain and roughness data are provided as a single vector map SerraSantaLuzia.map. To load data we will use pw.read_vectormap().

Loading terrain and roughness data

To differentiate between the terrain and roughness parts of the vector map we need to tell pw.read_vectormap() which map type it is reading along with the EPSG code. To display loaded data we convert the vector data to raster data using method to_rastermap() and passing resolution as the argument (res = 100 m). The location of the mast is shown with an “X”.

Creating topography map

The elevation and roughness maps are the two components needed to make a TopographyMap object in pywasp, which is used to calculate the terrain effects.

In PyWAsP roughness map’s are automatically converted to a land cover maps with land cover id’s instead of roughness values. The corresponding roughness values for each id are stored in a conversion table, which is returned alongside the land cover map.

[5]:
elev_map = wk.read_vector_map("data/SerraSantaLuzia.map", crs="EPSG:32629", map_type="elevation")
pw.vectormap_to_rastermap(elev_map, res = 100).plot(figsize=(14, 10))
plt.plot(bwc.west_east, bwc.south_north, marker="X", color="white", markersize=10)

landcover_map, landcover_table = wk.read_vector_map("data/SerraSantaLuzia.map", crs="EPSG:32629", map_type="roughness")
landcover_raster = pw.vectormap_to_rastermap(landcover_map, res = 100, lctable=landcover_table)
landcover_raster.plot(cmap="tab20", figsize=(14, 10))
plt.plot(bwc.west_east, bwc.south_north, marker="X", color="black", markersize=10)

topo_map = pw.wasp.TopographyMap(elev_map, landcover_map, landcover_table)
../_images/tutorials_tutorial1_nb_10_0.png
../_images/tutorials_tutorial1_nb_10_1.png

From the maps we can clearly observe the coast of the Atlantic ocean to the west as well the river Limia which dss into the ocean. Shown in red on the roughness map.

Output Locations

Next up we will take a look at the site effects for the future turbine location. We need to tell pywasp what location we would like to calculate the site effects at. This is done by creating a xarray.Dataset with dimensions (west_east, south_north, and height) with coordinates values where we want to calculate the effects. Since this is a common need in PyWAsP, a convience function create_dataset has been created for this purpose.

We tell pywasp that the projection of west_east, south_north coordinates are in a UTM zone 29 by adding the attribute epsg giving it the integer number 32629.

[6]:
output_locs =  wk.create_dataset(
    np.array([514013.1]),
    np.array([4617763]),
    np.array([50]),
    32629
)
print(output_locs)
<xarray.Dataset>
Dimensions:      (point: 1)
Coordinates:
    height       (point) int64 50
    south_north  (point) int64 4617763
    west_east    (point) float64 5.14e+05
    crs          int8 0
Dimensions without coordinates: point
Data variables:
    output       (point) float64 0.0
Attributes:
    Conventions:  CF-1.8
    history:      2024-06-05T12:32:35+00:00:\twindkit==0.7.1.dev52+gc85bfbd\t...

Site effects

Before generalizing the observed wind climate we read in earlier, we will take a look at the site effects at the mast location.

Let us first explore the roughness elements around the mast for each of 12 sectors. We can do this by using get_rou_rose method of topo_map, and pass in the location of the mast, note that roses are 2D so the height doesn’t matter in this case, and the number of sectors we would like to use. It requires a WAsP configuration object (pywasp.wasp.Config) that stores the parameters. We will simply use the default values.

PyWAsP will use the location information on the first object passed in (it must be a xr.Dataset), so we can use the bwc object to get the effects at the mast location.

[7]:
conf = pw.wasp.Config()
rou_rose, _ = topo_map.get_rou_rose(bwc, nsecs=12, conf=conf)
wk.plot.roughness_rose(rou_rose)
[8]:
conf = pw.wasp.Config()
output_locs_mast = wk.create_dataset(loc_x, loc_y, 0, 32629)

rou_rose, _ = topo_map.get_rou_rose(output_locs_mast, 12, conf)
wk.plot.roughness_rose(rou_rose)

The presence of the Atlantic ocean is clearly visible to the west, and we see that the rose corresponds well to what we would expect from the roughness elements we saw in the map.

To get the site effects we will use the get_site_effects method from our topo_map object.

get_site_effects expects the unique west_east and south_north values, the number of sectors, and the height above ground, and the configuration object.

pywasp note depending on the number of locations, get_site_effects may take some time to complete

[9]:
site_effects = topo_map.get_site_effects(bwc, nsecs=12, conf=conf)
print(site_effects)
<xarray.Dataset>
Dimensions:              (point: 1, sector: 12)
Coordinates:
    height               (point) float64 25.3
    crs                  int8 0
    west_east            (point) float64 5.147e+05
    south_north          (point) float64 4.621e+06
  * sector               (sector) float64 0.0 30.0 60.0 ... 270.0 300.0 330.0
    sector_ceil          (sector) float64 15.0 45.0 75.0 ... 285.0 315.0 345.0
    sector_floor         (sector) float64 345.0 15.0 45.0 ... 255.0 285.0 315.0
Dimensions without coordinates: point
Data variables: (12/15)
    z0meso               (sector, point) float32 0.137 0.2776 ... 0.009736
    slfmeso              (sector, point) float32 0.9432 1.0 1.0 ... 0.403 0.5834
    displ                (sector, point) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    flow_sep_height      (sector, point) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    user_def_speedups    (sector, point) float32 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0
    orographic_speedups  (sector, point) float32 1.156 1.206 ... 1.43 1.284
    ...                   ...
    orographic_turnings  (sector, point) float32 -2.929 6.098 ... -4.751 -7.953
    obstacle_turnings    (sector, point) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    roughness_turnings   (sector, point) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    site_elev            (point) float32 381.0
    rix                  (point) float32 0.04571
    dirrix               (sector, point) float32 0.001304 0.0005211 ... 0.06163
Attributes:
    Conventions:      CF-1.8
    history:          2024-06-05T12:32:31+00:00:\twindkit==0.7.1.dev52+gc85bf...
    wasp_header:      SerraSantaluzia
    Package name:     windkit
    Package version:  0.7.1.dev52+gc85bfbd
    Creation date:    2024-06-05T12:32:37+00:00
    Object type:      Topographic effects
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind Energy
    title:            WAsP site effects

The site effects dataset is made up of following variables:

  • z0meso: effective upstream roughness by sector

  • slfmeso: effective upstream sea-land-fraction by sector

  • displ: displacement height by sector

  • user_def_speedups: user defined speed-ups (as factor)

  • orographic_speedups: speed-ups due to orography (as factor)

  • obstacle_speedups: obstacle induced speed-ups (as factor)

  • roughness_speedups: speed-ups due to roughness changes (as factor)

  • user_def_turnings: user defined wind turning (in degrees)

  • orographic_turnings: wind turning due to orography (in degrees)

  • obstacle_turnings: wind turning due to obstracels (in degrees)

  • roughness_turnings: wind turning due to roughness changes (in degrees)

  • elev: surface elevation

  • rix: ruggedness index

  • dirrix: sector wise ruggedness index

The generalized wind climate calculation

Now it is time to actually generalize our observed wind climate and then downscale it to our chosen turbine location. The generalized wind climate is a site-independent characterization of the wind climate for the entire Waspdale area. Generalization involves calculating the wind speed for predefined surface roughnesses and heights above ground. So the first thing we need to do is chose these roughness and height classes. Typically, the wind climate contains up to 5 standard heights and 5 standard roughness classes.

In a generalized wind climate data set the wind observations have been ‘cleaned’ with respect to the site-specific conditions. The data sets are site-independent and the wind distributions have been referenced to five standard roughness classes and five standard height above ground level.

[10]:
gen_roughnesses = np.array([0.0, 0.001, 0.03, 0.1, 0.4, 1.5])  # 5 standard roughness classes
gen_heights = np.array([10.0, 50.0, 100.0, 150.0, 200.0])  # 5 standard height classes

The generalization is done with the generalize function from pywasp’s wasp module. It takes the binned wind climate (bwc), our topographic map (topo_map), the roughness and height classes, and finally the configuration object (conf).

[11]:
gwc = pw.wasp.generalize(bwc, topo_map, conf, gen_roughnesses=gen_roughnesses, gen_heights=gen_heights)
print(gwc)
<xarray.Dataset>
Dimensions:        (point: 1, sector: 12, gen_roughness: 6, gen_height: 5)
Coordinates:
    height         (point) float64 25.3
    south_north    (point) float64 4.621e+06
    west_east      (point) float64 5.147e+05
    crs            int8 0
  * sector         (sector) float64 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
    sector_ceil    (sector) float64 15.0 45.0 75.0 105.0 ... 285.0 315.0 345.0
    sector_floor   (sector) float64 345.0 15.0 45.0 75.0 ... 255.0 285.0 315.0
  * gen_roughness  (gen_roughness) float64 0.0 0.001 0.03 0.1 0.4 1.5
  * gen_height     (gen_height) float64 10.0 50.0 100.0 150.0 200.0
Dimensions without coordinates: point
Data variables:
    A              (sector, gen_height, gen_roughness, point) float32 5.69 .....
    k              (sector, gen_height, gen_roughness, point) float32 1.717 ....
    wdfreq         (sector, gen_height, gen_roughness, point) float32 0.06491...
    site_elev      (point) float32 381.0
Attributes:
    Conventions:      CF-1.8
    history:          2024-06-05T12:32:31+00:00:\twindkit==0.7.1.dev52+gc85bf...
    wasp_header:      SerraSantaluzia
    Package name:     windkit
    Package version:  0.7.1.dev52+gc85bfbd
    Creation date:    2024-06-05T12:32:39+00:00
    Object type:      Geostrophic Wind Climate
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind Energy
    title:            Generalized wind climate

The generalized wind climate contains weibull parameters (A and k) for each location, sector, and height and roughness class.

Estimating wind power

Now that we have calculated a generalised wind climate (GWC) with site-independent wind climate data, we can apply those data to prospect turbine sites. PyWAsP can adjust the data for the situation found at the turbine site and produce a prediction of the wind climate for the site itself.

The pywasp.wasp module provides two ways of doing this: downscale and get_wasp_down. The latter assumes site effects that have already been calculated, while the first one will start by calculating them. Since we have already calculated the site effects, we will use get_wasp_down. The function requires passing a genwc_interp method. Two choices are available: 'nearest' and 'natural'. This is important when several generalized wind climates are present and will tell pywasp how to interpolate them to the location. Since just one wind climate is present in this study, the same generalized wind climate will be used everywhere.

To estimate the the power production for a turbine positioned at the site we need to add a turbine wind turbine generator object that holds the power curve and other information about the characteristics of the turbine. We will read a Wind Turbine Generator dataset using pw.io.read_wtg. The turbine used here is the Bonus 1 MW wind turbine.

[12]:
wtg = wk.read_wtg("./data/Bonus_1_MW.wtg")
print(wtg)
<xarray.Dataset>
Dimensions:                        (mode: 1, wind_speed: 22)
Coordinates:
  * wind_speed                     (wind_speed) float64 4.0 5.0 ... 24.0 25.0
  * mode                           (mode) int64 0
Data variables:
    power_output                   (mode, wind_speed) float64 2.41e+04 ... 1e+06
    thrust_coefficient             (mode, wind_speed) float64 0.915 ... 0.161
    air_density                    (mode) float64 1.225
    stationary_thrust_coefficient  (mode) float64 0.161
    wind_speed_cutin               (mode) float64 4.0
    wind_speed_cutout              (mode) float64 25.0
    rated_power                    (mode) float64 1e+06
    name                           <U10 'Bonus 1 MW'
    rotor_diameter                 float64 54.2
    hub_height                     float64 50.0
    regulation_type                int64 2
Attributes:
    Conventions:      CF-1.8
    Package name:     windkit
    Package version:  0.7.1.dev52+gc85bfbd
    Creation date:    2024-06-05T12:32:39+00:00
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind Energy

We can see that the dataset holds the information about the WTG, including: - Power and thrust coefficient curves for varying wind speeds - Wind speed cutin and cutout - Rated power - Rotor diameter - Hub height

The WTG can have different modes with different values for different air density situations for example. This WTG has 1 mode for the standard air density of 1.225 kg / m^3.

To plot the power and thrust coefficient curves we can use the plotting function plot_wtg_curves from PyWAsP:

[13]:
wtg["power_output"].attrs["long_name"] = "Power Output"
wtg["power_output"].attrs["units"] = "kW"
wtg["wind_speed"].attrs["long_name"] = "Wind Speed"
wtg["wind_speed"].attrs["units"] = "m/s"
wtg["thrust_coefficient"].attrs["long_name"] = "Ct"
wtg["thrust_coefficient"].attrs["units"] = ""
wk.plot.power_ct_curves(wtg)

Now that we have the turbine data we can combine it with a location in our map to predict the power. We will add the location and hub height of our turbine to xr.Dataset we can create by using the built-in create_dataset in PyWAsP:

[14]:
x, y, hub_height = 514325, 4622375, 50.0
turbine_loc =  wk.create_dataset([x], [y], [hub_height], crs="EPSG:32629")

Lets plot the turbine location on the elevation map to see where it’s located.

[15]:
fig, ax = plt.subplots(1, 1, figsize=(12, 8.5))
pw.vectormap_to_rastermap(elev_map, res = 100).plot(ax=ax, vmin=0, vmax=550)
ax.plot(turbine_loc.west_east, turbine_loc.south_north, marker="X", markersize=8, color="black")
ax.plot(bwc.west_east, bwc.south_north, marker="X", color="white", markersize=8)
ax.set(xlabel="west_east [m]", ylabel="south_north [m]");
../_images/tutorials_tutorial1_nb_31_0.png

We can see that the wind turbine is located on a local peak on the hill to the west in the map near the coastline. Next up, we will calculate the power production for the turbine. To do this, we will first calculate the wind climate at the location of the turbine hub height using PyWAsP’s downscale function which takes in a generalized wind climate, the topography data, the turbine position(s) and hub height(s), the wasp configuration, and we also add the genwc_interp=nearest" keyword to tell PyWAsP we want it to use the nearest GWC to each turbine (in this case its always the same, since we only have one GWC).

[16]:
pwc = pw.wasp.downscale(gwc, topo_map, turbine_loc, conf, interp_method="nearest")
print(pwc)
<xarray.Dataset>
Dimensions:        (sector: 12, point: 1)
Coordinates:
  * sector         (sector) float64 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
    height         (point) float64 50.0
    south_north    (point) int64 4622375
    west_east      (point) int64 514325
    crs            int8 0
    sector_ceil    (sector) float64 15.0 45.0 75.0 105.0 ... 285.0 315.0 345.0
    sector_floor   (sector) float64 345.0 15.0 45.0 75.0 ... 255.0 285.0 315.0
Dimensions without coordinates: point
Data variables:
    A              (sector, point) float32 7.04 6.572 6.188 ... 7.326 7.172
    k              (sector, point) float32 2.041 2.385 2.498 ... 2.248 2.092
    wdfreq         (sector, point) float32 0.06301 0.0392 ... 0.09475 0.08239
    site_elev      (point) float32 460.0
    air_density    (point) float32 1.163
    wspd           (point) float32 6.998
    power_density  (point) float32 344.5
Attributes:
    Conventions:      CF-1.8
    history:          2024-06-05T12:32:39+00:00:\twindkit==0.7.1.dev52+gc85bf...
    title:            WAsP site effects
    Package name:     windkit
    Package version:  0.7.1.dev52+gc85bfbd
    Creation date:    2024-06-05T12:32:42+00:00
    Object type:      Met fields
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind Energy

We can now combine the predicted wind climate with the wind turbine generator to get the annual power production.

[17]:
aep = pw.wasp.gross_aep(pwc, wtg)
print(aep)
<xarray.Dataset>
Dimensions:      (point: 1)
Coordinates:
    height       (point) float64 50.0
    south_north  (point) int64 4622375
    west_east    (point) int64 514325
    crs          int8 0
Dimensions without coordinates: point
Data variables:
    gross_AEP    (point) float32 2.568
Attributes:
    Conventions:      CF-1.8
    history:          2024-06-05T12:32:39+00:00:\twindkit==0.7.1.dev52+gc85bf...
    title:            WAsP site effects
    Package name:     windkit
    Package version:  0.7.1.dev52+gc85bfbd
    Creation date:    2024-06-05T12:32:42+00:00
    Object type:      Anual Energy Production
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind Energy

As we can see, the predicted power production for the turbine at the site is 2.6 GWh/y

Estimating wind farm production

Now that we have estimated the power for a single turbine, we can look at the power production for a whole wind farm. The workflow is indentical to the one we just used for a single turbine. We have predefined some turbine positions in the data/turbine_positions.csv file, so lets create the location dataset and plot the turbines on our map.

[18]:
wtg_locs = pd.read_csv('./data/turbine_positions.csv')

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

fig, ax = plt.subplots(1, 1, figsize=(12, 8.5))
pw.vectormap_to_rastermap(elev_map, res = 100).plot(ax=ax, vmin=0, vmax=550)
ax.plot(output_locs.west_east, output_locs.south_north, ls="none", marker="X", markersize=8, color="black")
ax.plot(bwc.west_east, bwc.south_north, marker="X", color="white", markersize=8)

for i, (x, y) in enumerate(zip(output_locs.west_east.values, output_locs.south_north.values)):
    ax.annotate(i, (x, y+50))
ax.set(xlabel="west_east [m]", ylabel="south_north [m]", xlim=(510000, 519000), ylim=(4620700, 4625500))
[18]:
[Text(0.5, 0, 'west_east [m]'),
 Text(0, 0.5, 'south_north [m]'),
 (510000.0, 519000.0),
 (4620700.0, 4625500.0)]
../_images/tutorials_tutorial1_nb_38_1.png

We can see that the turbines are mostly on the (locally) highest points in the terrain.

To calculate the predicted wind climates at the 15 positions (at 50 m AGL) we will again use the downscale function.

[19]:
pwc = pw.wasp.downscale(gwc, topo_map, output_locs, conf, interp_method="nearest")

Now we can calculate the AEP

[20]:
aep_nowake = pw.wasp.gross_aep(pwc, wtg)["gross_AEP"]
print(aep_nowake)
print(aep_nowake.sum())
<xarray.DataArray 'gross_AEP' (point: 15)>
array([2.9344273, 2.794619 , 2.5766804, 2.6569452, 2.8252776, 2.9034677,
       2.7838368, 3.0851603, 3.18083  , 2.4066768, 2.8184903, 3.058005 ,
       3.0135605, 3.2057064, 3.4222102], dtype=float32)
Coordinates:
    height       (point) int64 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
    south_north  (point) int64 4622313 4622199 4622336 ... 4624252 4624142
    west_east    (point) int64 513914 514161 514425 ... 515808 516060 516295
    crs          int8 0
Dimensions without coordinates: point
Attributes:
    description:    Gross annual energy production, without considering wake ...
    long_name:      Gross annual energy production
    standard_name:  gross_aep
    units:          GWh
    grid_mapping:   crs
<xarray.DataArray 'gross_AEP' ()>
array(43.665897, dtype=float32)
Coordinates:
    crs      int8 0

We get back an AEP for each turbine. In total we estimate that our wind farm produces 43.7 GWh annually, without taking wake effects and other losses into account.

Lets now use wake effects in the AEP calculations to see how much these affect our annual power production. PyWAsP uses PyWake internally to take wake effects into account.

Several wake models can be used, but we will use the wake model that corresponds to the wake model found in the offical WAsP GUI software for Windows, called “PARK2” in PyWAsP.

[21]:
aep_park2 = pw.wasp.potential_aep(
    pwc[["A", "k", "wdfreq"]],
    wtg,
    wind_farm_model="PARK2_onshore",
    turbulence_intensity=0.0
)["potential_AEP_sector"]

# Sum over sectors
aep_park2 = aep_park2.sum(dim="sector")

print(aep_park2)
print(aep_park2.sum())
<xarray.DataArray 'potential_AEP_sector' (point: 15)>
array([2.83374745, 2.70305958, 2.46916798, 2.54701922, 2.71523867,
       2.86011705, 2.75593011, 3.03593558, 3.12073931, 2.33494791,
       2.76085495, 2.97103928, 2.86495357, 3.03652726, 3.3034571 ])
Coordinates:
    height       (point) int64 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
    south_north  (point) int64 4622313 4622199 4622336 ... 4624252 4624142
    west_east    (point) int64 513914 514161 514425 ... 515808 516060 516295
    crs          int8 0
  * point        (point) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
<xarray.DataArray 'potential_AEP_sector' ()>
array(42.31273502)
Coordinates:
    crs      int8 0

We can see that wake effects causes a reduction of the AEP from 43.7 to 42.3 GWh/y (a 3% reduction).

Let us visaulize the turbine-by-turbine AEP and reduction by wake effects.

[22]:
fig, axes = plt.subplots(1, 2, figsize=(15, 5))
aep_nowake.plot(ax=axes[0], marker="o", label=f"No wakes: {aep_nowake.sum().values:.1f} [GWh]")
aep_park2.plot(ax=axes[0], marker="o", label=f"PARK2 wakes: {aep_park2.sum().values:.1f} [GWh]")

(aep_park2/ aep_nowake).plot(ax=axes[1], marker="o", color="C1")

axes[0].legend(frameon=False, title="Wake model")
axes[0].set(ylabel="AEP [GWh]", ylim=(2.2, 4.4), title="")
axes[1].set(ylabel=r"AEP / AEP$_\mathrm{no wakes}$", title="");
../_images/tutorials_tutorial1_nb_46_0.png

The First plot shows the AEP for each turbine. It varies alot from turbine to turbine due to the differences in the local wind climate.

The second plot shows the relative decrease in AEP due to wake effects. We can see that turbine at index 12 and 13 (Index starts at 0) are affected (relatively) most by wake effects, with their AEP reduced by about 5%.

[23]:
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, conf, interp_method="natural")
pwc_grid
[23]:
<xarray.Dataset>
Dimensions:        (height: 1, south_north: 41, west_east: 41, sector: 12)
Coordinates:
  * height         (height) float64 50.0
  * south_north    (south_north) float64 4.62e+06 4.62e+06 ... 4.626e+06
  * west_east      (west_east) float64 5.117e+05 5.118e+05 ... 5.177e+05
  * sector         (sector) float64 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
    crs            int8 0
    sector_ceil    (sector) float64 15.0 45.0 75.0 105.0 ... 285.0 315.0 345.0
    sector_floor   (sector) float64 345.0 15.0 45.0 75.0 ... 255.0 285.0 315.0
Data variables:
    A              (sector, height, south_north, west_east) float32 5.162 ......
    k              (sector, height, south_north, west_east) float32 1.982 ......
    wdfreq         (sector, height, south_north, west_east) float32 0.06678 ....
    site_elev      (south_north, west_east) float32 15.3 22.06 ... 215.2 199.8
    air_density    (height, south_north, west_east) float32 1.216 ... 1.194
    wspd           (height, south_north, west_east) float32 5.545 5.387 ... 4.86
    power_density  (height, south_north, west_east) float32 200.1 ... 121.8
Attributes:
    Conventions:      CF-1.8
    history:          2024-06-05T12:32:58+00:00:\twindkit==0.7.1.dev52+gc85bf...
    Package name:     windkit
    Package version:  0.7.1.dev52+gc85bfbd
    Creation date:    2024-06-05T12:33:05+00:00
    Object type:      Met fields
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind Energy
    title:            WAsP site effects
[24]:
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,
    turbine_ids=np.arange(len(output_locs.west_east.values), dtype=int),
    group_ids=np.zeros_like(output_locs.west_east.values, dtype=int),
    wtg_keys=["bonus_1_mw"]*len(output_locs.west_east.values),
    crs="EPSG:32629",
)

flow_map = pw.wasp.wind_farm_flow_map(
    pwc_grid,
    wtg_dict,
    wind_turbines,
    output_grid,
    wind_farm_model="PARK2_onshore",
    turbulence_intensity=0.0,
    n_subsector=5,
)
print(flow_map)
<xarray.Dataset>
Dimensions:                          (sector: 12, height: 1, south_north: 41,
                                      west_east: 41, mode: 1)
Coordinates:
  * sector                           (sector) float64 0.0 30.0 ... 300.0 330.0
  * height                           (height) float64 50.0
  * south_north                      (south_north) float64 4.62e+06 ... 4.626...
  * west_east                        (west_east) float64 5.117e+05 ... 5.177e+05
    crs                              int8 0
  * mode                             (mode) int64 0
Data variables: (12/15)
    potential_AEP_sector             (height, south_north, west_east, sector) float64 ...
    gross_AEP_sector                 (height, south_north, west_east, sector, mode) float64 ...
    AEP_deficit_sector               (height, south_north, west_east, sector, mode) float64 ...
    wspd_sector                      (height, south_north, west_east, sector) float64 ...
    wspd_eff_sector                  (height, south_north, west_east, sector) float64 ...
    wspd_deficit_sector              (height, south_north, west_east, sector) float64 ...
    ...                               ...
    gross_AEP                        (height, south_north, west_east, mode) float64 ...
    AEP_deficit                      (height, south_north, west_east, mode) float64 ...
    wspd                             (height, south_north, west_east) float64 ...
    wspd_eff                         (height, south_north, west_east) float64 ...
    wspd_deficit                     (height, south_north, west_east) float64 ...
    turbulence_intensity_eff         (height, south_north, west_east) float64 ...
Attributes:
    Conventions:      CF-1.8
    history:          2024-06-05T12:32:58+00:00:\twindkit==0.7.1.dev52+gc85bf...
    Package name:     windkit
    Package version:  0.7.1.dev52+gc85bfbd
    Creation date:    2024-06-05T12:33:19+00:00
    Object type:      Wind Farm Flow Map
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind Energy
    title:            WAsP site effects
[25]:
fig, axes = plt.subplots(2, 3, figsize=(15, 7.5))
flow_map["wspd_sector"].sel(sector=240).plot(ax=axes[0, 0], levels=np.linspace(2.0, 10.0, 17))
flow_map["wspd_eff_sector"].sel(sector=240).plot(ax=axes[0, 1], levels=np.linspace(2.0, 10.0, 17))
(100.0 * flow_map["wspd_deficit_sector"]).sel(sector=240).plot(ax=axes[0, 2], levels=np.linspace(-10.0, 10.0, 21))
flow_map["wspd"].plot(ax=axes[1, 0], levels=np.linspace(2.0, 10.0, 17))
flow_map["wspd_eff"].plot(ax=axes[1, 1], levels=np.linspace(2.0, 10.0, 17))
(100.0 * flow_map["wspd_deficit"]).plot(ax=axes[1, 2], levels=np.linspace(-10.0, 10.0, 21))
for ax in axes.flat:
    ax.plot(output_locs.west_east, output_locs.south_north, ls="none", marker="X", markersize=8, color="black")
    ax.plot(bwc.west_east, bwc.south_north, marker="X", color="white", markersize=8)
    ax.set(title="")
plt.tight_layout()
../_images/tutorials_tutorial1_nb_50_0.png
[26]:
fig, axes = plt.subplots(2, 3, figsize=(15, 7.5))
flow_map["gross_AEP_sector"].sel(sector=240).plot(ax=axes[0, 0])
flow_map["potential_AEP_sector"].sel(sector=240).plot(ax=axes[0, 1])
(100.0 * flow_map["AEP_deficit_sector"]).sel(sector=240).plot(ax=axes[0, 2], levels=np.linspace(-20.0, 20.0, 21))
flow_map["gross_AEP"].plot(ax=axes[1, 0])
flow_map["potential_AEP"].plot(ax=axes[1, 1])
(100.0 * flow_map["AEP_deficit"]).plot(ax=axes[1, 2], levels=np.linspace(-20.0, 20.0, 21))
for ax in axes.flat:
    ax.plot(output_locs.west_east, output_locs.south_north, ls="none", marker="X", markersize=8, color="black")
    ax.plot(bwc.west_east, bwc.south_north, marker="X", color="white", markersize=8)
    ax.set(title="")
plt.tight_layout()
../_images/tutorials_tutorial1_nb_51_0.png
[ ]: