import numpy as np
from scipy.interpolate import RectBivariateSpline, SmoothBivariateSpline
from scipy.spatial import Delaunay
from logging import warning, debug
[docs]
def bilinear(domain_in,
nlon_in, nlat_in,
zlon_in, zlat_in,
zlon_out, zlat_out):
# 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:
debug("Do Delaunay triangularisation for interpolating")
points_ref = np.concatenate([zlon_ref, zlat_ref], axis=0).T
triangular = Delaunay(points_ref)
lon_out = np.asarray(zlon_out).flatten(order="F")[np.newaxis]
lat_out = np.asarray(zlat_out).flatten(order="F")[np.newaxis]
points_out = np.concatenate([lon_out, lat_out], axis=0).T
simplexes = triangular.find_simplex(points_out)
barycenters = np.matmul(
triangular.transform[simplexes, :2],
(points_out
- triangular.transform[simplexes, 2])[..., np.newaxis]
)[..., 0]
coeffs = np.concatenate(
[np.transpose(barycenters),
1 - barycenters.sum(axis=1)[np.newaxis]], axis=0).T
i_out = triangular.simplices[simplexes]
# Put the data into a dictionary
weights = {
"i": 0 * i_out,
"j": i_out,
"wgt": coeffs
}
# TODO: do something for unregular domains when data outside convex hull
return weights
