Site Object

For a given position, reference wind speed ws and reference wind direction wd or time, Site provides the local wind condition in terms of local wind speed WS, local wind direction WD, local turbulence intensity TI and the probability of each flow case (wd-ws combination or time). Furthermore, Site contains a Distance object responsible for calculating the down-wind, cross-wind and vertical distance between wind turbines. The distances can be either straight-line, terrain-elevation-following or streamline-following distances.

[1]:

import os
from IPython.display import HTML
from matplotlib.animation import FuncAnimation, PillowWriter
from numpy import newaxis as na
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
import numpy as np
from py_wake import NOJ, NOJLocal
from py_wake.examples.data.ParqueFicticio import ParqueFicticioSite, ParqueFicticio_path
from py_wake.examples.data.dtu10mw._dtu10mw import DTU10MW
from py_wake.examples.data.hornsrev1 import Hornsrev1Site, Hornsrev1WRFSite, V80
from py_wake.examples.data.iea37 import IEA37Site
from py_wake.flow_map import XYGrid
from py_wake.site import UniformSite, UniformWeibullSite, WaspGridSite, XRSite
from py_wake.site.distance import StraightDistance
from py_wake.site.shear import PowerShear
from py_wake.utils import layouts
from py_wake.utils import weibull
from py_wake.utils.maps import dk_coast
from py_wake.utils.plotting import setup_plot
from py_wake.utils.functions import mean_deg
import xarray as xr
xr.set_options(display_expand_data=True);

Example sites

PyWake has a few predefined demonstration sites:

IEA37Site: UniformSite (fix wind speed (9.8m/s), predefined wind sector probability).
Hornsrev1: UniformWeibullSite (Weibull distributed wind speed, predefined wind sector propability, uniform wind a over flat wind area).
ParqueFicticioSite: WaspGridSite (position-dependent Weibull distributed wind speed and sector probability. Terrain following distances over non-flat terrain). Loaded from a set of *.grd files exported from WAsP.
Hornsrev1WRF: Time and position dependent wind speed and direction. The site covers Hornsrev1, April 2020. Streamline-following distances

[2]:

sites = {"IEA37": IEA37Site(n_wt=16),
         "Hornsrev1": Hornsrev1Site(),
         "ParqueFicticio": ParqueFicticioSite(),
         "Hornsrev1WRF": Hornsrev1WRFSite()}

Downloading data from 'https://gitlab.windenergy.dtu.dk/api/v4/projects/3385/repository/files/wrf%2FHornsrev1%5F2020%5F04%2Enc/raw' to file '/root/PyWakeTestFiles/f57bffa22370637a3f5f3b67970fcc96-raw'.

Configuration saved to /root/.config/windkit/windkit_config.toml.

User-defined sites

You can define your own site using one of the Site classes:

UniformWeibullSite: Site with uniform sector-dependent Weibull distributed wind speed.
XRSite: The flexible general base class behind all Sites.
WaspGridSite: Site with gridded non-uniform inflow based on *.grd files exported from WAsP.
WRFSite: A subclass of XRSite which take a WRF netcdf as input containing time and position dependent WS, WD and TKE/TI

For more information on these classes, please see the API reference on the Site object.

XRSite

The XRSite is the most general and flexible Site. For the input dataset there are some required and optional data variables, such as:

Required data variables:
- P: probability of flow case(s)
or
- Weibull_A: Weibull scale parameter(s)
- Weibull_k: Weibull shape parameter(s)
- Sector_frequency: Probability of each wind direction sector
Optional data variables:
- WS: Wind speed, if not present, the reference wind speed ws is used
- Speedup: Factor multiplied to the wind speed
- Turning: Wind direction turning
- TI: Turbulence intensity
- xxx: Custom variables needed by the wind turbines to compute power, ct or loads
Each data variable may be constant or depend on a combination of the following inputs (Note, the input variables must be ordered according to the list, i.e. P(wd,ws) is ok, while P(ws,wd) is not):
- i: Wind turbine position (one position per wind turbine)
- x,y: Gridded 2d position
- x,y,h: Gridded 3d position
- time: Time
- wd: Refernce wind direction
- ws : Reference wind speed

[6]:

f = [0.036, 0.039, 0.052, 0.07, 0.084, 0.064, 0.086, 0.118, 0.152, 0.147, 0.1, 0.052]
A = [9.177, 9.782, 9.532, 9.91, 10.043, 9.594, 9.584, 10.515, 11.399, 11.687, 11.637, 10.088]
k = [2.393, 2.447, 2.412, 2.592, 2.756, 2.596, 2.584, 2.549, 2.471, 2.607, 2.627, 2.326]
wd = np.linspace(0, 360, len(f), endpoint=False)
ti = .1

