"""LINCOM Wind Level
A dataset that contains the 3 wind components, 9 derivatives, turnings, and
speedups of the wind from LINCOM at a specified height above ground level.
"""
import tempfile
import zipfile
from pathlib import Path
import numpy as np
import pyproj
import xarray as xr
from windkit.spatial import create_dataset
from windkit.metadata import _LINCOM_WIND_LEVEL_ATTRS, update_var_attrs, update_history
from ..core import Rvea0305
FORT_CALC_WL = Rvea0305.windlevel_mod.l1lev # pylint:disable=c-extension-no-member
[docs]
class WindLevel:
"""WindLevel object
The WindLevel contains the wind speed, flow inclination angle and derivatives that
can be used to create result data. In addition, this object can be converted to an
xarray Dataset to visualize some of the key layers as done in WAsP Engineering.
Results are prepared for the entire grid in the raster. If you wish to extract
results for certain points you should use get_wind_points after building your
WindLevel.
For WindLevel wind components, u is along the mean flow direction, and v is
cross flow.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
fou_space : FourierSpace
FourierSpace object FourierSpace object used to create wind level
du_terrain, du_roughness, du_total : numpy.ndarray
1D Array of wind derivatives due to orography, roughness, and the total
Contains 3 derivatives for each component of the 3D wind vector at each point
flow_inc : numpy.ndarray
1D Array of the flow inclination angle
u : numpy.ndarray
1D Array of all 3 components of the wind vector at each point
height : float
height above ground level of results
"""
[docs]
def __init__(
self, fou_space, du_terrain, du_roughness, du_total, flow_inc, u, height
):
# Grid desciption
self.nx_proj = fou_space.nx_proj
self.ny_proj = fou_space.ny_proj
self.nx_comp = fou_space.nx_comp
self.ny_comp = fou_space.ny_comp
self.crs = fou_space.crs
# Wind level results
self.du_terrain = du_terrain # dun
self.du_roughness = du_roughness # dur
self.du_total = du_total # dut
self.flow_inc = flow_inc # ta
self.u = u
self.height = height
# Copies of data from FourierSpace
self.orography = fou_space.orography
self.ustarp = fou_space.ustarp
self.rgh0 = fou_space.rgh0
self.terr_incl = fou_space.hx
self.iarr = fou_space.iarr
self.rarr = fou_space.rarr
self.wind = fou_space.wind
[docs]
@classmethod
def get_wind_level(cls, fou_space, height):
"""Calculate wind level data from fourier space
Parameters
----------
fouspace : FourierSpace
FourierSpace object used to calculate wind level data
height: float
Height to calculate wind level data
Returns
-------
WindLevel
Instance of WindLevel() class
"""
(
du_terrain,
du_roughness,
du_total,
flow_inc,
u,
) = FORT_CALC_WL( # pylint:disable=unused-variable
fou_space.ny_comp,
fou_space.bigl,
fou_space.expnsl,
fou_space.exprsl,
fou_space.un,
fou_space.aun,
fou_space.ur,
height,
fou_space.dhdx,
fou_space.dhdy,
fou_space.iarr,
fou_space.rarr,
)
return cls(fou_space, du_terrain, du_roughness, du_total, flow_inc, u, height)
[docs]
def to_dataset(self):
"""Create a Wind Level dataset to hold the wind information
Returns
-------
xarray.Dataset
Wind Level dataset to hold the wind information
"""
# Format 1D LINCOM arrays to multi-dimensional Numpy arrays
# NOTE: All WindLevel arrays are on the Computational domain
nx = self.nx_comp
ny = self.ny_comp
u_arr = np.reshape(self.u, (3, ny, nx))
flow_ang_arr = np.reshape(self.flow_inc, (ny, nx))
du_tot_arr = np.reshape(self.du_total, (9, ny, nx))
return _to_wl_dataset(u_arr, flow_ang_arr, du_tot_arr, self)
[docs]
def to_file(self, file_name):
"""Save to zip file that can be converted back to wind_level object
Parameters
----------
file_name : str or Path
Path for storing .zip file
"""
with zipfile.ZipFile(file_name, "w") as out_file:
with tempfile.NamedTemporaryFile() as npz_file:
# Wind and crs are special so don't handle them like normal
npz_objs = {
i[0]: i[1]
for i in self.__dict__.items()
if i[0] not in ["wind", "crs"]
}
npz_objs["crs"] = str(self.crs)
# Save all objects except wind to npz file and arcive in zip
np.savez(npz_file, **npz_objs)
out_file.write(npz_file.name, "wind_level.npy")
with tempfile.NamedTemporaryFile() as nc_file:
self.wind.to_netcdf(nc_file.name)
out_file.write(nc_file.name, "wind.nc")
[docs]
@classmethod
def read_file(cls, file_name):
"""Instantiate WindLevel class from .zip file
Parameters
----------
file_name : str or Path
zip file path
Returns
-------
WindLevel
Instance of WindLevel() class
Raises
------
ValueError
If zip file is not a Fourier space archive
"""
with zipfile.ZipFile(file_name, "r") as in_file:
if in_file.namelist() != ["wind_level.npy", "wind.nc"]:
print(in_file.namelist())
raise ValueError("file is not a fourier space archive.")
with tempfile.TemporaryDirectory() as tmpdirname:
tmpdir = Path(tmpdirname)
in_file.extractall(tmpdir)
# Build up dictionary based on numpy file
print(tmpdirname)
npz_data = np.load(tmpdir / "wind_level.npy", "r")
# Integer fields (FourierSpace Stub)
int_fields = [
"nx_proj",
"ny_proj",
"nx_comp",
"ny_comp",
]
fs_dict = {i: int(npz_data[i]) for i in int_fields}
# numpy fields (FourierSpace Stub)
fs_fields = [
"orography",
"ustarp",
"rgh0",
"terr_incl",
"iarr",
"rarr",
]
fs_dict.update({i: npz_data[i] for i in fs_fields})
# Add custom variables
fs_dict.update(
{
"srs": str(npz_data["crs"]),
}
)
# Add wind
wind_ds = xr.load_dataset(tmpdir / "wind.nc")
fs_dict.update({"wind": wind_ds})
fs_stub = _fs_data(fs_dict)
# numpy fields (WindLevel)
wl_fields = [
"du_terrain",
"du_roughness",
"du_total",
"flow_inc",
"u",
]
wl_dict = {i: npz_data[i] for i in wl_fields}
wl_dict.update(
{
"height": float(npz_data["height"]),
}
)
return cls(fs_stub, **wl_dict)
class _fs_data:
"""Dummy class for storing part of fourier space needed to restore from file"""
def __init__(self, from_file_dict):
# Grid desciption
self.nx_proj = from_file_dict["nx_proj"]
self.ny_proj = from_file_dict["ny_proj"]
self.nx_comp = from_file_dict["nx_comp"]
self.ny_comp = from_file_dict["ny_comp"]
self.crs = pyproj.CRS.from_user_input(from_file_dict["srs"])
# Copies of data from FourierSpace
self.orography = from_file_dict["orography"]
self.ustarp = from_file_dict["ustarp"]
self.rgh0 = from_file_dict["rgh0"]
self.hx = from_file_dict["terr_incl"]
self.iarr = from_file_dict["iarr"]
self.rarr = from_file_dict["rarr"]
self.wind = from_file_dict["wind"]
def _to_wl_dataset(u, flow_ang, du_tot, wind_level):
"""Create a Wind Level dataset to hold the wind information
Parameters
----------
u: array (3, ny, nx)
3 wind speed components
flow_ang: array (ny, nx)
Flow inclination angle
du_tot: array (9, ny, nx)
Total derivatives (9 derivatives wind components and dimension)
wind_level: WindLevel
WindLevel to get some of dimensional data
Returns
-------
xarray.Dataset
Dataset containing all the information for the wind level
"""
# TODO: Add Dynamic Roughness, Fetch, Friction Velocity, and Terrain Inclination
nx = wind_level.nx_proj
ny = wind_level.ny_proj
# Create dummy dataset
we = np.linspace(wind_level.rarr[0], wind_level.rarr[2], nx)
sn = np.linspace(wind_level.rarr[1], wind_level.rarr[3], ny)
out_ds = (
create_dataset(we, sn, wind_level.height, wind_level.crs)
.drop_vars("output")
.drop_dims("height")
)
rast_dims = ("south_north", "west_east")
# Add wind speed variables
wspd = np.sqrt(u[0, ...] ** 2 + u[1, ...] ** 2 + u[2, ...] ** 2)
out_ds["WS"] = (rast_dims, wspd[:ny, :nx])
out_ds["U"] = (rast_dims, u[0, :ny, :nx])
out_ds["V"] = (rast_dims, u[1, :ny, :nx])
out_ds["W"] = (rast_dims, u[2, :ny, :nx])
# Flow inclination angle
out_ds["flow_inclination"] = (rast_dims, flow_ang[:ny, :nx])
# Topographic data
out_ds["Z0"] = (rast_dims, np.reshape(wind_level.rgh0, (ny, nx)))
out_ds["terrain_inclination"] = (
rast_dims,
np.reshape(wind_level.terr_incl, (ny, nx)),
)
# Flow derivatives
out_ds["DU_DX"] = (rast_dims, du_tot[0, :ny, :nx])
out_ds["DV_DX"] = (rast_dims, du_tot[1, :ny, :nx])
out_ds["DW_DX"] = (rast_dims, du_tot[2, :ny, :nx])
out_ds["DU_DY"] = (rast_dims, du_tot[3, :ny, :nx])
out_ds["DV_DY"] = (rast_dims, du_tot[4, :ny, :nx])
out_ds["DW_DY"] = (rast_dims, du_tot[5, :ny, :nx])
out_ds["DU_DZ"] = (rast_dims, du_tot[6, :ny, :nx])
out_ds["DV_DZ"] = (rast_dims, du_tot[7, :ny, :nx])
out_ds["DW_DZ"] = (rast_dims, du_tot[8, :ny, :nx])
# Update metadata of created dataset
out_ds = update_var_attrs(out_ds, _LINCOM_WIND_LEVEL_ATTRS)
out_ds.attrs["title"] = "LINCOM wind level"
return update_history(out_ds)