Classical resource assessment step-by-step

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. 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 is installed with PyWAsP, so you should be ready to go. If not, you can install PyWake from a conda enviroment containing pywasp by 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

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.

• numpy is python’s main numerical array package https://www.numpy.org/

• pandas is python’s main package for manipulation of tabular data https://pandas.pydata.org/

• matplotlib is python’s main plotting package https://matplotlib.org/

• xarray is a powerful high level package for labelled multi-dimensional arrays http://xarray.pydata.org

• pywasp docummentation is found at http://docs.wasp.dk/

• py_wake docummentation is found at https://topfarm.pages.windenergy.dtu.dk/PyWake/

• py_wake_park docummentation is found at https://gitlab.windenergy.dtu.dk/TOPFARM/cuttingedge/pywake/pywake_park

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

%load_ext autoreload
%autoreload 2
import warnings

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

import windkit as wk
import pywasp as pw

warnings.filterwarnings("ignore")  # We will ignore warnings to avoid cluttering the notebook

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), better known as a simple latitude/longitude. When possible, the file readers will try to parse the coordinate reference system (CRS). We can tell pywasp this by explicitly passing the keyword argument crs=4326 to the open_bwc.

bwc = wk.read_bwc('./data/SerraSantaLuzia.omwc')
print(bwc)

<xarray.Dataset> Size: 4kB
Dimensions:       (point: 1, sector: 12, wsbin: 32)
Coordinates:
    height        (point) float64 8B 25.3
    south_north   (point) float64 8B 41.74
    west_east     (point) float64 8B -8.823
    crs           int8 1B 0
    wsceil        (wsbin) float64 256B 1.0 2.0 3.0 4.0 ... 29.0 30.0 31.0 32.0
    wsfloor       (wsbin) float64 256B 0.0 1.0 2.0 3.0 ... 28.0 29.0 30.0 31.0
    sector_ceil   (sector) float64 96B 15.0 45.0 75.0 ... 285.0 315.0 345.0
    sector_floor  (sector) float64 96B 345.0 15.0 45.0 ... 255.0 285.0 315.0
  * wsbin         (wsbin) float64 256B 0.5 1.5 2.5 3.5 ... 28.5 29.5 30.5 31.5
  * sector        (sector) float64 96B 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
Dimensions without coordinates: point
Data variables:
    wdfreq        (sector, point) float64 96B 0.05314 0.03321 ... 0.1148 0.0707
    wsfreq        (wsbin, sector, point) float64 3kB 0.02601 0.04219 ... 0.0 0.0
Attributes:
    Conventions:      CF-1.8
    history:          2025-07-04T13:07:08+00:00:\twindkit==1.0.0\tcreate_data...
    description:      SerraSantaluzia
    Package name:     windkit
    Package version:  1.0.0
    Creation date:    2025-07-04T13:07:08+00:00
    Object type:      Binned Wind Climate
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind

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 could convert and update the mast coordinates explicitly by using bwc = wk.spatial.reproject(bwc, "EPSG:32629"). This is not strictly needed usually, because the bwc will be automatically reprojected to the coordinate system of the map. However, here the coordiantes are converted because we will store them in in variables loc_x and loc_y for later use.

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.

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

wk.plot.histogram_lines(bwc)

Data type cannot be displayed: application/vnd.plotly.v1+json

