import numpy as np
from scipy.interpolate import SmoothBivariateSpline
from scipy.spatial import Delaunay
from logging import warning, debug
from .utils import \
_latlon_to_xyz, _spherical_delaunay_global, \
_spherical_delaunay_regional, _tangent_point
[docs]
def bilinear(domain_in,
nlon_in, nlat_in,
zlon_in, zlat_in,
zlon_out, zlat_out,
use_delaunay=True):
# Extend input domain with outside corners for structured domains
if not getattr(domain_in, "unstructured_domain", False):
zlon_ref = np.zeros((nlat_in + 2, nlon_in + 2))
zlon_ref[1:-1, 1:-1] = zlon_in
zlon_ref[1:-1, 0] = domain_in.zlon_side[0, :nlat_in]
zlon_ref[1:-1, -1] = domain_in.zlon_side[0, nlat_in:2 * nlat_in]
zlon_ref[0, 1:-1] = domain_in.zlon_side[
0, 2 * nlat_in:2 * nlat_in + nlon_in]
zlon_ref[-1, 1:-1] = domain_in.zlon_side[0, 2 * nlat_in + nlon_in:]
zlon_ref[0, 0] = domain_in.zlonc_side[0, 0]
zlon_ref[-1, 0] = domain_in.zlonc_side[nlat_in, 0]
zlon_ref[0, -1] = domain_in.zlonc_side[nlat_in + 1, 0]
zlon_ref[-1, -1] = domain_in.zlonc_side[2 * nlat_in + 1, 0]
zlat_ref = np.zeros((nlat_in + 2, nlon_in + 2))
zlat_ref[1:-1, 1:-1] = zlat_in
zlat_ref[1:-1, 0] = domain_in.zlat_side[0, :nlat_in]
zlat_ref[1:-1, -1] = domain_in.zlat_side[0, nlat_in:2 * nlat_in]
zlat_ref[0, 1:-1] = domain_in.zlat_side[0,
2 * nlat_in:2 * nlat_in + nlon_in]
zlat_ref[-1, 1:-1] = domain_in.zlat_side[0, 2 * nlat_in + nlon_in:]
zlat_ref[0, 0] = domain_in.zlatc_side[0, 0]
zlat_ref[-1, 0] = domain_in.zlatc_side[nlat_in, 0]
zlat_ref[0, -1] = domain_in.zlatc_side[nlat_in + 1, 0]
zlat_ref[-1, -1] = domain_in.zlatc_side[2 * nlat_in + 1, 0]
else:
zlon_ref = zlon_in
zlat_ref = zlat_in
# Keeping coordinates between -180 and +180
mask_out = zlon_out == 180
zlon_out = (np.asarray(zlon_out) + 180) % 360 - 180
zlon_out[mask_out] = 180
mask_in = zlon_ref == 180
zlon_ref = (zlon_ref + 180) % 360 - 180
zlon_ref[mask_in] = 180
# Unwraps zlon_ref
if not getattr(domain_in, "unstructured_domain", False):
zlon_ref = np.degrees(np.unwrap(np.radians(zlon_ref), axis=1))
zlon_ref -= 360 * np.round(zlon_ref.mean() / 360)
# For regular domains
if not getattr(domain_in, "unstructured_domain", False):
lon_cyclic = getattr(domain_in, "lon_cyclic", False)
# Fractional indexes along each axis; the domain sides are given
# the indexes -1 and n, which are clipped to the outermost centers
lat_axis = zlat_ref[:, 0]
ind_lat = np.arange(-1, nlat_in + 1)
if np.all(np.diff(lat_axis) < 0):
lat_axis = lat_axis[::-1]
ind_lat = ind_lat[::-1]
# Cyclic domains are closed with the centers of the first and last
# columns on the other side of the periodic boundary
if lon_cyclic:
lon_centers = zlon_ref[0, 1:-1]
lon_axis = np.concatenate(
[[lon_centers[-1] - 360], lon_centers, [lon_centers[0] + 360]])
else:
lon_axis = zlon_ref[0]
ind_lon = np.arange(-1, nlon_in + 1)
# Bring target longitudes in the same revolution as the input domain
lon_out = np.asarray(zlon_out).flatten(order="F")
lat_out = np.asarray(zlat_out).flatten(order="F")
lon_out = lon_axis[0] + (lon_out - lon_axis[0]) % 360
meshi_out = np.interp(lat_out, lat_axis, ind_lat)
meshj_out = np.interp(lon_out, lon_axis, ind_lon)
# Lower corners of the interpolation cells, upper corners being +1
meshi_out = np.clip(meshi_out, 0, nlat_in - 1)
imin = np.minimum(np.floor(meshi_out), nlat_in - 2).astype(int)
if lon_cyclic:
jmin = np.minimum(np.floor(meshj_out), nlon_in - 1).astype(int)
else:
meshj_out = np.clip(meshj_out, 0, nlon_in - 1)
jmin = np.minimum(np.floor(meshj_out), nlon_in - 2).astype(int)
alpha_imax = meshi_out - imin
alpha_jmax = meshj_out - jmin
# Put the data into a dictionary
weights = {
"i": np.concatenate(
(imin[:, np.newaxis],
imin[:, np.newaxis],
imin[:, np.newaxis] + 1,
imin[:, np.newaxis] + 1),
axis=1),
"j": np.concatenate(
(jmin[:, np.newaxis],
jmin[:, np.newaxis] + 1,
jmin[:, np.newaxis] + 1,
jmin[:, np.newaxis]),
axis=1),
"wgt": np.concatenate(
((1 - alpha_imax[:, np.newaxis]) * (1 - alpha_jmax[:, np.newaxis]),
(1 - alpha_imax[:, np.newaxis]) * alpha_jmax[:, np.newaxis],
alpha_imax[:, np.newaxis] * alpha_jmax[:, np.newaxis],
alpha_imax[:, np.newaxis] * (1 - alpha_jmax[:, np.newaxis])),
axis=1)
}
# Puts NaNs for data outside the domain for regular domains
inside = (lat_out >= lat_axis[0]) & (lat_out <= lat_axis[-1])
if not lon_cyclic:
inside = inside & (lon_out <= lon_axis[-1])
weights = {
"i": weights["i"][inside],
"j": weights["j"][inside],
"wgt": weights["wgt"][inside],
"filtered": np.where(inside)[0],
"non_filtered": np.where(~inside)[0]
}
# For cyclic domains in longitude, use modulo of j-index
if lon_cyclic:
weights["j"] = weights["j"] % nlon_in
# For unstructured domains
else:
if domain_in.projection == "xy":
debug("Do Delaunay triangularisation for interpolating")
points_ref = np.concatenate([zlon_ref, zlat_ref], axis=0).T
n_ref = points_ref.shape[0]
if getattr(domain_in, "lon_cyclic", False):
# Duplicate the reference points one full turn either side
# so the triangulation wraps seamlessly across the periodic
# boundary instead of treating it as the edge of the hull;
# indices into the duplicated copies are folded back to the
# original reference index (mod n_ref) below.
points_ref = np.concatenate([
points_ref + [-360, 0],
points_ref,
points_ref + [360, 0],
], axis=0)
triangular = Delaunay(points_ref)
lon_out = np.asarray(zlon_out).flatten(order="F")
lat_out = np.asarray(zlat_out).flatten(order="F")
points_out = np.stack([lon_out, lat_out], axis=1)
simplexes = triangular.find_simplex(points_out)
# Points falling outside the reference convex hull (e.g. poles
# not covered by the reference grid) get find_simplex == -1;
# filter them out rather than silently wrapping to whatever
# triangle happens to sit last in the array.
inside = simplexes >= 0
simplexes_in = simplexes[inside]
barycenters = np.matmul(
triangular.transform[simplexes_in, :2],
(points_out[inside]
- triangular.transform[simplexes_in, 2])[..., np.newaxis]
)[..., 0]
coeffs = np.concatenate(
[barycenters, 1 - barycenters.sum(axis=1)[:, np.newaxis]],
axis=1)
i_out = triangular.simplices[simplexes_in] % n_ref
# Put the data into a dictionary
weights = {
"i": 0 * i_out,
"j": i_out,
"wgt": coeffs,
"filtered": np.where(inside)[0],
"non_filtered": np.where(~inside)[0]
}
else:
debug("Do spherical Delaunay triangulation for interpolating")
lon_ref = np.asarray(zlon_ref).flatten()
lat_ref = np.asarray(zlat_ref).flatten()
points_ref_3d = _latlon_to_xyz(lon_ref, lat_ref)
lon_out = np.asarray(zlon_out).flatten(order="F")
lat_out = np.asarray(zlat_out).flatten(order="F")
points_out_3d = _latlon_to_xyz(lon_out, lat_out)
# A single gnomonic chart cannot cover more than a hemisphere,
# so the reference grid decides the algorithm: regional grids
# get the chart plus ordinary planar Delaunay (better
# triangulation quality, exact out-of-domain rejection at the
# data boundary), global grids get the 3D convex hull (no
# chart, hence no antimeridian seam and no pole degeneracy).
center, _ = _tangent_point(points_ref_3d)
if center is not None:
debug("Reference grid is regional: gnomonic chart")
i_out, coeffs, inside = _spherical_delaunay_regional(
points_ref_3d, points_out_3d, center)
else:
debug("Reference grid is global: spherical convex hull")
i_out, coeffs, inside = _spherical_delaunay_global(
points_ref_3d, points_out_3d)
if not inside.all():
debug(
"{} of {} output points fall outside the reference "
"grid and are excluded".format(
int((~inside).sum()), inside.size)
)
# Put the data into a dictionary
weights = {
"i": 0 * i_out,
"j": i_out,
"wgt": coeffs,
"filtered": np.where(inside)[0],
"non_filtered": np.where(~inside)[0]
}
# print(__file__)
# import code
# code.interact(local=dict(locals(), **globals()))
return weights
