"""Python interface to LINCOM's FourierSpace
The FourierSpace created by LINCOM is quite expensive to calculation, but only depends
on the topographic data, and the input wind, so it can be stored and re-used to calculate
different WindLevels using vertical extrapolation.
"""
import logging
import tempfile
import zipfile
from pathlib import Path
import numpy as np
import pyproj
import xarray as xr
from windkit import spatial as wksp
# import fortran libraries
from windkit.raster_map import _landcover_to_roughness
from ..core import Rvea0305
from .._cfg import _config as usr_config
FORT_SET_DOMAIN = Rvea0305.dim_mod.setdim1
FORT_RUN_LINCOM = Rvea0305.fourier_mod.l1fou
[docs]
class FourierSpace:
"""FourierSpace object
The FourierSpace contains the result from l1fou as well as some of the
input data that was requested. This is used to create wind level at
specified heights above ground level.
Parameters
----------
nx_proj, ny_proj : int
Number of points in the x and y dimension of the project space
nx_comp, ny_comp : int
Number of points in the x and y dimensions of the computational space. Currently
these are 1.5 times the project dimensions.
bigl : numpy.ndarray
OuterLengthScale (1/2 comp domain size)
expnsl : numpy.ndarray
Exponential height factor (1/2 comp domain size)
exprsl : numpy.ndarray
Exponential roughness factor (1/2 comp domain size)
un : numpy.ndarray
Terrain elevation perturbation (3 times comp domain size)
aun : numpy.ndarray
Terrain elevation perturbation (1/2 comp domain size)
ur : numpy.ndarray
Roughness perturbation (3 times comp domain size)
ustarp : numpy.ndarray
Friction Velocity Grid (comp domain size)
orography : numpy.ndarray
Terrain grid (Project domain size)
hx : numpy.ndarray
Terrain inclination wind speed direction in degrees (project domain size)
dhdx : numpy.ndarray
Terrain inclination in x-direction (comp domain size)
dhdy : numpy.ndarray
Terrain inclination in y-direction (comp domain size)
rgh0 : numpy.ndarray
Roughness updated with water values (project domain size)
iarr : numpy.ndarray
Array of integer settings
rarr : numpy.ndarray
Array of floating-point settings
wind : xr.Dataset
Input wind used to create fourier space
srs : CRS object
PyProj.CRS of projection used in simulation
Returns
-------
FourierSpace
Object with the requested settings
"""
[docs]
def __init__(
self,
nx_proj,
ny_proj,
nx_comp,
ny_comp,
bigl,
expnsl,
exprsl,
un,
aun,
ur,
ustarp,
orography,
hx,
dhdx,
dhdy,
rgh0,
iarr,
rarr,
wind,
srs,
):
self.nx_proj = nx_proj # nx_proj
self.ny_proj = ny_proj # ny_proj
self.nx_comp = nx_comp # nx_comp
self.ny_comp = ny_comp # ny_comp
self.bigl = bigl # OuterLengthScale (1/2 comp domain size)
self.expnsl = expnsl # Exponential height factor (1/2 comp domain size)
self.exprsl = exprsl # Exponential roughness factor (1/2 comp domain size)
self.un = un # Terrain elevation perturbation (3 times comp domain size)
self.aun = aun # Terrain elevation perturbation (1/2 comp domain size)
self.ur = ur # Roughness perturbation (3 times comp domain size)
self.ustarp = ustarp # Friction Velocity Grid (comp domain size)
self.orography = orography # Terrain grid (Project domain size)
self.hx = hx # Terrain inclination wind speed direction in degrees (project domain size)
self.dhdx = dhdx # Terrain inclination in x-direction (comp domain size)
self.dhdy = dhdy # Terrain inclination in y-direction (comp domain size)
self.rgh0 = rgh0 # Roughness updated with water values (project domain size)
self.iarr = iarr # Array of integer settings (23)
self.rarr = rarr # Array of floating-point settings (36)
self.wind = wind # Input wind used to create fourier space
self.crs = pyproj.CRS.from_user_input(
srs
) # PyProj.CRS of projection used in simulation
[docs]
@classmethod
def create_fourier_space(
cls, ter_grid, lc_grid, lctable, fetch_grid, wind, const_z0=0.0
):
"""Creates a Fourier space object
.. warning::
This function is experimental and its signature may change.
Parameters
----------
ter_grid : xarray.DataArray
Terrain windkit.raster_map
lc_grid : xarray.DataArray
landcover data as a rastermap
lctable : LandCoverTable
Mapping of landcovers to roughness
fetch_grid : grid.GridMap
Fetch grid windkit.raster_map
wind : wind.Wind
Wind description
Notes
-----
All GridMaps need to have the same dimensions and spatial extent
"""
### Check size of terr_grid
usr_config.check_point_count(wksp.count_spatial_points(ter_grid))
# Toggle for different calculation approaches
ifterra = 1 # Perform terrain calculation?
ifrough = 1 # Use provided z0 map? false uses a constant z0 map.
ifwtrgh = 1 # Generate water roughness table?
ifbackg = 1 # Adds an additional Ustar0
ifderiv = 1 # Doesn't seem to be used
ilinfld = 2 # Also not seemed to be used
# Setup the domain dimensions. Project domain is from the terrain grid, but the
# computational domain has a hard-coded 50% buffer in WEng.
nx_proj = ter_grid.sizes["west_east"]
ny_proj = ter_grid.sizes["south_north"]
nx_comp, ny_comp = FORT_SET_DOMAIN(nx_proj, ny_proj)
x = ter_grid["west_east"].values
y = ter_grid["south_north"].values
dx = x[1] - x[0]
dy = y[1] - y[0]
rou_grid = _landcover_to_roughness(lc_grid, lctable, "z0")
# Call LINCOM to create the Fourier space
(
bigl,
expnsl,
exprsl,
un,
aun,
ur,
ustarp,
hx,
dhdx,
dhdy,
rgh0,
iarr,
rarr,
ierror,
) = FORT_RUN_LINCOM(
nx_comp,
ny_comp,
nx_comp * ny_comp,
nx_proj,
ny_proj,
ifterra,
ifrough,
ifwtrgh,
ifbackg,
ifderiv,
ilinfld,
ter_grid.west_east.min(),
ter_grid.south_north.min(),
dx,
dy,
wind.latitude,
const_z0,
wind.type,
wind.wind_speed,
wind.sector,
wind.wind_height,
wind.z0,
wind.west_east,
wind.south_north,
ter_grid.values.flatten(),
rou_grid.flatten(), # rghmap
fetch_grid.values.flatten(),
)
if ierror != 0:
raise ValueError("License Error: ", ierror)
# Create a dictionary representation of the fourier space
return cls(
nx_proj,
ny_proj,
nx_comp,
ny_comp,
bigl,
expnsl,
exprsl,
un,
aun,
ur,
ustarp,
ter_grid.values.flatten(),
hx,
dhdx,
dhdy,
rgh0,
iarr,
rarr,
wind,
wksp.get_crs(ter_grid),
)
[docs]
def to_file(self, file_name):
"""Save to zip file that can be converted back to fs object
.. warning::
This function is experimental and its signature may change.
Parameters
----------
file_name : str or Path
Path for storing Furier space object az a .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, "fs.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 from_file(cls, file_name):
"""Instantiate FourierSpace class from a .zip file
.. warning::
This function is experimental and its signature may change.
Parameters
----------
file_name : str or Path
Path of a .zip file which contains Fourier space object
Returns
-------
FourierSpace
Instantiated FourierSpace class
Raises
------
ValueError
If a .zip file does not containing FourierSpace object
"""
logging.info("Reading FourierSpace from %s", file_name)
with zipfile.ZipFile(file_name, "r") as in_file:
if in_file.namelist() != ["fs.npy", "wind.nc"]:
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
logging.debug("ZipFile Extracted to %s.", tmpdirname)
npz_data = np.load(tmpdir / "fs.npy", "r")
# Integer fields
int_fields = [
"nx_proj",
"ny_proj",
"nx_comp",
"ny_comp",
]
fs_dict = {i: int(npz_data[i]) for i in int_fields}
# numpy fields
np_fields = [
"bigl",
"expnsl",
"exprsl",
"un",
"aun",
"ur",
"ustarp",
"orography",
"hx",
"dhdx",
"dhdy",
"rgh0",
"iarr",
"rarr",
]
fs_dict.update({i: npz_data[i] for i in np_fields})
# Add crs
fs_dict.update({"srs": str(npz_data["crs"])})
# Add wind
wind_ds = xr.load_dataset(tmpdir / "wind.nc")
fs_dict.update({"wind": wind_ds})
return cls(**fs_dict)