"""Routines to interpolate generalized extreme wind climate
Only Natural Neighbor interpolation is implemented at this time, and as it is an almost
pure python implemention it is quite slow. It is recommended that it is only used on
relatively small domains.
"""
import numpy as np
from windkit import spatial as wksp
from windkit.metadata import _GEWC_ATTRS, update_var_attrs, update_history
import xarray as xr
from ..io import _nat_neighbor_interp
_CORE_DIMS = ("sector", "year")
def _natn_interp_gewc(gewc, lut):
"""Natural neighbor interpolation based on GWC approach
Parameters
----------
gewc : xarray.Dataset
PyWAsP Generalized Extreme Wind Climate Dataset
lut : xarray.Dataset
PyWAsP Extreme wind lookup table Dataset (Only used for spatial information)
Returns
-------
xarray.Dataset
gewc interpolated to the spatial region of the lut, including possible reproj
"""
natn = _nat_neighbor_interp
# Prepare GWC object for processing
out_crs = wksp.get_crs(lut)
gewc_stk = wksp.spatial_stack(gewc, out_crs, False)
# Create a dummy output dataset with only x,y coordinates
out_locs = lut[["south_north", "west_east"]]
# Locate the natural neighbors
members, delauny_tri, tri_info = natn.find_natural_neighbors(gewc_stk, out_locs)
# Drop members that don't have any neighbors
gen = [(key, val) for (key, val) in members.items() if len(val) != 0]
# Create 2d list of points in output grid [x, y]
out_locs_pt = wksp.spatial_stack(out_locs)
grid_points = np.vstack([out_locs_pt.west_east, out_locs_pt.south_north]).T
# Add output arrays
out_dims = ("point", "sector", "year")
out_shp = (
out_locs_pt.sizes["point"],
gewc_stk.sizes["sector"],
gewc_stk.sizes["year"],
)
out_locs_pt["max_wspd"] = (out_dims, np.full(out_shp, -999.0, np.float32))
out_locs_pt["max_wdir"] = (out_dims, np.full(out_shp, -999.0, np.float32))
# out_locs_pt["max_wdir_mean"] = (out_dims, np.full(out_shp, -999., np.float32))
out_locs_pt = out_locs_pt.assign_coords(
{
"sector": gewc_stk.sector,
"year": gewc_stk.year,
}
)
# Preprocess values for speed
gewc_stk = gewc_stk.transpose("point", "sector", "year")
gwc_max_wspd = gewc_stk.max_wspd.values
gwc_max_wdir = gewc_stk.max_wdir.values
# loop over each member that has neighbors
for k_member, v_member in gen:
# Calculate the weighting from the natural neighbors
# (This is super expensive, anything to speed it up would help)
wgts, nn_pts = natn.nn_weights(
gewc_stk.west_east,
gewc_stk.south_north,
grid_points[k_member],
delauny_tri,
v_member,
tri_info,
)
gwc_neigh_wspd = gwc_max_wspd[nn_pts]
# gwc_neigh_wdir = gwc_wdir_cmplx[nn_pts]
near_pt = nn_pts[wgts.argmax()]
wgts = np.reshape(wgts, (wgts.size, 1, 1))
# Calcualte weighted average
mean_wspd = (gwc_neigh_wspd * wgts).sum(0) / wgts.sum()
# mean_wdir = np.mod(
# np.angle((gwc_neigh_wdir * wgts).sum(0) / wgts.sum(), True), 360
# )
out_locs_pt["max_wspd"][k_member, ...] = mean_wspd
out_locs_pt["max_wdir"][k_member, ...] = gwc_max_wdir[near_pt]
# out_locs_pt["max_wdir_mean"][k_member, ...] = mean_wdir
out_locs_pt = wksp.spatial_unstack(out_locs_pt)
# Update metadata
update_var_attrs(out_locs_pt, _GEWC_ATTRS)
out_locs_pt.attrs["title"] = "Generalized Extreme Wind Climate"
return out_locs_pt
[docs]
def interpolate_gewc(gewc, target, resampling="natural"):
"""Interpolate GEWC file
.. warning::
This function is experimental and its signature may change.
Parameters
----------
gewc : xarray.Dataset
PyWAsP Generalized Extreme Wind Climate Dataset
target : xarray.Dataset
PyWAsP dataset containing target south_north, west_east, and height coordinates
resampling : str
Only natural neighbors ("natural") is implemented at this time. You can set this
to None if you do not need interpolation because the grids match.
Returns
-------
xarray.Dataset
gewc resampled to the target spatial structure
"""
if resampling.lower() == "natural":
res = _natn_interp_gewc(gewc, target)
return update_history(res)
elif resampling is not None:
raise NotImplementedError(
"Only natural neighbors interpolation is implemented."
)
if not wksp.are_spatially_equal(
wksp.to_stacked_point(gewc), wksp.to_stacked_point(target)
):
raise ValueError(
"gewc and target grids are different, "
+ "so you need to specify a resampling method."
)
return update_history(gewc)