Uniform site

[7]:

# Site with constant wind speed, sector frequency, constant turbulence intensity and power shear
uniform_site = XRSite(
    ds=xr.Dataset(data_vars={'WS': 10, 'P': ('wd', f), 'TI': ti},
                  coords={'wd': wd}),
    shear=PowerShear(h_ref=100, alpha=.2))

Wind direction dependent weibull distributed wind speed

[8]:

# Site with wind direction dependent weibull distributed wind speed
uniform_weibull_site = XRSite(
    ds=xr.Dataset(data_vars={'Sector_frequency': ('wd', f), 'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k), 'TI': ti},
                  coords={'wd': wd}))

Wind turbine dependent speedup and turning

[9]:

# Site with a speedup and a turning value per WT
x_i, y_i = np.arange(5) * 100, np.zeros(5)  # WT positions

complex_fixed_pos_site = XRSite(
    ds=xr.Dataset(
        data_vars={'Speedup': ('i', np.arange(.8, 1.3, .1)),
                   'Turning': ('i', np.arange(-2, 3)),
                   'P': ('wd', f)},
        coords={'i': np.arange(5), 'wd': wd}),
    initial_position=np.array([x_i, y_i]).T)

Position dependent speedup

[10]:

# Site with gridded speedup information
complex_grid_site = XRSite(
    ds=xr.Dataset(
        data_vars={'Speedup': (['x', 'y'], np.arange(.8, 1.4, .1).reshape((3, 2))),
                   'P': ('wd', f)},
        coords={'x': [0, 500, 1000], 'y': [0, 500], 'wd': wd}))

WS dependent speedup, custom (wd,ws)-dependent propability

[11]:

# Site with ws dependent speedup and wd- and ws distributed probability
P_ws = weibull.cdf(np.array([3, 5, 7, 9, 11, 13]), 10, 2) - weibull.cdf(np.array([0, 3, 5, 7, 9, 11]), 10, 2)
P_wd_ws = P_ws[na, :] * np.array(f)[:, na]

complex_ws_site = XRSite(
    ds=xr.Dataset(
        data_vars={'Speedup': (['ws'], np.arange(.8, 1.4, .1)),
                   'P': (('wd', 'ws'), P_wd_ws), 'TI': ti},
        coords={'ws': [1.5, 4, 6, 8, 10, 12], 'wd': wd}))

Position and time dependent wind resources

This example creates a site with 2D non-uniform time series resources for wind speed (ws), wind direction (wd), and turbulence intensity (ti). The data is provided as arrays with dimensions corresponding to spatial coordinates (x, y) and time.

[12]:

site_x, site_y = np.meshgrid(np.arange(0.1, 1000, 100), np.arange(0.1, 2000, 100))
site_x, site_y = site_x.flatten(), site_y.flatten()
site_time = np.arange(100)
site_ws = np.random.uniform(3.0, 21.0, (len(site_x), len(site_y), len(site_time)))
site_wd = np.random.uniform(0.0, 360.0, (len(site_x), len(site_y), len(site_time)))
ds = xr.Dataset(
    data_vars=dict(
        WS=(["x", "y", "time"], site_ws),
        WD=(["x", "y", "time"], site_wd),
        TI=(["x", "y", "time"], np.ones_like(site_ws) * 0.1),  # hardcoded TI=0.1
        P=1,  # deterministic wind resource
    ),
    coords=dict(
        x=("x", site_x),
        y=("y", site_y),
        time=("time", site_time),
    ),
)
non_uniform_ts_site = XRSite(ds)
ws = site_ws.mean((0,1)) # spatial mean wind speed
wd = mean_deg(site_wd, (0,1)) # spatial mean wind direction (angular mean)
wss_at_mean_loc = non_uniform_ts_site.local_wind(
    site_x.mean(), site_y.mean(), ws=ws, wd=wd, time=site_time
)["WS_ilk"]
print(f"Mean wind speed at the mean location: {wss_at_mean_loc.mean():.2f} m/s")

# check the map of the mean wind speed at the site
mean_resource = ds.WS.mean(dim="time").values
plt.contourf(site_x, site_y, mean_resource)
plt.colorbar()
plt.title("Mean wind speed [m/s]")
plt.xlabel("x [m]")
_ = plt.ylabel("y [m]")

Mean wind speed at the mean location: 11.84 m/s

Local wind speed vs. position dependent weibull distribution

Local wind speed, WS, can be specified in absolute values WS(x,y,h,ws) or as local a Speedup(x,y,h) factor that is multiplied to the reference wind speed ws. Using this approach, wind turbines will see the specified local wind speed, WS subtracted the deficit from other wind turbines.

An alternative tempting approach is to use position-dependent weibull parameters. This means that wind turbines will see the reference wind speed, ws subtracted the deficit from other wind turbines, but the weight of each flow case when summing up the AEP will be position dependent. This approach, however, is mixing up flow cases and should be avoided, see example below.

Especially, it is important to avoid using both approaches as the local variation will be taken into accound twice.

Assume we have the site with a constant reference wind speed, ws=9m/s and wind direction wd=270deg and two wind turbine as shown below.

The local wind speed without wakes, WS, is seen to vary due to lcoal Speedup and the deficit from the upstream turbine subtracts 3.1m/s from the effective wind speed, WS_eff of the downstream turbine.

[24]:

levels=np.linspace(0,15,31)
plt.figure(figsize=(6,3))
site = WaspGridSite.from_wasp_grd(ParqueFicticio_path, distance=StraightDistance())
sim_res = NOJLocal(site, V80())([263700, 264000], [6505650, 6505650], ws=9, wd=270, TI=0.1)
sim_res.flow_map(XYGrid(extend=0.3)).plot_wake_map(levels=levels, plot_windturbines=False)
sim_res.windFarmModel.windTurbines.plot(sim_res.x, sim_res.y,wt_number=False, wd=270)
def get_text(sim_res):
    return [f'''ws: {sr.ws.item():.1f}m/s
WS: {sr.WS.item():.1f}m/s
deficit: {(sr.WS-sr.WS_eff).item():.1f}m/s
WS_eff: {sr.WS_eff:.1f}m/s
Power: {sr.Power.item()*1e-6:.1f}MW
''' for sr in [sim_res.isel(wt=i,wd=0,ws=0) for i in [0,1]]]
s = get_text(sim_res)
plt.annotate(s[0], (sim_res.x[0]+10, sim_res.y[0]))
plt.annotate(s[1], (sim_res.x[1]+10, sim_res.y[1]));

If this flowcase was modelled with distributed local flow-case probabilities instead of local speedup, the local flow-case probabilities would mean that the upstream wind turbine sees the reference wind speed, ws=9.7m/s (left plot below), while the downstream wind turbine sees the refernce wind speed, ws=13.0m/s (right plot below).

The local speedup is thereby already taken into account, and the local wind speed WS equals the reference wind speed. The deficit from the upstream turbine seeing 9.7m/s is therefore also 3.1m/s as it should be. So far, everything matches the situation seen above.

The power of the downstream turbine, however, is wrong as it is extracted from the situation in the right plot below. The problem is that the deficit from the upstream wind turbine is calculated for an inflow wind speed of 13m/s instead of 9.7m/s. In other words, the results of the individual wind turbines can be modelled correctly with the distributed probability approach, but the interaction becomes wrong.

[25]:

axes = plt.subplots(1,2, figsize=(12,3))[1]
for i, (ax, ws) in enumerate(zip(axes, sim_res.WS.values)):
    sim_res_uni = NOJLocal(UniformSite(), V80())(sim_res.x, sim_res.y, ws=ws, wd=270, TI=0.1)
    sim_res_uni.flow_map(XYGrid(extend=0.3)).plot_wake_map(ax=ax, levels=levels, plot_windturbines=False)
    sim_res_uni.windFarmModel.windTurbines.plot(sim_res.x, sim_res.y,wt_number=False, wd=270, ax=ax)
    s = get_text(sim_res_uni)
    ax.annotate(s[0], (sim_res.x[0]+10, sim_res.y[0]), alpha=[1,0.4][i])
    ax.annotate(s[1], (sim_res.x[1]+10, sim_res.y[1]), alpha=[0.4,1][i]);

Distance

The wake effects depends on the downwind, crosswind and vertical distances between wind turbines. These distances are calculated by a Distance object, which is an input to the Site. There are currently three different distance classes

StraightDistance: Calculate the straight-line distances between wind turbines. Default distance object for most flat sites, e.g. IEA37Site and the Hornsrev1Site
TerrainFollwingDistance: Calculate the down-wind distance along a line with constant height above ground. Default distance object for WaspGridSite e.g. ParqueFicticioSite.
StreamlineDistance: Calcualte the downwind, crosswind and vertical distances along streamlines. Default for WRFSite (if the streamline argument is set to True), e.g. Hornsrev1WRFSite

[26]:

wd = [0, 30,90] # wind direction at source

for name, site in sites.items():
    print ("------- %s -------"%name)
    wt_x, wt_y = site.initial_position[0]

    dw_ijlk, cw_ijlk, dh_ijlk = site.distance(src_x_ilk=[wt_x], src_y_ilk=[wt_y], src_h_ilk=[70], wd_l=wd,
                                              time=[1281, 1281, 1281], # all 3 wd calculated at same time step
                                              dst_xyh_jlk=([wt_x], [wt_y-2000], [90]))

    print ('Wind direction: \t\t%d deg\t\t%d deg\t\t%d deg'%tuple(wd))
    print ('Down wind distance [m]: \t%.1f\t\t%.1f\t\t%.1f'%tuple(dw_ijlk[0,0,:,0]))
    print ('Cross wind distance [m]: \t%.1f\t\t%.1f\t\t%.1f'%tuple(cw_ijlk[0,0,:,0]))
    dh_ijlk = np.broadcast_to(dh_ijlk, dw_ijlk.shape)
    print ('Height difference [m]: \t\t%.1f\t\t%.1f\t\t%.1f'%tuple(dh_ijlk[0,0,:,0]))
    print()

------- IEA37 -------
Wind direction:                 0 deg           30 deg          90 deg
Down wind distance [m]:         2000.0          1732.1          0.0
Cross wind distance [m]:        0.0             1000.0          2000.0
Height difference [m]:          20.0            20.0            20.0

------- Hornsrev1 -------
Wind direction:                 0 deg           30 deg          90 deg
Down wind distance [m]:         2000.0          1732.1          0.0
Cross wind distance [m]:        0.0             1000.0          2000.0
Height difference [m]:          20.0            20.0            20.0

------- ParqueFicticio -------
Wind direction:                 0 deg           30 deg          90 deg
Down wind distance [m]:         nan             nan             nan
Cross wind distance [m]:        0.0             1000.0          2000.0
Height difference [m]:          20.0            20.0            20.0

------- Hornsrev1WRF -------
Wind direction:                 0 deg           30 deg          90 deg
Down wind distance [m]:         2177.2          2997.4          0.0
Cross wind distance [m]:        -857.1          -1441.5         2000.0
Height difference [m]:          20.0            20.0            20.0

Streamline distance

The StreamlineDistance takes a vector field or function as input and computes the streamline from e.g. a wind turbine and downstream. The downwind distance is then calculated as the distance along the streamline until the point where the streamline intersects with the rotor plane of the downstream wind turbine, and the cross-wind and vertical distance as the distance from the intersection point to the center of the rotor.

[27]:

site = sites['Hornsrev1WRF']
t = 1281
ds = site.ds.sel(time=t).interp(h=119)
sx, sy = slice(5,-6), slice(2,-4)
X,Y = np.meshgrid(ds.x.values[sx], ds.y.values[sy])
theta = np.deg2rad(270-ds.WD[sx,sy])
WS = ds.WS[sx,sy]
Vx,Vy = np.cos(theta)*WS, np.sin(theta)*WS
ds.WS[sx,sy].plot.contourf(levels=50, x='x', cmap='Blues_r')
plt.quiver(X, Y, Vx.T, Vy.T, alpha=1, width=.003)

wt_x, wt_y = site.initial_position[[0]].T

stream_lines = site.distance.vectorField.stream_lines(wd=np.full_like(wt_x,0), start_points=np.array([wt_x, wt_y, np.full_like(wt_x,wf_h)]).T,
                                                      time=np.array([t]), dw_stop=np.full_like(wt_x,10000))
for sl in stream_lines:
    plt.plot(sl[:, 0], sl[:, 1], 'r', lw=1)
wts.plot([wt_x,wt_x],[wt_y,wt_y-2000])
setup_plot(grid=False)

Wind resource distribution plots

The Site object has a few plot function to visualize its properties, mainly the wind resource given by the wind rose and the probability functions.

[28]:

site = sites['Hornsrev1']

Plotting wind rose.

[29]:

_ = site.plot_wd_distribution(n_wd=12, ws_bins=[0,5,10,15,20,25])

Plotting probability density function for the four sectors studied.

[30]:

_ = site.plot_ws_distribution(wd=[0,90,180,270])

Plotting probablity density function for the four sector studied \(\pm\) 45 degrees.

If include_wd_distribution=true, the wind speed probability distributions are multiplied by the wind direction probability.

The sector size is set to 360 / len(wd). This only makes sense if the wd array is evenly distributed

[31]:

_ = site.plot_ws_distribution(wd=[0,90,180,270], include_wd_distribution=True)