"""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 pywasp._windkit_internal import (
_GEWC_ATTRS,
_spatial_dims_in_order,
_update_history,
_update_var_attrs,
)
from pywasp._batch import batch_layout
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
"""
# find_natural_neighbors requires a point/stacked_point source; gewc's
# own structure is never restored, so to_point's flattening is enough.
out_crs = wksp.get_crs(lut)
gewc_stk = wksp.reproject(wksp.to_point(gewc), out_crs)
# Create a dummy output dataset with only x,y coordinates. Never height:
# natural-neighbor interpolation is purely horizontal, and west_east/
# south_north never depend on height regardless of lut's own structure.
out_locs = lut[["south_north", "west_east"]]
# Locate the natural neighbors
members, delauny_tri, tri_info = _nat_neighbor_interp.find_natural_neighbors(
gewc_stk, out_locs
)
# Filter to only members that have neighbors
valid_members = [(key, val) for key, val in members.items() if val]
# The anchor is opened on out_locs's own (horizontal-only) structure, so
# restore() below hands the result back in lut's native structure
# directly -- no unstack needed.
ctx = batch_layout(
out_locs,
batch_dims=_spatial_dims_in_order(out_locs),
core_dims=(),
required_coords=("west_east", "south_north"),
)
# `members` is keyed by find_natural_neighbors's own
# wksp.to_stacked_point(out_locs), which deduplicates repeated horizontal
# locations in order of first occurrence (e.g. a point structure with
# several heights per location). Reproduce that same order-preserving
# dedup here so ctx's full, possibly duplicated points line up with those
# keys; a no-op when out_locs has no duplicate locations to begin with.
we, sn = ctx.flat_coords(dtype=np.float64)
dedup_xy = list(dict.fromkeys(zip(we, sn)))
dedup_index = {xy: i for i, xy in enumerate(dedup_xy)}
to_dedup = np.array([dedup_index[xy] for xy in zip(we, sn)])
grid_points = np.array(dedup_xy)
# Add output arrays
out_shp = (len(dedup_xy), gewc_stk.sizes["sector"], gewc_stk.sizes["year"])
max_wspd = np.full(out_shp, -999.0, np.float32)
max_wdir = np.full(out_shp, -999.0, np.float32)
# 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
# Interpolate each member that has neighbors
for k_member, v_member in valid_members:
wgts, nn_pts = _nat_neighbor_interp.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]
near_pt = nn_pts[wgts.argmax()]
wgts = wgts.reshape(-1, 1, 1)
# Calculate weighted average
mean_wspd = (gwc_neigh_wspd * wgts).sum(0) / wgts.sum()
max_wspd[k_member, ...] = mean_wspd
max_wdir[k_member, ...] = gwc_max_wdir[near_pt]
result = ctx.to_dataset(
{
"max_wspd": ctx.restore_variable(
np.moveaxis(max_wspd[to_dedup], 0, -1), core_dims=_CORE_DIMS
),
"max_wdir": ctx.restore_variable(
np.moveaxis(max_wdir[to_dedup], 0, -1), core_dims=_CORE_DIMS
),
},
coords={"sector": gewc_stk["sector"], "year": gewc_stk["year"]},
)
# Update metadata
_update_var_attrs(result, _GEWC_ATTRS)
result.attrs["title"] = "Generalized Extreme Wind Climate"
return result
[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)
