"""Lookup Tables for Extreme Wind Estimation
When calculating extreme winds using PyWAsP it is common to pre-calculate a lookup table
for a few wind speeds and then apply these to the Generalized Extreme Wind Climate to
get the localized extreme winds. This will reduce the calculation load, and provides
reasonable results. Here you will find functions to both create and apply the lookup
table.
"""
import xarray as xr
from windkit.spatial import (
create_dataset,
get_crs,
spatial_stack,
spatial_unstack,
interpolate_to_grid,
)
from windkit.metadata import _LINCOM_V50_ATTRS, update_var_attrs, update_history
from .fetch import calculate_fetchmap, roughness_to_landmask
from .fourier_space import FourierSpace
from .wind import create_wind
from .wind_level import WindLevel
from .wind_point import get_wind_points
from ..core import Rvea0293
DEFAULT_Z0 = 0.05
DEFAULT_HGT = 10.0
def _run_lincom(wind, out_ds, oro_grid, rgh_grid, mask_rast, mask_vect):
"""Calculate downscaled wind speed and direction using LINCOM for a single wind
Parameters
----------
wind : xarray.Dataset
LINCOM wind dataset
out_ds : xarray.Dataset
Dataset describing the spatial structure of the requested LINCOM results
oro_grid, rgh_grid : xarray.DataArray
Orographic and roughness grids as a DataArray
mask_rast : xarray.DataArray
RasterMap of water (0) and land (1) mask
mask_vect : VectorMap
VectorMap of water (0) and land (1) mask
Returns
-------
xarray.Dataset
Wind speed and direction from LINCOM results for a single wind
"""
# Calculate the fetch
fetch_grid = calculate_fetchmap(mask_vect, mask_rast, wind)
# Calculate the Fourier Space
fspc = FourierSpace.create_fourier_space(
oro_grid, rgh_grid[0], rgh_grid[1], fetch_grid, wind
)
# Calculate wind level
hgt_list = []
for hgt in out_ds.height:
wind_level = WindLevel.get_wind_level(fspc, hgt)
wp_ds = get_wind_points(wind_level, out_ds.sel(height=hgt))
wp_ds["height"] = [hgt]
hgt_list.append(wp_ds[["WS", "WD", "height"]])
return xr.merge(hgt_list)
[docs]
def create_lut(nsec, wind_speeds, out_ds, oro_grid, rgh_rast):
"""Run LINCOM to create lookup table for extreme wind estimation
Parameters
----------
nsec : int
Number of sectors to include in lookup table
wind_speeds : list or numpy.ndarray
List of wind speeds to include in lookup table
out_ds: xarray.Dataset
Dataset containing the spatial structure to calculate
oro_grid : xarray.DataArray
Raster of orographic heights as a DataArray
rgh_rast : RasterMap
Raster of roughness lengths as a RasterMap
Returns
-------
out : xarray.Dataset
Extreme wind lookup table
Notes
-----
This routine is setup to run LINCOM for creating lookup tables. It will run
for the provided wind speeds and sectors. Outputing only the downscaled wind
speed and direction.
"""
wind = create_wind(
"generalized",
1.0,
DEFAULT_HGT,
DEFAULT_Z0,
0.0,
0.0,
get_crs(oro_grid),
nsec=nsec,
)
# Convert roughness RasterMap to LINCOM grid, and calculate land mask
lc_grid, lctable = rgh_rast
mask_rast, mask_vect = roughness_to_landmask(lc_grid, lctable)
# Run LINCOM for each wind speed and sector
rslt_list = []
for wspd in wind_speeds:
in_wind = wind.copy()
in_wind["wind_speed"] = in_wind.wind_speed * wspd
sec_list = []
for sec in in_wind.sector:
sec_wind = in_wind.sel(sector=sec)
lc_rslt = _run_lincom(
sec_wind, out_ds, oro_grid, rgh_rast, mask_rast, mask_vect
)
# Add sector dimension before appending to list
sec_list.append(lc_rslt.expand_dims({"sector": [sec]}))
# Aggregate sectors and add input wind speed as dimension
sec_ds = xr.merge(sec_list)
rslt_list.append(sec_ds.expand_dims({"gen_wspd": [wspd]}))
out = xr.merge(rslt_list)
out.attrs = out_ds.attrs
return update_history(out)
def _lut2d(lut, gewc):
"""Python wrapper for Rvea0293.calc_points_lut2d
Parameters
----------
lut : xarray.Dataset
Dataset containing an extreme wind lookup table
gewc : xarray.Dataset
Generalized Extreme Wind Climate dataset
Returns
-------
out_ds : xarray.Dataset
Downscaled extreme winds using the lookup table
"""
lut2d_fort = Rvea0293.calc_points_lut2d
lut_pts = spatial_stack(lut)
lut_pts = lut_pts.stack(conditions=["gen_wspd", "sector"]).reset_index("conditions")
lut_pts = lut_pts.transpose("point", "conditions")
gewc_pts = spatial_stack(gewc)
gewc_pts = gewc_pts.transpose("sector", "year", "point")
result = lut2d_fort(
lut_pts.gen_wspd,
lut_pts.sector,
lut_pts.WS,
lut_pts.WD,
30.0,
gewc_pts.max_wspd,
gewc_pts.max_wdir,
)
##########################
# Build the output dataset
##########################
result_names = ("max_wspd", "max_wdir", "ierror")
# Combine the results into a dataset
out_ds = create_dataset(
gewc_pts["west_east"],
gewc_pts["south_north"],
gewc_pts["height"],
get_crs(gewc_pts),
"point",
).drop_vars("output")
out_ds = out_ds.assign_coords(year=gewc_pts.coords["year"])
for i, name in enumerate(result_names):
out_ds[name] = xr.DataArray(result[i], dims=["sector", "year", "point"])
out_ds[name].attrs["_pwio_data_is_2d"] = False
# Add attributes from gewc
out_ds.attrs = {**out_ds.attrs, **lut_pts.attrs}
out_ds.attrs["title"] = "Predicted annual extreme wind"
# Convert to original projection and raster if necessary
out_ds = spatial_unstack(out_ds)
res = update_var_attrs(out_ds, _LINCOM_V50_ATTRS)
return update_history(res)
[docs]
def apply_lut(lut, gewc, gewc_interp="nearest"):
"""Apply lookup table to Generalized Extreme Wind Atlas
.. warning::
This function is experimental and its signature may change.
Parameters
----------
lut : xarray.Dataset
Dataset containing an extreme wind lookup table
gewc : xarray.Dataset
Generalized Extreme Wind Climate dataset
gewc_interp : str, optional
Method to interpolate gewc, by default "nearest"
Returns
-------
xarray.Dataset
Local downscaled extreme winds each location in the lookup table
"""
####################################################
# Interpolate gewc to lut, which is our output grid
####################################################
lut_3d = spatial_stack(lut)
gewc = interpolate_to_grid(gewc, lut_3d, gewc_interp)
pewc = _lut2d(lut_3d, gewc)
pewc.attrs = lut_3d.attrs
pewc = spatial_unstack(pewc)
res = pewc.max("sector", keep_attrs=True)
return update_history(res)