[1]:
%load_ext autoreload
%autoreload 2
[2]:
import warnings
warnings.filterwarnings('ignore') # We will ignore warnings to avoid cluttering the notebook
import numpy as np
import windkit as wk
import pywasp as pw
from tutorial_helpers import *
Evaluating roughness maps at Østerild
In this tutorial you will use several new features related to the treatment of roughness maps; these functionalities are not available in the WAsP GUI, but have been developed using the pywasp
interface. Data from the Danish national test center for large wind turbines at Østerild will be used: there are two meteorological masts located to the north and south of a row of turbines (red and black points in figure below, respectively). The masts have measurements up to 244 m height, and are
instrumented with both cup and sonic anemometers. The wind data used in this exercise are quality-controlled: periods where wakes and icing influenced the measurements have been removed. In addition, sectors that were influenced by the wakes of nearby turbines have been filtered out: this causes some sectors to have no observations (tripping a warning when reading the observed wind climates).
Østerild wind data from Dec. 13th of 2016 to May 24th of 2018 (71802 10-min measurements) are used, with the following data and processing details: * data available from both northern and southern masts (red triangles in figure below); * icing episodes removed by selecting data with standard deviation \sigma_{U}>10^{-5}m/s; * flow distortion correction via “Flow distortion on boom-mounted cup
anemometers” (P. Lindelöw-Marsden et al.); * wind directions linearly interpolated between 40 and 244 m; * data from sectors with wind turbine wakes were removed; * 12-sector histogram files were already generated using pywasp
.
Reading Østerild wind data
First we will read mast data from the two masts; these have been stored as multiple tab files, and are read using the “helper” function read_tabfiles_oesterild
made for this tutorial.
[3]:
tabs = read_tabfiles_oesterild()
print(tabs)
<xarray.Dataset>
Dimensions: (point: 16, sector: 12, wsbin: 39)
Coordinates:
height (point) float64 10.0 40.0 70.0 106.0 ... 178.0 210.0 244.0
crs int8 0
* wsbin (wsbin) float64 0.5 1.5 2.5 3.5 4.5 ... 35.5 36.5 37.5 38.5
wsceil (wsbin) float64 1.0 2.0 3.0 4.0 5.0 ... 36.0 37.0 38.0 39.0
wsfloor (wsbin) float64 0.0 1.0 2.0 3.0 4.0 ... 35.0 36.0 37.0 38.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
west_east (point) float64 4.928e+05 4.928e+05 ... 4.928e+05 4.928e+05
south_north (point) float64 6.323e+06 6.323e+06 ... 6.327e+06 6.327e+06
Dimensions without coordinates: point
Data variables:
wdfreq (sector, point) float64 0.0 0.0 0.0 ... 0.0407 0.04364 0.04703
wsfreq (wsbin, sector, point) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
Conventions: CF-1.8
history: 2024-06-03T15:08:26+00:00:\twindkit==0.7.1.dev53+gcced3...
wasp_header: Mast 2 10 m from timeseries
Package name: windkit
Package version: 0.7.1.dev53+gcced321
Creation date: 2024-06-03T15:08:26+00:00
Object type: Binned Wind Climate
author: Bjarke Tobias Olsen
author_email: btol@dtu.dk
institution: DTU Wind Energy
header: Binned wind climates (bwc) from Østerild from the North...
We note in the Dimensions
of the xarray.Dataset
that was read as tabs
shows 16 point
s; the first 8 points correspond to the South mast, while the last 8 points correspond to the Northern mast.
These are binned wind climate files, whose format was explained earlier in the resource grid tutorial. I.e., there are no mean speeds, but rather probabilities (frequencies of occurence) per wind speed bin (wsbin
) and directional sector
stored in wsfreq
(with directional probabilities over all speeds in wdfreq
).
Preparing map files
Traditionally WAsP has only been able to use vector maps for roughness length. In pywasp
it is very easy to work with both raster (grid) and vector (contour) maps, with simple functions to convert between formats.
In this tutorial we’ll load raster map files from a number of sources, in the following order: * Danish Digital Terrain-elevation Map (DTM). For the description of the terrain height we use the digital elevation model from Denmark’s Geodatastyrelsen, at a resolution of 100 m (the full dataset was obtained via airplane-based lidar scans, and much finer resolutions can also be downloaded). * CORINE landcover. The European Union’s CORINE initiative (‘COoRdination of INformation on the Environment’) is popular for wind resource assessments in Europe, due to its relatively fine resolution of 100 m (e.g. it is used for land cover in the New European Wind Atlas). * Sentinel satellite data. Through the INNOWIND project, which aims to use new information from satellites to create more accurate wind resource assessments, a new land cover product has been created. It has 20 m resolution and provides not just landcover classes, but also tree height and leaf area index in forested areas.
[4]:
elev_ras = wk.read_raster_map(os.path.join('data/maps','Oesterild_DTM_100m.grd'), map_type='elevation',crs=32632)
After reading the raster map of terrain elevation above, we now convert it to a vector map using the pywasp.io.rastermap_to_vectormap
function previously outlined and used in the tutorial on resource grids.
[5]:
elev_map = pw.rastermap_to_vectormap(elev_ras)
To load a raster map where each grid-cell (pixel) represents a roughness length, we can specify that we are loading a map of type 'roughness'
within the read_rastermap
function. Internally pywasp
always employs a vectormap where each contour has two integer values representing the land use classes on each side of the contour, along with a roughness look-up table that provides the z_0 associated with each class. This allows more detailed formulation of landcover, as we will
see later in this tutorial; it also allows one to easily modify the roughness lengths in a project without operating on the map (as discussed in the earlier tutorial on resource grids). > pywasp note: it is a neccessary step for roughness data in raster format to be converted to vector representation (contours of roughness change) for the flow-perturbation calculations. pywasp
can work with raster maps directly, but it is recommended to convert them to vectors first.
[6]:
corine_z0_ras, corine_lctable = wk.read_raster_map(os.path.join('data/maps','corine.grd'), map_type='roughness',crs=32632)
corine_z0_map, corine_lctable = pw.rastermap_to_vectormap(corine_z0_ras, lctable=corine_lctable)
Using landcover (LC) maps
Above we loaded raster mapfiles with commonly-used elevation and roughness length data, and transformed them to vector format; along with reading the roughness lengths we converted them to landcover classes, also getting a corresponding lookup table.
But what if we have LC maps, and not z_0? This is actually the case for the Sentinel data.
Below we see how to load a raster map of landcover and corresponding landcover tables, using the pw.LandCoverTable.read_json
and pw.read_rastermap
methods (the cells in the map represent landcover classes, which correspond to the z_0 in the lookup table).
This is done for both the CORINE data and the Sentinel data, along with converting each map to vector format:
[7]:
corine_lctable = pw.LandCoverTable.read_json(os.path.join('data/maps','corine_DTU.json'))
corine_ras = wk.read_raster_map(os.path.join('data/maps','corine_lu.grd'), map_type='landcover',crs=32632)
corine_map, _ = pw.rastermap_to_vectormap(corine_ras, lctable=corine_lctable)
sentinel_lctable = pw.LandCoverTable.read_json(os.path.join('data/maps','scadis_sentinel.json'))
sentinel_ras = wk.read_raster_map(os.path.join('data/maps','sentinel.tif'), map_type='landcover', crs=32632)
sentinel_map, _ = pw.rastermap_to_vectormap(sentinel_ras, lctable=sentinel_lctable)
Roughness and landcover around masts
Let’s explore the maps that have been loaded, using the ‘plot()’ method of xarray
(recalling that 'corine_ras'
is an xarray DataArray, and find the approximate roughness and landcover classes around the masts. We also use red dots to show the locations of the masts, which were contained in the 'tabs'
variable created by reading the wind data at the beginning of this tutorial.
[8]:
x_mastS = tabs['west_east'][0]
x_mastN = tabs['west_east'][15]
y_mastS = tabs['south_north'][0]
y_mastN = tabs['south_north'][15]
[9]:
corine_ras.plot(figsize=(10,7))
plt.plot(x_mastN, y_mastN, 'or', ms=5)
plt.plot(x_mastS, y_mastS, 'or', ms=5)
plt.show()
we can zoom in more (limiting the plot area using xlim
and ylim
) and make the color bar discrete (noting there are 44 LC types in 'corine_ras'
and changing the color-map), to better see the LC’s around the masts:
[10]:
rzoom = 2500
my_cmap=plt.get_cmap('terrain', 44)
corine_ras.plot(figsize=(10,7), cmap=my_cmap)
plt.plot(x_mastN, y_mastN, 'or', ms=5)
plt.plot(x_mastS, y_mastS, 'or', ms=5)
plt.xlim([x_mastS - rzoom, x_mastS + rzoom])
plt.ylim([y_mastS - rzoom, y_mastN + rzoom])
plt.show()
From the CORINE map we can find out that much of the LandCoverID
around the mast locations is 24, which according to the CORINE look-up table corresponds to “Coniferous forest” or z0 = 1.2
. On the other hand, from the Sentinel map LandCoverID
for the mast locations is 1, which corresponds to “non-forest (cropland grassland other)” or z0 = 0.1
.
[11]:
sentinel_ras.plot(figsize=(10,7), cmap=plt.get_cmap('rainbow', 44))
plt.plot(x_mastN, y_mastN, 'or', ms=10)
plt.plot(x_mastS, y_mastS, 'or', ms=10)
plt.xlim([x_mastS - rzoom, x_mastS + rzoom])
plt.ylim([y_mastS - rzoom, y_mastN + rzoom])
plt.show()
As seen above, the Sentinel data in this case appears to include much more detail than CORINE; e.g. we can see line-like features where the forests have been cut. (Note: water is class 5 in the Sentinel lookup-table, where the classes are not ordered according to z_0. We have picked a different color table to accomodate this and make it blue. You can check the LC(z_0) by simply looking at sentinel_lctable
.)
Now we combine the vector elevation map obtained earlier with the CORINE and Sentinel landcover vector maps and CORINE z0 vector map, and make a TopographyMap
for each pair. Such combination of elevation and roughness/landcover maps is done, because they are both needed for the flow modelling in WAsP/pywasp
. Here we store them in a dictionary ('allmaps'
) for further processing:
[12]:
allmaps = {}
allmaps['Sentinel landcover'] = pw.wasp.TopographyMap(elev_map, sentinel_map, sentinel_lctable)
allmaps['CORINE landcover']= pw.wasp.TopographyMap(elev_map, corine_map, corine_lctable)
allmaps['CORINE z0']= pw.wasp.TopographyMap(elev_map, corine_z0_map, corine_lctable)
[13]:
corine_ras.plot()
plt.show()
Exploring roughness around a point: the roughness rose
WAsP’s internal boundary layer (IBL) model computes speedup effects due to (a limited number of) the most significant roughness changes in each directional sector; this was discussed in the earlier tutorial on (py)WAsP’s response to roughness changes. A roughness rose can be plotted, which displays the sector-wise z_0 as a function of distance.
But to use WAsP’s roughness-change model, and for displaying z_0-roses, we must set up pywasp
’s configuration. This was also briefly covered/used in the resource_grid tutorial.
Here we use the default WAsP configuration, but additionally employ the new landcover roughness model; the latter adds the ability to deal with displacement heights, i.e. treating flow that is “displaced” upward over forests based on mean tree height and leaf-area density (as found e.g. via Sentinel data).
pywasp note: WAsP has about 60 parameters and constants that can be changed to control the behaviour of its various submodels, most importly those involving terrain-induced perturbations, vertical extrapolation, and getting/applying generalized wind statistics.
To get an overview of the settings of the various parameters, you can print the conf.terrain
and conf.climate
configurations. In practice we occasionally adjust only a small fraction of parameters from their default values, just as in WAsP, though further research and use of pywasp
may change this; also, only some are actually visible to end-users. In-depth docummentation of them (beyond the European Wind Atlas, 1989) is currently under way, though one can use getpardesc(pno)
for
parameter number pno
to get some description (done below).
[14]:
conf = set_config()
conf.terrain
[14]:
Index Parameter Current Default Description
value value
===== ========================== ============ =========== ==================================================
31 decay_length 10000.000000 10000 Decay-length for importance of roughness relative to distance from center
40 upper_kink_ibl_profile 0.300000 0.3 Factor for upper kink in roughness change wind profile in the ibl
41 lower_kink_ibl_profile 0.090000 0.09 Factor for lower kink in roughness change wind profile in the ibl
67 max_no_roughness_changes 10.000000 10 Max number of roughness changes/sector. See \cite{Floors2018b} for description and sensitivity. This array is dimensioned in the variable mrchs in dimpar.f90.
68 max_rms_error 0.300000 0.3 Maximum root-mean square (rms) error in log(roughness) analysis of roughness changes for a given sector. Setting a lower error threshold will allow a higher number of roughness changes to be taken into consideration, for each sector. Choosing 0 causes this parameter to become inactive, with the maximum number of changes then set by parameter~67.
69 use_displacement_height 1.000000 1 add displacements to orographic grid if = 1, and reduce displacement slope near origin to parameter~70
82 roughness_analysis 1.000000 1 rou analysis version: 0=old, =1 use spiderweb with old rou maps
97 water_roughness_length 0.000200 0.0002 ($2\times10^{-4}$) water roughness length used in program
[15]:
conf.terrain.get_desc(67)
================================================================================
max_no_roughness_changes
================================================================================
standard_name : max_no_of_roughness_changes
short_name : max_no_roughness_changes
long_name : max_no_of_roughness_changes_per_sector
GUID : MAXRC
Rvea : Rvea0287
model : IBL/roughness-change model
index : 67
valid_min : 1
valid_max : 20
valid_range : [1, 20]
default_value : 10
definition:
Max number of roughness changes/sector. See \\cite{Floors2018b} for
description and sensitivity. This array is dimensioned in the variable mrchs
in dimpar.f90.
================================================================================
[16]:
conf.terrain.get_desc(31)
================================================================================
decay_length
================================================================================
standard_name : roughness_decay_length
short_name : decay_length
long_name : Decay-length for roughness area size
GUID : RDECAYL
Rvea : Rvea0287
model : WAsP IBZ flow modelling
index : 31
units : m
valid_min : 1000
valid_max : 1000000
valid_range : [1000, 1000000]
default_value : 10000
definition:
Decay-length for importance of roughness relative to distance from center
reference : European Wind Atlas---\cite{Troen1989}, around equation~8.24
notes:
("Used to get 'meso'-roughness: determines relative weight to roughnesses "
'depending on distance from wind calculation point (center).
================================================================================
To see the roughness roses of all maps at a certain location we can do the following: * loop over all maps stored in the dictionary 'allmaps'
* use the sel
method from xarray to select the first point in tabs
(South mast, height = 10 m via point=0
) * [also possible: select the site factors using the ``to_point`` method] * plot a roughness rose for each map
NOTE: we arbitrarily select from 10 m height, because the roughness rose does not depend on height
[17]:
tabs0 = tabs.sel(point=0)
south_mast_loc = wk.create_dataset(
tabs0.west_east,
tabs0.south_north,
tabs0.height,
wk.spatial.get_crs(tabs0)
)
for key, topo_map in allmaps.items():
print(key)
rou_rose, _ = topo_map.get_rou_rose(south_mast_loc, nsecs=12, conf=conf)
wk.plot.roughness_rose(rou_rose)
Sentinel landcover
CORINE landcover
CORINE z0
Comparison of observed and modelled wind, per roughness maps
We can now proceed with a comparison between the observed and modelled wind speed (and its distribution), using the WAsP model to vertically extrapolate the measurements from a certain height to another; the wind profile and thus this extrapolation depends upon z_0.
To start we first need to generalize the observed wind climate, after “cleaning” local terrain effects from the wind observations; then after downscaling the generalized wind climate, we re-apply local effects for another height or location. The processes are performed by the generalize
and downscale
methods, as shown in the earlier tutorial on resource maps.
Below we select observations from a certain height and mast and create a generalized wind climate (gwc) from them: * the first argument is input binned wind climate (remember first 8 points correspond to the South mast, remaining to the North mast) * the second argument is the mapfile that we want to use to generalize the observations * the third argument is the WAsP model configuration
[18]:
bwc = tabs.sel(point=3, drop=False) # Selecting wind data from the south mast at 106 m above ground level
gwc = pw.wasp.generalize(
bwc,
allmaps['Sentinel landcover'],
conf
)
[19]:
print(gwc)
<xarray.Dataset>
Dimensions: (point: 1, sector: 12, gen_roughness: 5, gen_height: 5)
Coordinates:
height (point) float64 106.0
south_north (point) float64 6.323e+06
west_east (point) float64 4.928e+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.03 0.1 0.4 1.5
* gen_height (gen_height) int64 10 25 50 100 250
Dimensions without coordinates: point
Data variables:
A (sector, gen_height, gen_roughness, point) float32 6.858 ....
k (sector, gen_height, gen_roughness, point) float32 1.998 ....
wdfreq (sector, gen_height, gen_roughness, point) float32 0.00775...
site_elev (point) float32 10.15
Attributes:
Conventions: CF-1.8
history: 2024-06-03T15:08:26+00:00:\twindkit==0.7.1.dev53+gcced3...
wasp_header: Mast 2 10 m from timeseries
Package name: windkit
Package version: 0.7.1.dev53+gcced321
Creation date: 2024-06-03T15:08:40+00:00
Object type: Geostrophic Wind Climate
author: Bjarke Tobias Olsen
author_email: btol@dtu.dk
institution: DTU Wind Energy
header: Binned wind climates (bwc) from Østerild from the North...
title: Generalized wind climate
Now that we have a generalized wind climate, we can reapply this wind climate for any nearby location. In this case we will create a horizontal transect of locations along a line connecting the meteorological masts, from 10 km south to 10 km north, at 106 m above ground level. The coordinates of the mast are taken from bwc
. The output of the downscale function is a Weibull wind climate ('wwc'
), i.e. sector-wise Weibull distributions at our transect.
[20]:
loc_x = float(bwc.west_east)
loc_y = float(bwc.south_north)
half_ext = 10000.0 # half extent in x and y in m
points_y = np.arange(loc_y-10000,loc_y+10000,2000)
transect = wk.create_dataset(np.repeat(loc_x,len(points_y)),
points_y,
np.repeat(bwc.height.values,len(points_y)),
32632)
wwc = pw.wasp.downscale(
gwc,
allmaps['Sentinel landcover'],
transect,
conf,
interp_method='nearest'
)
[21]:
plot_transect_z0(tabs,transect,allmaps['CORINE landcover'],conf)
[22]:
plot_transect_displ(tabs,transect,allmaps['Sentinel landcover'],conf)
[23]:
plot_transect_ws(tabs,wwc)
[24]:
plot_transect_pd(tabs,wwc)
Vertical profile of power density
Instead of a horizontal transect, we can also easily plot a vertical profile of the power density. For this a we can use the function wk.plot_vertical_profile
. It takes the same arguments as the plot_transect_pd
function, i.e. the dataset of binned wind climates and a Weibull wind climate. To do this we will:
Use the
create_grid
function to create the output locations for the flow modellingPerform the downscaling step as described in the section “Validation of roughness maps”. We can use the same generalized wind climate as before.
Call the
wk.plot_vertical_profile
function using the objects that have been created in the previous two steps
Let’s do this first for the South mast:
[25]:
print(tabs)
<xarray.Dataset>
Dimensions: (point: 16, sector: 12, wsbin: 39)
Coordinates:
height (point) float64 10.0 40.0 70.0 106.0 ... 178.0 210.0 244.0
crs int8 0
* wsbin (wsbin) float64 0.5 1.5 2.5 3.5 4.5 ... 35.5 36.5 37.5 38.5
wsceil (wsbin) float64 1.0 2.0 3.0 4.0 5.0 ... 36.0 37.0 38.0 39.0
wsfloor (wsbin) float64 0.0 1.0 2.0 3.0 4.0 ... 35.0 36.0 37.0 38.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
west_east (point) float64 4.928e+05 4.928e+05 ... 4.928e+05 4.928e+05
south_north (point) float64 6.323e+06 6.323e+06 ... 6.327e+06 6.327e+06
Dimensions without coordinates: point
Data variables:
wdfreq (sector, point) float64 0.0 0.0 0.0 ... 0.04364 0.04703
wsfreq (wsbin, sector, point) float64 0.0 0.0 0.0 ... 0.0 0.0 0.0
air_density (point) float64 1.225 1.225 1.225 1.225 ... 1.225 1.225 1.225
wspd (point) float64 4.159 6.075 7.363 8.468 ... 10.22 10.62 10.93
power_density (point) float64 99.46 240.3 387.4 ... 1.148e+03 1.263e+03
Attributes:
Conventions: CF-1.8
history: 2024-06-03T15:08:26+00:00:\twindkit==0.7.1.dev53+gcced3...
wasp_header: Mast 2 10 m from timeseries
Package name: windkit
Package version: 0.7.1.dev53+gcced321
Creation date: 2024-06-03T15:08:44+00:00
Object type: Met fields
author: Bjarke Tobias Olsen
author_email: btol@dtu.dk
institution: DTU Wind Energy
header: Binned wind climates (bwc) from Østerild from the North...
[26]:
loc_x = tabs.west_east.values[0] # selecting the South mast x coordinate
loc_y = tabs.south_north.values[0] # selecting the South mast y coordinate
heights = tabs.height.values[:8] # selecting the South mast heights
bwc_smast= tabs.isel(point=range(0,8)) # South mast data
transect = wk.create_dataset(np.repeat(loc_x,len(heights)),
np.repeat(loc_y,len(heights)),
heights,
32632)
wwc_south = pw.wasp.downscale(
gwc, allmaps['Sentinel landcover'], transect, conf, interp_method='nearest' )
[27]:
wk.plot.vertical_profile(bwc_smast.power_density,wwc_south.power_density)