import numpy as np
from scipy.interpolate import RectBivariateSpline, 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]
meshj, meshi = np.meshgrid(range(nlon_in + 2), range(nlat_in + 2))
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)
# Shifting first input columns to be West of target domain
zlon_ref = np.unwrap(zlon_ref, axis=1, period=360)
if np.min(zlon_ref) > np.min(zlon_out):
zlon_ref += 360 * \
np.floor((np.min(zlon_out) - np.min(zlon_ref)) / 360)
# Expand cyclic input domains to overlap with output
if not getattr(domain_in, "lon_cyclic", False) and not getattr(domain_in, "unstructured_domain", False):
# Extending zlon_ref zonally
ind_start = np.where(
np.all((zlon_ref + 360) > zlon_ref.max(), axis=0))[0][0]
zlon_ref = np.concatenate(
[zlon_ref, zlon_ref[:, ind_start:] + 360], axis=1)
meshj = np.concatenate([meshj, meshj[:, ind_start:] + nlon_in], axis=1)
meshi = np.concatenate([meshi, meshi[:, ind_start:]], axis=1)
# For regular domains
if not getattr(domain_in, "unstructured_domain", False):
# Reverse latitudes if not increasing
lat_in = zlat_ref[:, 0]
if np.all(np.diff(lat_in) < 0):
lat_in = lat_in[::-1]
meshi = nlat_in + 2 - meshi - 1
interpi = RectBivariateSpline(
zlon_ref[0],
lat_in,
meshi.T,
bbox=[None, None, None, None],
)
interpj = RectBivariateSpline(
zlon_ref[0],
lat_in,
meshj.T,
bbox=[None, None, None, None],
)
# Carry out the interpolation
lon_out = np.asarray(zlon_out).flatten(order="F")
lat_out = np.asarray(zlat_out).flatten(order="F")
meshi_out = interpi(lon_out, lat_out, grid=False)
meshj_out = interpj(lon_out, lat_out, grid=False)
# Shift meshi_out and meshj_out to account for extension of input domain
meshi_out -= 1
meshi_out = np.maximum(meshi_out, 0)
meshi_out = np.minimum(meshi_out, nlat_in - 2)
meshj_out -= 1
meshj_out = np.maximum(meshj_out, 0)
meshj_out = np.minimum(meshj_out, nlon_in - 2)
# Compute weights
imin = np.floor(meshi_out).astype(int)
jmin = np.floor(meshj_out).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 >= zlat_ref.min()) & (lat_out <= zlat_ref.max())
if not getattr(domain_in, "lon_cyclic", False):
inside = inside \
& (lon_out >= zlon_ref.min()) & (lon_out <= zlon_ref.max())
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 getattr(domain_in, "lon_cyclic", False):
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
