WAsP Flow Model
WAsP includes flow models for horizontal and vertical extrapolation of wind climates. The models includes effects on the flow from terrain inhomogeneity and from the structure of the atmospheric surface layer.
Illustration of the generalization (up arrow) and downscaling (down arrow) which is the key idea of extrapolation of wind climates in WAsP
- Core to the flow extrapolation in WAsP is the two steps illustrated in the figure above:
Generalization
Downscaling
Generalization is the process of “removing” local terrain and surface layer effects from the observed wind climate and instead “generalizing” the wind climate to a number of “standard” conditions (of surface roughness and height above ground).
Downscaling is the reverse process where local terrain and surface layer effects are “added” to the generalized wind climate to obtain a new predicted wind climate at the desired location.
A generalized wind climate tries to be representative of a larger area, an observed and prediced wind climate is representative only of a specific location.
To illustrate the flow model, we will import pywasp, windkit, and some other commonly used libraries.
In [1]: import numpy as np
In [2]: import xarray as xr
In [3]: import matplotlib.pyplot as plt
In [4]: import windkit as wk
In [5]: import pywasp as pw
We will also use terrain and wind data from the Santa Luzia site.
In [6]: bwc = wk.read_bwc("./source/tutorials/data/SerraSantaLuzia.omwc", crs="EPSG:4326")
In [7]: bwc = wk.spatial.reproject(bwc, to_crs="EPSG:32629")
In [8]: elev_map = wk.read_vector_map(
...: "source/tutorials/data/SerraSantaLuzia.map",
...: map_type="elevation",
...: crs="EPSG:32629",
...: )
...:
In [9]: lc_map, lc_tbl = wk.read_vector_map(
...: "source/tutorials/data/SerraSantaLuzia.map",
...: map_type="roughness",
...: crs="EPSG:32629",
...: )
...:
In [10]: topo_map = pw.wasp.TopographyMap(
....: elev_map,
....: lc_map,
....: lc_tbl
....: )
....:
To get a sense of the topography at the site, we will plot the elevation and roughness maps. It easiest to visualize the roughness data as a raster map, because the vector map version is roughness-change lines, which can be hard to detect patterns from. So, we will convert to raster data first, and then plot it. First the elevation:
In [11]: elev_raster = pw.io.vectormap_to_rastermap(elev_map, res=200)
In [12]: elev_raster.plot()
Second, the surface roughness:
In [13]: z0_map = wk._vectormap_helpers._landcover_to_z0(lc_map, lc_tbl)
In [14]: z0_raster = pw.io.vectormap_to_rastermap(z0_map, res=200)
In [15]: np.log1p(z0_raster).plot() # log(z0+1) to better see patterns
Note
The function to convert a land cover map to a roughness map wk._vectormap_helpers._landcover_to_z0 is considered experiemental for now, and is not currently part of the offical WindKit API.
From the maps we see flat land with low roughness along the western edge covered by ocean. Inland, to the northwest, and to the northeast we see two hills with lower roughness on the top. In the south, low surface roughness and flat elevation reveals a river wedged into the land, stretching from west to east.
Terrain effects
Now that we have created a TopographyMap and inspected the topography at the site,
let us start modeling terrain effects using the WAsP model.
First, we will define output locations by using windkit.spatial.create_dataset().
In this case, we define a 41x41x1 cuboid-shaped output grid (1681 total points).
Second, we call TopographyMap.get_site_effects() to get the terrain effects at the
output locations. We also have to tell how many sectors we want to use (we use 12 in this case).
In [16]: dx, dy = 300, 300
In [17]: nx, ny = 51, 51
In [18]: x0, y0 = 510500, 4613500
In [19]: steps = np.arange(nx)*dx
In [20]: x = x0 + steps
In [21]: y = y0 + steps
In [22]: xmin, xmax = x.min(), x.max()
In [23]: ymin, ymax = y.min(), y.max()
In [24]: output_locs = wk.create_dataset(
....: west_east=x,
....: south_north=y,
....: height=[50],
....: crs="EPSG:32629",
....: struct="cuboid",
....: )
....:
In [25]: site_effects = topo_map.get_site_effects(output_locs, nsecs=12)
In [26]: print(site_effects)
<xarray.Dataset>
Dimensions: (height: 1, south_north: 51, west_east: 51, sector: 12)
Coordinates:
* height (height) int64 50
* south_north (south_north) int64 4613500 4613800 ... 4628200 4628500
* west_east (west_east) int64 510500 510800 ... 525200 525500
* sector (sector) float64 0.0 30.0 60.0 ... 270.0 300.0 330.0
crs int8 0
sector_ceil (sector) float64 15.0 45.0 75.0 ... 285.0 315.0 345.0
sector_floor (sector) float64 345.0 15.0 45.0 ... 255.0 285.0 315.0
Data variables: (12/14)
z0meso (sector, south_north, west_east) float32 0.001509 .....
slfmeso (sector, south_north, west_east) float32 0.2933 ... 1.0
displ (sector, south_north, west_east) float32 0.0 ... 0.0
user_def_speedups (sector, height, south_north, west_east) float32 1.0...
orographic_speedups (sector, height, south_north, west_east) float32 0.9...
obstacle_speedups (sector, height, south_north, west_east) float32 1.0...
... ...
orographic_turnings (sector, height, south_north, west_east) float32 -0....
obstacle_turnings (sector, height, south_north, west_east) float32 0.0...
roughness_turnings (sector, height, south_north, west_east) float32 0.0...
site_elev (south_north, west_east) float32 0.1 0.1 ... 781.0
rix (south_north, west_east) float32 0.0 0.0 ... 0.1406
dirrix (sector, south_north, west_east) float32 0.0 ... 0.0
Attributes:
Conventions: CF-1.8
history: 2023-07-26T08:48:15:\twindkit==0.6.3\t wk.create_datase...
Package name: windkit
Package version: 0.6.3
Creation date: 2023-07-26T08:48:17
Object type: Topographic effects
author: Neil Davis
author_email: neda@dtu.dk
institution: DTU Wind
title: WAsP site effects
Note
pywasp.wasp.TopographyMap.get_site_effects() uses 1x runs in your PyWAsP subscription
The site effects are a xarray.Dataset with many variables (see e.g. Data Structures for a full list),
including speed-up effects and turnings from the orographic model and the roughness model.
It also holds variables “z0meso”, the effective upstream sector-wise roughness at each location,
and “slfmeso”, the effective sector-wise sea-to-land fraction (1 being land, 0 being ocean).
The variable “rix”, is the ruggedness index, which is a number between 0 and 1, signifying the
orographic complexity of the site at each point (0 being not very complex).
We can plot the terrain effects to get a better idea of their characteristics at the site. First, we will plot the orographic speed-up for one wind direction (210 degrees in this case). To show the wind direction in the map, a useful plotting function is defined for plotting the direction as an arrow on the map.
In [27]: def plot_arrow(x=0.03, y=0.03, ws=0.08, wd=210, ax=None):
....: """ helper function to plot the wind direction as an arrow """
....: if ax is None:
....: ax = plt.gca()
....: ax.arrow(
....: x,
....: y,
....: *wk.wind_vectors(0.08, 210),
....: transform=ax.transAxes,
....: head_width=0.04,
....: head_length=0.1,
....: alpha=0.8,
....: fc='k',
....: ec='k'
....: )
....: return ax
....:
In [28]: fig, ax = plt.subplots()
In [29]: ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax));
In [30]: (
....: site_effects["orographic_speedups"]
....: .sel(sector=210, height=50) # choose sector and height to plot
....: .plot.contourf(vmin=0.5, vmax=1.5, cmap="RdBu_r", ax=ax, add_colorbar=True)
....: );
....:
In [31]: plot_arrow(x=0.03, y=0.03, ws=0.08, wd=210, ax=ax);
In [32]: elev_raster.plot.contour(levels=16, cmap="gray", ax=ax, add_colorbar=False);
As expected, we see speed-ups on up-slope regions and slow-downs on down-slope regions. Next, we can take a look at the orographic turning:
In [33]: fig, ax = plt.subplots(1, 1)
In [34]: ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax));
In [35]: (
....: site_effects["orographic_turnings"]
....: .sel(sector=210, height=50)
....: .plot.contourf(center=0, ax=ax, add_colorbar=True)
....: );
....:
In [36]: plot_arrow(x=0.03, y=0.03, ws=0.08, wd=210, ax=ax);
In [37]: elev_raster.plot.contour(levels=16, cmap="gray", ax=ax, add_colorbar=False);
Now, let us plot the roughness speed-up:
In [38]: fig, ax = plt.subplots(1, 1)
In [39]: ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax));
In [40]: (
....: site_effects["roughness_speedups"]
....: .sel(sector=210, height=50)
....: .plot.contourf(vmin=0.8, vmax=1.2, cmap="RdBu_r", ax=ax, add_colorbar=True)
....: );
....:
In [41]: plot_arrow(x=0.03, y=0.03, ws=0.08, wd=210, ax=ax);
In [42]: np.log1p(z0_raster).plot.contour(levels=16, cmap="gray", ax=ax, add_colorbar=False);
We can see large areas of slow-down of the flow as it transitions from low roughness over the sea to higher surface roughness over land.
Generalization
To generalize a observed/binned wind climate, generalize() is used:
In [43]: gwc = pw.wasp.generalize(bwc, topo_map)
In [44]: print(gwc)
<xarray.Dataset>
Dimensions: (point: 1, sector: 12, gen_roughness: 5, gen_height: 5)
Coordinates:
height (point) float64 25.3
south_north (point) float64 4.621e+06
west_east (point) float64 5.147e+05
crs int8 0
sector_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
* sector (sector) float64 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.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 5.715 ....
k (sector, gen_height, gen_roughness, point) float32 1.721 ....
wdfreq (sector, gen_height, gen_roughness, point) float32 0.06491...
Attributes:
Conventions: CF-1.8
history: 2023-07-26T08:48:13:\twindkit==0.6.3\tcreate_dataset(we...
wasp_header: SerraSantaluzia
Package name: windkit
Package version: 0.6.3
Creation date: 2023-07-26T08:48:13
Object type: Binned Wind Climate
author: Neil Davis
author_email: neda@dtu.dk
institution: DTU Wind
title: Generalized wind climate
The result is a Generalized Wind Climate, which contains weibull parameters
for a number of surface roughnesses and heights
(prefixed with _gen to signify that these are conditions in the “generalization space”).
Note
The gen_heights and gen_roughnesses classes should be set to capture the characteristics at the site
as well as possible. It is recommended to select heights that cover all heights one wishes to predicted, and to include
specific heights of interest as standalone classes. The same goes for the surface roughness classes.
Note
pywasp.wasp.generalize() uses 1x runs in your PyWAsP subscription
Downscaling
To downscale a GWC to the desired output locations, using pre-calculated site_effects, downscale() is used:
In [45]: pwc = pw.wasp.get_wasp_down(gwc, site_effects, genwc_interp="nearest")
In [46]: print(pwc)
<xarray.Dataset>
Dimensions: (height: 1, south_north: 51, west_east: 51, sector: 12)
Coordinates:
* height (height) int64 50
* south_north (south_north) int64 4613500 4613800 ... 4628200 4628500
* west_east (west_east) int64 510500 510800 511100 ... 525200 525500
* sector (sector) float64 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
crs int8 0
sector_ceil (sector) float64 15.0 45.0 75.0 105.0 ... 285.0 315.0 345.0
sector_floor (sector) float64 345.0 15.0 45.0 75.0 ... 255.0 285.0 315.0
Data variables:
A (sector, height, south_north, west_east) float32 6.606 ... ...
k (sector, height, south_north, west_east) float32 1.998 ... ...
wdfreq (sector, height, south_north, west_east) float32 0.06397 .....
site_elev (south_north, west_east) float32 0.1 0.1 0.1 ... 773.6 781.0
Attributes:
Conventions: CF-1.8
history: 2023-07-26T08:48:15:\twindkit==0.6.3\t wk.create_datase...
Package name: windkit
Package version: 0.6.3
Creation date: 2023-07-26T08:48:19
Object type: Weibull Wind Climate
author: Neil Davis
author_email: neda@dtu.dk
institution: DTU Wind
title: WAsP site effects
The prediced wind climate is a weibull wind climate, containing sector-wise weibull parameters for the desired locations.
To downscale without pre-calculating the site effects, pw.wasp.downscale()
can be used instead:
In [47]: pwc = pw.wasp.downscale(gwc, topo_map, output_locs, genwc_interp="nearest")
Note
pywasp.wasp.downscale() uses 1x runs in your PyWAsP subscription
One-step extrapolation
For direct horizontal and vertical extrapolation of wind climates,
wherein the generalization and downscaling is done in a single call, the
generalize_and_downscale() can be used.
In [48]: pwc = pw.wasp.generalize_and_downscale(output_locs, bwc, topo_map)
Note
pywasp.wasp.generalize_and_downscale() uses 2x runs in your PyWAsP subscription
WAsP Model Configuration
When running WAsP through PyWAsP, there are a number of parameters that can be set through the
conf argument to the functions: TopographyMap.get_site_effects(),
generalize(), get_wasp_down(), and downscale().
The conf input is a Config object contains the model parameters, and is set
by default to correspond to the WAsP GUI version 12.8, using the latest version of WAsP
and a set of pre-defined parameters for the terrain and climate models.
To use a previous version of WAsP or to fine tune the model parameters, we can create an instance
of the Config class, which then can be passed as an argument to the functions:
In [49]: conf = pw.wasp.Config("WAsP_12.7") # Use a parameter-set from a previous WAsP version.
In [50]: pwc = pw.wasp.generalize_and_downscale(
....: output_locs,
....: bwc,
....: topo_map,
....: generalize_kwargs={"conf": conf},
....: downscale_kwargs={"conf": conf},
....: )
....:
The Config object can be initialized with a str, as above, which uses
parameters from a specific version of WAsP:
“WAsP_11.4”: WAsP 11 profile model using c_g=1.65 and powerlaw estimation for reversal height. Classic terrain analysis.
“WAsP_12.6”: WAsP 12 profile model with c_g solved by iteration and dA/dmu set such that it corresponds with the old profile model. Classic terarin analysis
“WAsP_12.7”: WAsP 12 profile model, as version 12.6 but including baroclinicity. Spider terrain analysis.
“WAsP_12.8”: WAsP 12 profile model, including sectorwise stability and baroclinicity as described in mesoclimate class (Profile model 3); terrain analysis 1
Initializing the object with no parameters, means the default, “WAsP_12.8”, is used. Once initialized, the object can also be further modified by manually chaging specific parameters.
The Config class sets parameters for the terrain and the climate WAsP models. An
instance of Config has two properties, climate and terrain, which are
instances of Climate and Terrain.
The detail of the set values can be check as follows:
config = pywasp.wasp.Config()
config.climate
"""
Index Parameter Current Default Description
value value
===== =============== ============ =========== ==================================================
10 a0 1.800000 1.8 A0 parameter in geostrophic drag law under neutral conditions
11 b0 4.500000 4.5 B0 parameter in neutral geostrophic drag law
56 rms_hfx_land 100.000000 100 Root mean square of heat flux over land
57 rms_hfx_water 30.000000 30 Rms Heatflux over water
58 off_hfx_land -40.000000 -40 Offset (mean) heatflux over land
59 off_hfx_water -8.000000 -8 Offset (mean) heatflux over water
60 stab_lim_land 1.000000 1 Factor on geostrophic height for stab envelope limiting
77 stab_lim_water 1.000000 1 Factor on geostrophic height (sea) for stability envelope limiting"""
config.terrain
"""
Index Parameter Current Default Description
value value
===== =============== ============ =========== ==================================================
31 decay_length 10000.000000 10000 Decay-length for importance of roughness relative to distance from center
67 max_no_rgh 10.000000 10 Max number of roughness changes per sector."""
The configuration values stored as key - value pairs and can be accessed or set individually by index
config.climate[10] = 3.6
config.climate
"""
Index Parameter Current Default Description
value value
===== =============== ============ =========== ==================================================
10 a0 3.600000 1.8 A0 parameter in geostrophic drag law under neutral conditions
11 b0 4.500000 4.5 B0 parameter in neutral geostrophic drag law
56 rms_hfx_land 100.000000 100 Root mean square of heat flux over land
57 rms_hfx_water 30.000000 30 Rms Heatflux over water
58 off_hfx_land -40.000000 -40 Offset (mean) heatflux over land
59 off_hfx_water -8.000000 -8 Offset (mean) heatflux over water
60 stab_lim_land 1.000000 1 Factor on geostrophic height for stab envelope limiting
77 stab_lim_water 1.000000 1 Factor on geostrophic height (sea) for stability envelope limiting"""
It is also possible to set a combination of values to define a climate profile model and a terrain analysis.
climate.set_profile_model(0)
terrain.set_terrain_analysis(1)
For further details of the above, check the API Reference Config.