Source code for pywasp.lincom.wind_level

"""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)