As expected the prevailing wind direction at the site is from west (from 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

Conveniently the terrain and roughness data are provided as a single vector map SerraSantaLuzia.map. To load data we will use wk.read_roughness_map() whereas for the elevation data we can use wk.read_elevation_map().

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 and roughness (z0) and displacement (d) values.

elev_map = wk.read_elevation_map("data/SerraSantaLuzia.map")

landcover_map = wk.read_roughness_map(
    "data/SerraSantaLuzia.map",
    crs="EPSG:32629"
)

topo_map = pw.wasp.TopographyMap(elev_map, landcover_map)

Landcover maps in WAsP are characterized by polygons, that denotes the properties of a certain type of landcover. The properties used are roughness (z0), and optionally a displacement height (d) and description (desc).

landcover_map

	id	geometry	z0	d	desc
0	0.0	POLYGON ((520129.5 4608000, 520157.1 4608018.5...	0.03	0.0
1	0.0	POLYGON ((520157.1 4608018.5, 520129.5 4608000...	0.03	0.0
2	0.0	POLYGON ((522602.4 4608000, 522613.5 4608058, ...	0.03	0.0
3	1.0	POLYGON ((522613.5 4608058, 522602.4 4608000, ...	0.60	0.0
4	1.0	POLYGON ((518024.6 4608000, 518062.3 4608074.5...	0.60	0.0
...	...	...	...	...	...
338	5.0	POLYGON ((525157.1 4637000, 525155.9 4636996.5...	0.10	0.0
339	1.0	POLYGON ((526041.6 4637000, 526024.6 4636974, ...	0.60	0.0
340	0.0	POLYGON ((512276.3 4637000, 512271.8 4636986, ...	0.03	0.0
341	1.0	POLYGON ((524201 4636832.5, 524152.8 4636863, ...	0.60	0.0
342	4.0	POLYGON ((509780.5 4637000, 509768.8 4636924, ...	0.00	0.0

343 rows × 5 columns

Displaying terrain and roughness data

To display the elevation roughness we can call wk.plot.elevation_map and wk.plot.landcover_map. These are plotting functions form windkit that use some default color schemes to visualize the elevation and roughness. The roughness is colored by typical topographical map land cover types (red=city, light green=agricultural, dark green=forest). We can add the location of the mast by converting the bwc object (a xarray dataset, ds) to a geopandas GeoDataFrame (gdf) by using the function ds_to_gdf. Geopandas is the engine used for processing vector data in windkit/pywasp. By adding this GeoDataFrame to the plotting axes, the mast location is shown with an “X”.

ax = wk.plot.elevation_map(elev_map, zorder=0)
wk.spatial.ds_to_gdf(bwc).plot(ax=ax, zorder=1, color="red", marker="X", markersize=40)

<Axes: title={'center': 'Elevation contour map'}, xlabel='Easting (m)', ylabel='Northing (m)'>

ax_lc = wk.plot.landcover_map(landcover_map, zorder=0)
wk.spatial.ds_to_gdf(bwc).plot(ax=ax_lc, zorder=1, color="black", marker="X", markersize=40)

<Axes: xlabel='Easting (m)', ylabel='Northing (m)'>

From the maps we can clearly observe the coast of the Atlantic ocean to the west as well as the river Limia which flows into the ocean. The geopandas dataframe also has a nice function called explore(), that opens the map in interactive window where you can zoom in on the features of the map and hover over areas.

landcover_map.explore("z0")

Make this Notebook Trusted to load map: File -> Trust Notebook

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.

output_locs =  wk.spatial.create_point(
    west_east=np.array([514013.1]),
    south_north=np.array([4617763]),
    height=np.array([50]),
    crs=32629
)
print(output_locs)

<xarray.Dataset> Size: 33B
Dimensions:      (point: 1)
Coordinates:
    height       (point) int64 8B 50
    south_north  (point) int64 8B 4617763
    west_east    (point) float64 8B 5.14e+05
    crs          int8 1B 0
Dimensions without coordinates: point
Data variables:
    output       (point) float64 8B 0.0
Attributes:
    Conventions:  CF-1.8
    history:      2025-07-04T13:07:15+00:00:\twindkit==1.0.0\t"point")

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.

rou_rose, _ = topo_map.get_rou_rose(bwc, n_sectors=12)
wk.plot.roughness_rose(rou_rose)

Data type cannot be displayed: application/vnd.plotly.v1+json

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

site_effects = topo_map.get_site_effects(bwc, n_sectors=12)
print(site_effects)
print(site_effects.data_vars)

<xarray.Dataset> Size: 945B
Dimensions:              (point: 1, sector: 12)
Coordinates:
    height               (point) float64 8B 25.3
    crs                  int8 1B 0
    west_east            (point) float64 8B 5.147e+05
    south_north          (point) float64 8B 4.621e+06
    sector_ceil          (sector) float64 96B 15.0 45.0 75.0 ... 315.0 345.0
    sector_floor         (sector) float64 96B 345.0 15.0 45.0 ... 285.0 315.0
  * sector               (sector) float64 96B 0.0 30.0 60.0 ... 300.0 330.0
Dimensions without coordinates: point
Data variables: (12/15)
    z0meso               (sector, point) float32 48B 0.137 0.2776 ... 0.009736
    slfmeso              (sector, point) float32 48B 0.9432 1.0 ... 0.403 0.5834
    displ                (sector, point) float32 48B 0.0 0.0 0.0 ... 0.0 0.0 0.0
    flow_sep_height      (sector, point) float32 48B 0.0 0.0 0.0 ... 0.0 0.0 0.0
    user_def_speedups    (sector, point) float32 48B 1.0 1.0 1.0 ... 1.0 1.0 1.0
    orographic_speedups  (sector, point) float32 48B 1.156 1.206 ... 1.43 1.284
    ...                   ...
    orographic_turnings  (sector, point) float32 48B -2.929 6.098 ... -7.953
    obstacle_turnings    (sector, point) float32 48B 0.0 0.0 0.0 ... 0.0 0.0 0.0
    roughness_turnings   (sector, point) float32 48B 0.0 0.0 0.0 ... 0.0 0.0 0.0
    site_elev            (point) float32 4B 381.0
    rix                  (point) float32 4B 0.04571
    dirrix               (sector, point) float32 48B 0.001304 ... 0.06163
Attributes:
    Conventions:      CF-1.8
    history:          2025-07-04T13:07:08+00:00:\twindkit==1.0.0\tcreate_data...
    description:      SerraSantaluzia
    Package name:     windkit
    Package version:  1.0.0
    Creation date:    2025-07-04T13:07:17+00:00
    Object type:      Topographic effects
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind
    title:            WAsP site effects
Data variables:
    z0meso               (sector, point) float32 48B 0.137 0.2776 ... 0.009736
    slfmeso              (sector, point) float32 48B 0.9432 1.0 ... 0.403 0.5834
    displ                (sector, point) float32 48B 0.0 0.0 0.0 ... 0.0 0.0 0.0
    flow_sep_height      (sector, point) float32 48B 0.0 0.0 0.0 ... 0.0 0.0 0.0
    user_def_speedups    (sector, point) float32 48B 1.0 1.0 1.0 ... 1.0 1.0 1.0
    orographic_speedups  (sector, point) float32 48B 1.156 1.206 ... 1.43 1.284
    obstacle_speedups    (sector, point) float32 48B 1.0 1.0 1.0 ... 1.0 1.0 1.0
    roughness_speedups   (sector, point) float32 48B 1.008 1.026 ... 0.891
    user_def_turnings    (sector, point) float32 48B 0.0 0.0 0.0 ... 0.0 0.0 0.0
    orographic_turnings  (sector, point) float32 48B -2.929 6.098 ... -7.953
    obstacle_turnings    (sector, point) float32 48B 0.0 0.0 0.0 ... 0.0 0.0 0.0
    roughness_turnings   (sector, point) float32 48B 0.0 0.0 0.0 ... 0.0 0.0 0.0
    site_elev            (point) float32 4B 381.0
    rix                  (point) float32 4B 0.04571
    dirrix               (sector, point) float32 48B 0.001304 ... 0.06163

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.

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

gwc = pw.wasp.generalize(bwc, topo_map, gen_roughnesses=gen_roughnesses, gen_heights=gen_heights)
print(gwc)

<xarray.Dataset> Size: 5kB
Dimensions:        (point: 1, sector: 12, gen_height: 5, gen_roughness: 6)
Coordinates:
    height         (point) float64 8B 25.3
    south_north    (point) float64 8B 4.621e+06
    west_east      (point) float64 8B 5.147e+05
    crs            int8 1B 0
    sector_ceil    (sector) float64 96B 15.0 45.0 75.0 ... 285.0 315.0 345.0
    sector_floor   (sector) float64 96B 345.0 15.0 45.0 ... 255.0 285.0 315.0
  * sector         (sector) float64 96B 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
  * gen_roughness  (gen_roughness) float64 48B 0.0 0.001 0.03 0.1 0.4 1.5
  * gen_height     (gen_height) float64 40B 10.0 50.0 100.0 150.0 200.0
Dimensions without coordinates: point
Data variables:
    A              (sector, gen_height, gen_roughness, point) float32 1kB 5.6...
    k              (sector, gen_height, gen_roughness, point) float32 1kB 1.7...
    wdfreq         (sector, gen_height, gen_roughness, point) float32 1kB 0.0...
    site_elev      (point) float32 4B 381.0
Attributes:
    Conventions:      CF-1.8
    history:          2025-07-04T13:07:08+00:00:\twindkit==1.0.0\tcreate_data...
    description:      SerraSantaluzia
    Package name:     windkit
    Package version:  1.0.0
    Creation date:    2025-07-04T13:07:18+00:00
    Object type:      Geostrophic Wind Climate
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind
    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.

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

<xarray.Dataset> Size: 704B
Dimensions:                        (mode: 1, wind_speed: 22)
Coordinates:
  * wind_speed                     (wind_speed) float64 176B 4.0 5.0 ... 25.0
  * mode                           (mode) int64 8B 0
Data variables:
    power_output                   (mode, wind_speed) float64 176B 2.41e+04 ....
    thrust_coefficient             (mode, wind_speed) float64 176B 0.915 ... ...
    air_density                    (mode) float64 8B 1.225
    stationary_thrust_coefficient  (mode) float64 8B 0.161
    wind_speed_cutin               (mode) float64 8B 4.0
    wind_speed_cutout              (mode) float64 8B 25.0
    rated_power                    (mode) float64 8B 1e+06
    name                           <U10 40B 'Bonus 1 MW'
    manufacturer                   <U16 64B 'Bonus Energy A/S'
    rotor_diameter                 float64 8B 54.2
    hub_height                     float64 8B 50.0
    regulation_type                int64 8B 2
Attributes:
    Conventions:      CF-1.8
    Package name:     windkit
    Package version:  1.0.0
    Creation date:    2025-07-04T13:07:19+00:00
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind

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.

Nonetheless, the WTG power curve will, by default, get corrected for the site-specific air density when estiamting the AEP. A higher than 1.225 kg / m^3 air density will result in an upwards shift of the power curve, while a lower air density will result in a downwards shift.

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

wk.plot.power_ct_curves(wtg)

Data type cannot be displayed: application/vnd.plotly.v1+json

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:

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

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

ax = wk.plot.elevation_map(elev_map, zorder=0)
wk.spatial.ds_to_gdf(bwc).plot(ax=ax, color="red", marker="X", markersize=40, zorder=1)
wk.spatial.ds_to_gdf(turbine_loc).plot(ax=ax, color="black", marker="X", markersize=40, zorder=2)

<Axes: title={'center': 'Elevation contour map'}, xlabel='Easting (m)', ylabel='Northing (m)'>

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

pwc = pw.wasp.downscale(gwc, topo_map, turbine_loc, interp_method="nearest")
print(pwc)

<xarray.Dataset> Size: 473B
Dimensions:        (point: 1, sector: 12)
Coordinates:
    height         (point) float64 8B 50.0
    south_north    (point) int64 8B 4622375
    west_east      (point) int64 8B 514325
    crs            int8 1B 0
    sector_ceil    (sector) float64 96B 15.0 45.0 75.0 ... 285.0 315.0 345.0
    sector_floor   (sector) float64 96B 345.0 15.0 45.0 ... 255.0 285.0 315.0
  * sector         (sector) float64 96B 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
Dimensions without coordinates: point
Data variables:
    A              (sector, point) float32 48B 7.041 6.572 6.188 ... 7.326 7.172
    k              (sector, point) float32 48B 2.041 2.385 2.498 ... 2.248 2.092
    wdfreq         (sector, point) float32 48B 0.063 0.0392 ... 0.09475 0.08239
    site_elev      (point) float32 4B 460.0
    air_density    (point) float32 4B 1.163
    wspd           (point) float32 4B 6.998
    power_density  (point) float32 4B 344.5
Attributes:
    Conventions:      CF-1.8
    history:          2025-07-04T13:07:19+00:00:\twindkit==1.0.0\t"point")\n2...
    title:            WAsP site effects
    Package name:     windkit
    Package version:  1.0.0
    Creation date:    2025-07-04T13:07:21+00:00
    Object type:      Met fields
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind

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

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

<xarray.Dataset> Size: 365B
Dimensions:           (point: 1, sector: 12)
Coordinates:
    height            (point) float64 8B 50.0
    south_north       (point) int64 8B 4622375
    west_east         (point) int64 8B 514325
    crs               int8 1B 0
    sector_ceil       (sector) float64 96B 15.0 45.0 75.0 ... 285.0 315.0 345.0
    sector_floor      (sector) float64 96B 345.0 15.0 45.0 ... 255.0 285.0 315.0
  * sector            (sector) float64 96B 0.0 30.0 60.0 ... 270.0 300.0 330.0
Dimensions without coordinates: point
Data variables:
    gross_aep         (point) float32 4B 2.466
    gross_aep_sector  (point, sector) float32 48B 0.1235 0.06057 ... 0.1672
Attributes:
    Conventions:      CF-1.8
    history:          2025-07-04T13:07:19+00:00:\twindkit==1.0.0\t"point")\n2...
    title:            WAsP site effects
    Package name:     windkit
    Package version:  1.0.0
    Creation date:    2025-07-04T13:07:21+00:00
    Object type:      Anual Energy Production
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind

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.

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

wind_turbines = wk.create_wind_turbines_from_arrays(
    west_east=wtg_locs["Easting"].values,
    south_north=wtg_locs["Northing"].values,
    height=wtg_locs["Hub height"].values,
    wtg_keys=["Bonus_1_MW"]*wtg_locs.shape[0],
    crs="EPSG:32629"
)

wind_turbines

<xarray.Dataset> Size: 1kB
Dimensions:      (point: 15)
Coordinates:
    height       (point) int64 120B 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
    south_north  (point) int64 120B 4622313 4622199 4622336 ... 4624252 4624142
    west_east    (point) int64 120B 513914 514161 514425 ... 516060 516295
    crs          int8 1B 0
    turbine_id   (point) int64 120B 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
    group_id     (point) float64 120B 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    wtg_key      (point) <U10 600B 'Bonus_1_MW' 'Bonus_1_MW' ... 'Bonus_1_MW'
Dimensions without coordinates: point
Data variables:
    *empty*
Attributes:
    Conventions:  CF-1.8
    history:      2025-07-04T13:07:21+00:00:\twindkit==1.0.0\tcreate_dataset(

xarray.Dataset

• Dimensions:
  • point: 15
• Coordinates: (7)
  • height
    (point)
    int64
    50 50 50 50 50 ... 50 50 50 50 50

    URI :

    axis :

    Z

    description :

    It is the height above ground either related to the measurements or to the hub of a turbine.

    long_name :

    Height above ground

    standard_name :

    height

    units :

    m

    array([50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50])

  • south_north
    (point)
    int64
    4622313 4622199 ... 4624252 4624142

    URI :

    axis :

    Y

    description :

    Positive northward direction in projected coordinates (meters).

    long_name :

    Projected south to north coordinate

    standard_name :

    projection_y_coordinate

    units :

    m

    array([4622313, 4622199, 4622336, 4622733, 4622683, 4622617, 4621743,
           4622164, 4622558, 4622605, 4623662, 4624074, 4624146, 4624252,
           4624142])

  • west_east
    (point)
    int64
    513914 514161 ... 516060 516295

    URI :

    axis :

    X

    description :

    Positive eastward direction in projected coordinates (meters).

    long_name :

    Projected west to east coordinate

    standard_name :

    projection_x_coordinate

    units :

    m

    array([513914, 514161, 514425, 513537, 513295, 513058, 515340, 515892,
           515988, 514808, 515711, 515550, 515808, 516060, 516295])

  • crs
    ()
    int8
    0

    crs_wkt :

    PROJCRS["WGS 84 / UTM zone 29N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 29N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-9,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Navigation and medium accuracy spatial referencing."],AREA["Between 12°W and 6°W, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Côte D'Ivoire (Ivory Coast). Faroe Islands. Guinea. Ireland. Jan Mayen. Liberia, Mali. Mauritania. Morocco. Portugal. Sierra Leone. Spain. United Kingdom (UK). Western Sahara."],BBOX[0,-12.01,84.01,-6]],ID["EPSG",32629]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    projected_crs_name :

    WGS 84 / UTM zone 29N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -9.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    long_name :

    CRS definition

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 29N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32629"]]

    array(0, dtype=int8)

  • turbine_id
    (point)
    int64
    0 1 2 3 4 5 6 7 8 9 10 11 12 13 14

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])

  • group_id
    (point)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

  • wtg_key
    (point)
    <U10
    'Bonus_1_MW' ... 'Bonus_1_MW'

    array(['Bonus_1_MW', 'Bonus_1_MW', 'Bonus_1_MW', 'Bonus_1_MW',
           'Bonus_1_MW', 'Bonus_1_MW', 'Bonus_1_MW', 'Bonus_1_MW',
           'Bonus_1_MW', 'Bonus_1_MW', 'Bonus_1_MW', 'Bonus_1_MW',
           'Bonus_1_MW', 'Bonus_1_MW', 'Bonus_1_MW'], dtype='<U10')

• Data variables: (0)
• Indexes: (0)
• Attributes: (2)

  Conventions :

  CF-1.8

  history :

  2025-07-04T13:07:21+00:00: windkit==1.0.0 create_dataset(

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

ax = wk.plot.elevation_map(elev_map, zorder=0)
wk.spatial.ds_to_gdf(bwc).plot(ax=ax, color="red", marker="X", markersize=40, zorder=1)
wk.spatial.ds_to_gdf(wind_turbines).plot(ax=ax, color="black", marker="X", markersize=40, zorder=1)
ax.set(xlim=(511000, 518000), ylim=(4620000, 4626000))

[(511000.0, 518000.0), (4620000.0, 4626000.0)]

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

pwc = pw.wasp.downscale(gwc, topo_map, wind_turbines, interp_method="nearest")

Now we can calculate the AEP

aep_nowake = pw.gross_aep(pwc, wtg)

aep_nowake = aep_nowake["gross_aep"]
print(aep_nowake)
print(aep_nowake.sum())

<xarray.DataArray 'gross_aep' (point: 15)> Size: 60B
array([2.8252254, 2.687702 , 2.4737535, 2.5607786, 2.727687 , 2.807531 ,
       2.6752436, 2.969323 , 3.060742 , 2.307256 , 2.6991627, 2.9334679,
       2.889413 , 3.0800772, 3.2932563], dtype=float32)
Coordinates:
    height       (point) int64 120B 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
    south_north  (point) int64 120B 4622313 4622199 4622336 ... 4624252 4624142
    west_east    (point) int64 120B 513914 514161 514425 ... 516060 516295
    crs          int8 1B 0
    turbine_id   (point) int64 120B 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
    group_id     (point) float64 120B 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    wtg_key      (point) <U10 600B 'Bonus_1_MW' 'Bonus_1_MW' ... 'Bonus_1_MW'
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' ()> Size: 4B
array(41.99062, dtype=float32)
Coordinates:
    crs      int8 1B 0

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

Notice that if we would not be correcting for air density, the 1.225 kg/m^3 reference power curve would be used, achieving a virtually less accurate result of 43.7 GWh annually. This can be checked by setting the air_density_correction = "none".

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.

aep_park2 = pw.potential_aep(
    pwc,
    wtg,
    wind_farm_model="PARK2_onshore",
    turbulence_intensity=0.0
)

aep_park2 = aep_park2["potential_aep"]
print(aep_park2)
print(aep_park2.sum())

<xarray.DataArray 'potential_aep' (point: 15)> Size: 120B
array([2.7278782 , 2.59931447, 2.36980802, 2.45435784, 2.62128109,
       2.76544897, 2.64794008, 2.92111536, 3.00184219, 2.23764838,
       2.64323762, 2.84922645, 2.74544182, 2.91545043, 3.17742836])
Coordinates:
    height       (point) int64 120B 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
    south_north  (point) int64 120B 4622313 4622199 4622336 ... 4624252 4624142
    west_east    (point) int64 120B 513914 514161 514425 ... 516060 516295
    crs          int8 1B 0
    turbine_id   (point) int64 120B 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
    group_id     (point) float64 120B 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    wtg_key      (point) <U10 600B 'Bonus_1_MW' 'Bonus_1_MW' ... 'Bonus_1_MW'
    mode         int64 8B 0
  * point        (point) int64 120B 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
Attributes:
    description:    Potential annual energy production, considering wake effe...
    long_name:      Potential annual energy production
    standard_name:  potential_aep
    units:          GWh
    grid_mapping:   crs
<xarray.DataArray 'potential_aep' ()> Size: 8B
array(40.67741931)
Coordinates:
    crs      int8 1B 0
    mode     int64 8B 0

We can see that wake effects causes a reduction of the AEP from 42 to 40.7 GWh/y (a 3.3 % reduction).

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

fig, axes = plt.subplots(1, 2, figsize=(15, 5))
aep_nowake.plot(ax=axes[0], marker="o", color="blue", label=f"No wakes: {aep_nowake.sum().values:.1f} [GWh]")
aep_park2.plot(ax=axes[0], marker="*", markersize=8, label=f"Wakes: {aep_park2.sum().values:.1f} [GWh]", color="green")

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

# Add grid to both plots
axes[0].grid(True, linestyle="--", alpha=0.7)
axes[1].grid(True, linestyle="--", alpha=0.7)

# Add subtitles to each plot
axes[0].set_title("Gross vs Potential AEP", fontsize=14)
axes[1].set_title("Ratio of Potential to Gross AEP", fontsize=14)

# Set x-axis and y-axis labels
axes[0].set(xlabel="Wind Turbine id", ylabel="AEP [GWh]", ylim=(2, 4))
axes[1].set(xlabel="Wind Turbine id", ylabel=r"AEP / AEP$_\mathrm{no wakes}$")

axes[0].legend(frameon=False, title="Wake model")

plt.show()

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

output_grid = wk.spatial.create_cuboid(
    west_east=np.linspace(511700.0, 517700.0, 41),
    south_north=np.linspace(4619500.0, 4625500.0, 41),
    height=[50.0],
    crs="EPSG:32629",
)

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

<xarray.Dataset> Size: 270kB
Dimensions:        (sector: 12, height: 1, south_north: 41, west_east: 41)
Coordinates:
    crs            int8 1B 0
    sector_ceil    (sector) float64 96B 15.0 45.0 75.0 ... 285.0 315.0 345.0
    sector_floor   (sector) float64 96B 345.0 15.0 45.0 ... 255.0 285.0 315.0
  * height         (height) float64 8B 50.0
  * south_north    (south_north) float64 328B 4.62e+06 4.62e+06 ... 4.626e+06
  * west_east      (west_east) float64 328B 5.117e+05 5.118e+05 ... 5.177e+05
  * sector         (sector) float64 96B 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
Data variables:
    A              (sector, height, south_north, west_east) float32 81kB 5.16...
    k              (sector, height, south_north, west_east) float32 81kB 1.98...
    wdfreq         (sector, height, south_north, west_east) float32 81kB 0.06...
    site_elev      (south_north, west_east) float32 7kB 15.3 22.06 ... 199.8
    air_density    (height, south_north, west_east) float32 7kB 1.216 ... 1.194
    wspd           (height, south_north, west_east) float32 7kB 5.545 ... 4.86
    power_density  (height, south_north, west_east) float32 7kB 200.1 ... 121.8
Attributes:
    Conventions:      CF-1.8
    history:          2025-07-04T13:07:45+00:00:\twindkit==1.0.0\t"cuboid")\n...
    Package name:     windkit
    Package version:  1.0.0
    Creation date:    2025-07-04T13:07:50+00:00
    Object type:      Met fields
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind
    title:            WAsP site effects

xarray.Dataset

• Dimensions:
  • sector: 12
  • height: 1
  • south_north: 41
  • west_east: 41
• Coordinates: (7)
  • crs
    ()
    int8
    0

    crs_wkt :

    PROJCRS["WGS 84 / UTM zone 29N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 29N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-9,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Navigation and medium accuracy spatial referencing."],AREA["Between 12°W and 6°W, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Côte D'Ivoire (Ivory Coast). Faroe Islands. Guinea. Ireland. Jan Mayen. Liberia, Mali. Mauritania. Morocco. Portugal. Sierra Leone. Spain. United Kingdom (UK). Western Sahara."],BBOX[0,-12.01,84.01,-6]],ID["EPSG",32629]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    projected_crs_name :

    WGS 84 / UTM zone 29N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -9.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    long_name :

    CRS definition

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 29N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32629"]]

    array(0, dtype=int8)

  • sector_ceil
    (sector)
    float64
    15.0 45.0 75.0 ... 315.0 345.0

    description :

    The sector ceiling is the bin edge value in clockwise direction.

    long_name :

    Wind direction binned into sectors (clockwise edge value)

    standard_name :

    wind_from_direction

    units :

    degree

    array([ 15.,  45.,  75., 105., 135., 165., 195., 225., 255., 285., 315.,
           345.])

  • sector_floor
    (sector)
    float64
    345.0 15.0 45.0 ... 285.0 315.0

    description :

    The sector floor is the bin edge value in counter-clockwise direction.

    long_name :

    Wind direction binned into sectors (counter-clockwise edge value)

    standard_name :

    wind_from_direction

    units :

    degree

    array([345.,  15.,  45.,  75., 105., 135., 165., 195., 225., 255., 285.,
           315.])

  • height
    (height)
    float64
    50.0

    URI :

    axis :

    Z

    description :

    It is the height above ground either related to the measurements or to the hub of a turbine.

    long_name :

    Height above ground

    standard_name :

    height

    units :

    m

    array([50.])

  • south_north
    (south_north)
    float64
    4.62e+06 4.62e+06 ... 4.626e+06

    URI :

    axis :

    Y

    description :

    Positive northward direction in projected coordinates (meters).

    long_name :

    Projected south to north coordinate

    standard_name :

    projection_y_coordinate

    units :

    m

    array([4619500., 4619650., 4619800., 4619950., 4620100., 4620250., 4620400.,
           4620550., 4620700., 4620850., 4621000., 4621150., 4621300., 4621450.,
           4621600., 4621750., 4621900., 4622050., 4622200., 4622350., 4622500.,
           4622650., 4622800., 4622950., 4623100., 4623250., 4623400., 4623550.,
           4623700., 4623850., 4624000., 4624150., 4624300., 4624450., 4624600.,
           4624750., 4624900., 4625050., 4625200., 4625350., 4625500.])

  • west_east
    (west_east)
    float64
    5.117e+05 5.118e+05 ... 5.177e+05

    URI :

    axis :

    X

    description :

    Positive eastward direction in projected coordinates (meters).

    long_name :

    Projected west to east coordinate

    standard_name :

    projection_x_coordinate

    units :

    m

    array([511700., 511850., 512000., 512150., 512300., 512450., 512600., 512750.,
           512900., 513050., 513200., 513350., 513500., 513650., 513800., 513950.,
           514100., 514250., 514400., 514550., 514700., 514850., 515000., 515150.,
           515300., 515450., 515600., 515750., 515900., 516050., 516200., 516350.,
           516500., 516650., 516800., 516950., 517100., 517250., 517400., 517550.,
           517700.])

  • sector
    (sector)
    float64
    0.0 30.0 60.0 ... 270.0 300.0 330.0

    ancillary_variables :

    sector_ceil sector_floor

    description :

    A sector identifies a given range of wind directions binned around a central value , in a wind rose, to classify where the wind blows from.

    long_name :

    Wind direction binned into sectors (center value)

    standard_name :

    wind_from_direction

    units :

    degree

    array([  0.,  30.,  60.,  90., 120., 150., 180., 210., 240., 270., 300., 330.])

• Data variables: (7)
  • A
    (sector, height, south_north, west_east)
    float32
    5.162 5.056 5.091 ... 5.384 5.319

    grid_mapping :

    crs

    description :

    Scale parameter in the Weibull distribution

    long_name :

    Weibull-scale parameter of wind distribution

    standard_name :

    weibull_scale_parameter

    units :

    m s-1

    array([[[[5.1624684, 5.0558662, 5.0908823, ..., 4.2981634, 4.52969  ,
              4.747555 ],
             [5.1882067, 5.135966 , 5.099238 , ..., 4.4492817, 4.590198 ,
              4.7172446],
             [5.1291513, 4.9322596, 4.595451 , ..., 4.8833656, 4.6424494,
              4.640777 ],
             ...,
             [4.8705606, 4.8142333, 4.854361 , ..., 5.3254066, 5.485707 ,
              5.333031 ],
             [5.0088205, 4.840428 , 4.8644543, ..., 5.2730603, 5.1686344,
              5.2502704],
             [5.0491605, 4.885099 , 4.9733753, ..., 5.2756877, 5.176761 ,
              5.10841  ]]],
    
    
           [[[3.9357603, 3.7330105, 3.9756203, ..., 5.127866 , 5.2260065,
              5.1887474],
             [3.982159 , 3.9275534, 4.064692 , ..., 5.2337103, 5.267398 ,
              5.1806703],
             [3.9391415, 3.793193 , 3.7050915, ..., 5.3742394, 5.1398063,
    ...
              5.508224 ],
             [6.7245917, 6.6038337, 6.334707 , ..., 6.246952 , 5.694558 ,
              5.2830067],
             [6.912048 , 6.9017944, 6.3230877, ..., 6.540055 , 5.4443455,
              5.1552677]]],
    
    
           [[[6.0952516, 6.1244087, 6.2940702, ..., 3.5448468, 3.7971709,
              3.9457622],
             [6.057892 , 6.0485373, 6.1828423, ..., 3.5704257, 3.763005 ,
              3.842969 ],
             [5.995926 , 5.8673897, 5.7738233, ..., 4.0823226, 3.7506669,
              3.7756565],
             ...,
             [5.5731897, 5.372512 , 5.146339 , ..., 5.609263 , 5.7461915,
              5.702161 ],
             [5.718336 , 5.360762 , 5.038505 , ..., 5.5603123, 5.4257703,
              5.501515 ],
             [5.857077 , 5.449075 , 5.093078 , ..., 5.808359 , 5.3837748,
              5.318997 ]]]], shape=(12, 1, 41, 41), dtype=float32)

  • k
    (sector, height, south_north, west_east)
    float32
    1.982 1.975 1.994 ... 2.143 2.154

    grid_mapping :

    crs

    description :

    Form parameter in the Weibull distribution identify the shape of the distribution. This can be connected to characteristics of theterrain surrounding themeasuring masts

    long_name :

    Weibull-shape parameter of wind distribution

    standard_name :

    weibull_shape_parameter

    array([[[[1.9824219, 1.9746094, 1.9941406, ..., 1.9980469, 2.0214844,
              2.0214844],
             [1.9863281, 1.9941406, 2.0058594, ..., 1.9941406, 2.0097656,
              2.0175781],
             [1.9824219, 1.9863281, 2.0019531, ..., 2.0292969, 2.0136719,
              2.0214844],
             ...,
             [2.0058594, 2.0058594, 2.0097656, ..., 2.0722656, 2.0839844,
              2.0605469],
             [2.0136719, 2.0136719, 2.0097656, ..., 2.0683594, 2.0839844,
              2.0722656],
             [2.0253906, 2.0136719, 2.0097656, ..., 2.0800781, 2.0839844,
              2.0722656]]],
    
    
           [[[2.3378906, 2.2832031, 2.3261719, ..., 2.2792969, 2.2871094,
              2.2910156],
             [2.3378906, 2.3066406, 2.3300781, ..., 2.2753906, 2.2792969,
              2.2792969],
             [2.3417969, 2.3378906, 2.3691406, ..., 2.3066406, 2.2832031,
    ...
              2.3144531],
             [2.1738281, 2.1777344, 2.1777344, ..., 2.2832031, 2.3066406,
              2.3144531],
             [2.1855469, 2.1699219, 2.1777344, ..., 2.2832031, 2.3105469,
              2.3183594]]],
    
    
           [[[2.0449219, 2.0527344, 2.0410156, ..., 2.1269531, 2.1308594,
              2.1308594],
             [2.0488281, 2.0605469, 2.0488281, ..., 2.0917969, 2.1113281,
              2.1074219],
             [2.0488281, 2.0527344, 2.0253906, ..., 2.1152344, 2.0878906,
              2.0996094],
             ...,
             [2.0214844, 2.0253906, 2.0292969, ..., 2.1230469, 2.1269531,
              2.1542969],
             [2.0214844, 2.0253906, 2.0332031, ..., 2.1230469, 2.1347656,
              2.1542969],
             [2.0175781, 2.0214844, 2.0292969, ..., 2.1152344, 2.1425781,
              2.1542969]]]], shape=(12, 1, 41, 41), dtype=float32)

  • wdfreq
    (sector, height, south_north, west_east)
    float32
    0.06677 0.06735 ... 0.08941 0.09854

    grid_mapping :

    crs

    description :

    Frequency of occurance for each wind direction bin/sector

    long_name :

    Relative frequency per wind_direction over time and wind_speed

    short_name :

    wd_bin_freq

    units :

    `1`

    array([[[[0.06677344, 0.06735355, 0.06429566, ..., 0.07531703,
              0.07711753, 0.07903025],
             [0.06821653, 0.07000158, 0.06643819, ..., 0.0792757 ,
              0.07905921, 0.0824925 ],
             [0.06841841, 0.0673871 , 0.06175123, ..., 0.08146112,
              0.08399127, 0.08369356],
             ...,
             [0.06391951, 0.06419186, 0.06408658, ..., 0.05682411,
              0.06040112, 0.06437391],
             [0.06434254, 0.06300002, 0.06335249, ..., 0.05644369,
              0.06141147, 0.0661058 ],
             [0.06344649, 0.06073633, 0.06282016, ..., 0.0549212 ,
              0.06376923, 0.06569072]]],
    
    
           [[[0.03778108, 0.03675805, 0.03589549, ..., 0.0665926 ,
              0.06245298, 0.05904483],
             [0.03832453, 0.0385283 , 0.03725101, ..., 0.06977993,
              0.06499024, 0.06153809],
             [0.03823056, 0.03740584, 0.03485667, ..., 0.06137894,
    ...
              0.09529064, 0.09680722],
             [0.09163941, 0.09293695, 0.08814103, ..., 0.09231574,
              0.09553654, 0.09412713],
             [0.09763809, 0.09579047, 0.08747481, ..., 0.10002708,
              0.0937321 , 0.09409786]]],
    
    
           [[[0.09683747, 0.10490197, 0.10275893, ..., 0.06394778,
              0.06743169, 0.07070699],
             [0.09740739, 0.1007514 , 0.10036812, ..., 0.06418205,
              0.06655344, 0.06996541],
             [0.09812856, 0.10043584, 0.10039401, ..., 0.06995085,
              0.06867456, 0.07023737],
             ...,
             [0.07496112, 0.07486518, 0.07389177, ..., 0.07450692,
              0.08205397, 0.1030598 ],
             [0.0784072 , 0.07647292, 0.07182025, ..., 0.07321104,
              0.08583056, 0.10216461],
             [0.08414437, 0.07578106, 0.0704685 , ..., 0.07674919,
              0.08941462, 0.09853814]]]], shape=(12, 1, 41, 41), dtype=float32)

  • site_elev
    (south_north, west_east)
    float32
    15.3 22.06 47.66 ... 215.2 199.8

    description :

    height above sea level of a point on the surface of an orographic element e.g. hill, mountain

    long_name :

    Terrain Elevation

    standard_name :

    ground_level_altitude

    units :

    m

    grid_mapping :

    crs

    array([[ 15.302768,  22.05551 ,  47.655792, ...,  35.67707 ,  33.988403,
             32.347656],
           [ 17.60524 ,  28.119587,  52.118202, ...,  44.222023,  39.728123,
             33.013466],
           [ 18.31766 ,  27.830061,  46.473217, ...,  63.516724,  41.586285,
             32.80675 ],
           ...,
           [ 39.313004,  43.326805,  51.603905, ..., 266.78992 , 257.91586 ,
            228.58943 ],
           [ 45.808395,  50.034653,  58.51156 , ..., 258.95557 , 230.7904  ,
            212.71527 ],
           [ 49.376396,  62.463497,  70.96958 , ..., 250.40358 , 215.18497 ,
            199.7564  ]], shape=(41, 41), dtype=float32)

  • air_density
    (height, south_north, west_east)
    float32
    1.216 1.215 1.212 ... 1.192 1.194

    source_file :

    10yr_mean_rhoinp_ERA5v1.nc

    grid_mapping :

    crs

    description :

    The density of air or atmospheric density, is the mass per unit volume of Earth's atmosphere.

    long_name :

    Air Density

    standard_name :

    air_density

    units :

    kg m-3

    array([[[1.2157065, 1.2148893, 1.2117953, ..., 1.2132416, 1.2134461,
             1.2136444],
            [1.2154275, 1.2141559, 1.2112566, ..., 1.2122098, 1.2127527,
             1.2135637],
            [1.2153418, 1.2141908, 1.2119379, ..., 1.2098808, 1.212528 ,
             1.2135886],
            ...,
            [1.2128024, 1.2123182, 1.2113181, ..., 1.1855681, 1.1866215,
             1.190108 ],
            [1.2120183, 1.2115076, 1.2104847, ..., 1.1864989, 1.1898451,
             1.1919981],
            [1.2115873, 1.210008 , 1.2089827, ..., 1.1875134, 1.1917033,
             1.1935421]]], shape=(1, 41, 41), dtype=float32)

  • wspd
    (height, south_north, west_east)
    float32
    5.545 5.387 5.6 ... 5.132 4.86

    description :

    Speed is the magnitude of velocity. Wind is defined as a two-dimensional (horizontal) air velocity vector, with no vertical component. (Vertical motion in the atmosphere has the standard name upward_air_velocity.) The wind speed is the magnitude of the wind velocity.

    long_name :

    Horizontal Wind Speed

    standard_name :

    wind_speed

    units :

    m s-1

    grid_mapping :

    crs

    array([[[5.545009 , 5.387219 , 5.6002545, ..., 4.6217294, 4.759212 ,
             4.80211  ],
            [5.5077157, 5.392588 , 5.53437  , ..., 4.618456 , 4.7228193,
             4.699772 ],
            [5.4340277, 5.2657948, 5.1993637, ..., 4.8671308, 4.60993  ,
             4.56762  ],
            ...,
            [5.6410613, 5.539159 , 5.469485 , ..., 5.960794 , 5.752547 ,
             5.1193304],
            [5.7139297, 5.5881677, 5.549388 , ..., 5.9944935, 5.3082542,
             4.9745317],
            [5.6829877, 5.763244 , 5.718164 , ..., 5.974998 , 5.131796 ,
             4.8600335]]], shape=(1, 41, 41), dtype=float32)

  • power_density
    (height, south_north, west_east)
    float32
    200.1 183.5 200.5 ... 140.0 121.8

    description :

    Power density is the amount of power available in the flow per unit surface

    long_name :

    Mean power density of the wind speed

    units :

    W m-2

    grid_mapping :

    crs

    array([[[200.06717, 183.51259, 200.50394, ..., 122.31127, 132.90495,
             135.83174],
            [196.03445, 184.0504 , 193.67836, ..., 126.27257, 132.30315,
             131.11972],
            [187.6766 , 168.87952, 158.93053, ..., 142.58034, 125.56654,
             121.37914],
            ...,
            [200.03613, 188.8858 , 181.25926, ..., 215.68974, 194.21101,
             144.09125],
            [206.57803, 193.27618, 191.0774 , ..., 220.35808, 153.65952,
             132.21971],
            [202.02501, 211.55939, 209.46297, ..., 216.7571 , 139.98589,
             121.83944]]], shape=(1, 41, 41), dtype=float32)

• Indexes: (4)
  • height
    PandasIndex

    PandasIndex(Index([50.0], dtype='float64', name='height'))

  • south_north
    PandasIndex

    PandasIndex(Index([4619500.0, 4619650.0, 4619800.0, 4619950.0, 4620100.0, 4620250.0,
           4620400.0, 4620550.0, 4620700.0, 4620850.0, 4621000.0, 4621150.0,
           4621300.0, 4621450.0, 4621600.0, 4621750.0, 4621900.0, 4622050.0,
           4622200.0, 4622350.0, 4622500.0, 4622650.0, 4622800.0, 4622950.0,
           4623100.0, 4623250.0, 4623400.0, 4623550.0, 4623700.0, 4623850.0,
           4624000.0, 4624150.0, 4624300.0, 4624450.0, 4624600.0, 4624750.0,
           4624900.0, 4625050.0, 4625200.0, 4625350.0, 4625500.0],
          dtype='float64', name='south_north'))

  • west_east
    PandasIndex

    PandasIndex(Index([511700.0, 511850.0, 512000.0, 512150.0, 512300.0, 512450.0, 512600.0,
           512750.0, 512900.0, 513050.0, 513200.0, 513350.0, 513500.0, 513650.0,
           513800.0, 513950.0, 514100.0, 514250.0, 514400.0, 514550.0, 514700.0,
           514850.0, 515000.0, 515150.0, 515300.0, 515450.0, 515600.0, 515750.0,
           515900.0, 516050.0, 516200.0, 516350.0, 516500.0, 516650.0, 516800.0,
           516950.0, 517100.0, 517250.0, 517400.0, 517550.0, 517700.0],
          dtype='float64', name='west_east'))

  • sector
    PandasIndex

    PandasIndex(Index([0.0, 30.0, 60.0, 90.0, 120.0, 150.0, 180.0, 210.0, 240.0, 270.0, 300.0,
           330.0],
          dtype='float64', name='sector'))

• Attributes: (10)

  Conventions :

  CF-1.8

  history :

  2025-07-04T13:07:45+00:00: windkit==1.0.0 "cuboid") 2025-07-04T13:07:45+00:00: windkit==1.0.0 to_point(source) 2025-07-04T13:07:45+00:00: windkit==1.0.0 _spatial_stack(output_locs, self.elev_map.crs).load() 2025-07-04T13:07:48+00:00: windkit==1.0.0 to_raster(source, True) 2025-07-04T13:07:48+00:00: windkit==1.0.0 _spatial_unstack(work_locs) 2025-07-04T13:07:48+00:00: windkit==1.0.0 topo_map.get_site_effects(output_locs, conf, n_sectors=n_sectors) 2025-07-04T13:07:48+00:00: windkit==1.0.0 to_point(source) 2025-07-04T13:07:48+00:00: windkit==1.0.0 wksp._spatial_stack(site_effects) 2025-07-04T13:07:49+00:00: windkit==1.0.0 to_raster(source, True) 2025-07-04T13:07:49+00:00: windkit==1.0.0 wksp._spatial_unstack(wc) 2025-07-04T13:07:50+00:00: windkit==1.0.0 add_met_fields(wc) 2025-07-04T13:07:50+00:00: windkit==1.0.0 func(*args, **kwargs) 2025-07-04T13:07:50+00:00: windkit==1.0.0 func(*args, **kwargs)

  Package name :

  windkit

  Package version :

  1.0.0

  Creation date :

  2025-07-04T13:07:50+00:00

  Object type :

  Met fields

  author :

  Bjarke Tobias Olsen

  author_email :

  btol@dtu.dk

  institution :

  DTU Wind

  title :

  WAsP site effects

flow_map = pw.wind_farm_flow_map(
    pwc_grid,
    wtg,
    wind_turbines,
    output_grid,
    wind_farm_model="PARK2_onshore",
    turbulence_intensity=0.0,
    n_subsector=5,
)
print(flow_map)

<xarray.Dataset> Size: 1MB
Dimensions:                          (sector: 12, height: 1, south_north: 41,
                                      west_east: 41)
Coordinates:
    crs                              int8 1B 0
    mode                             int64 8B 0
  * sector                           (sector) float64 96B 0.0 30.0 ... 330.0
  * height                           (height) float64 8B 50.0
  * south_north                      (south_north) float64 328B 4.62e+06 ... ...
  * west_east                        (west_east) float64 328B 5.117e+05 ... 5...
Data variables: (12/15)
    potential_aep_sector             (sector, height, south_north, west_east) float64 161kB ...
    gross_aep_sector                 (sector, height, south_north, west_east) float64 161kB ...
    potential_aep_deficit_sector     (sector, height, south_north, west_east) float64 161kB ...
    wspd_sector                      (sector, height, south_north, west_east) float64 161kB ...
    wspd_eff_sector                  (sector, height, south_north, west_east) float64 161kB ...
    wspd_deficit_sector              (sector, height, south_north, west_east) float64 161kB ...
    ...                               ...
    gross_aep                        (height, south_north, west_east) float64 13kB ...
    potential_aep_deficit            (height, south_north, west_east) float64 13kB ...
    wspd                             (height, south_north, west_east) float64 13kB ...
    wspd_eff                         (height, south_north, west_east) float64 13kB ...
    wspd_deficit                     (height, south_north, west_east) float64 13kB ...
    turbulence_intensity_eff         (height, south_north, west_east) float64 13kB ...
Attributes:
    Conventions:      CF-1.8
    history:          2025-07-04T13:07:45+00:00:\twindkit==1.0.0\t"cuboid")\n...
    Package name:     windkit
    Package version:  1.0.0
    Creation date:    2025-07-04T13:08:02+00:00
    Object type:      Wind Farm Flow Map
    author:           Bjarke Tobias Olsen
    author_email:     btol@dtu.dk
    institution:      DTU Wind
    title:            WAsP site effects

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(wind_turbines.west_east, wind_turbines.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()

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["potential_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["potential_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()