import numpy as np
import scipy
[docs]
def vertical_interp(
pres_in, dpres_in,
pres_out, dpres_out,
cropstrato,
vinterp_type="weight", weights_nsubsteps=20
):
"""Compute the interpolation coefficients to apply
a linear vertical interpolation"""
nvertin, lenin = pres_in.shape
nvertout, lenout = pres_out.shape
if lenin != lenout:
raise Exception(
"The input mesh has a different length than the target"
)
# Initializing index mesh for scipy.griddata
meshin = np.array((nvertin + 2) * [list(range(lenin))])
meshout = np.array(nvertout * [list(range(lenout))])
levmeshin = np.array((lenin) * [list(range(-1, nvertin + 1))]).T
levin = levmeshin.flatten()
# Replacing infinite by max/min +/- eps
pin_max = pres_in[(pres_in > -np.inf) & (pres_in < np.inf)].max()
pin_min = pres_in[(pres_in > -np.inf) & (pres_in < np.inf)].min()
pout_max = pres_out[(pres_out > -np.inf) & (pres_out < np.inf)].max()
pout_min = pres_out[(pres_out > -np.inf) & (pres_out < np.inf)].min()
eps = max(np.nanmax(np.abs(np.diff(pres_in, axis=0))),
np.nanmax(np.abs(np.diff(pres_out, axis=0))))
pmin = min(pin_min, pout_min) - eps
pmax = max(pin_max, pout_max) + eps
pres_in[pres_in == -np.inf] = pmin
pres_in[pres_in == np.inf] = pmax
pres_out[pout_min == -np.inf] = pmin
pres_out[pout_min == np.inf] = pmax
# Assuring press_in is in ascending order.
flip_in = False
if not np.all(np.diff(pres_in, axis=0) > 0):
pres_in = np.flip(pres_in, axis=0)
dpres_in = np.flip(dpres_in, axis=0)
levmeshin = np.flip(levmeshin, axis=0)
flip_in = True
# Assuring press_out is in ascending order.
flip_out = False
if not np.all(np.diff(pres_out, axis=0) > 0):
pres_out = np.flip(pres_out, axis=0)
dpres_out = np.flip(dpres_out, axis=0)
flip_out = True
# Initializing pressure mesh for scipy.griddata
# For input pressures (those of the model),
# extending the roof and floor to detect ak points
# outside of the hull of model pressures
pres_in = np.concatenate(
[
pmin * np.ones((1, lenin)),
pres_in,
pmax * np.ones((1, lenout)),
],
axis=0,
)
# Checking assertion made by numpy.interp method.
assert np.all((np.diff(pres_in, axis=0) > 0) | (pres_in[:-1] == 0))
# Convert to layer interface if pressure-weighted interpolation
if vinterp_type == "weight":
# Update input pressures to interfaces
pres_in = np.concatenate([
np.log(np.where(
np.exp(pres_in[1:-1]) - dpres_in / 2 <= 0,
max(1e-10, np.exp(pmin)),
np.exp(pres_in[1:-1]) - dpres_in / 2
)),
np.log(np.where(
np.exp(pres_in[[-2]]) + dpres_in[[-2]] / 2 <= 0,
max(1e-10, np.exp(pmin)),
np.exp(pres_in[[-2]]) + dpres_in[[-2]] / 2)),
pmax * np.ones((1, lenout))
])
pres_in = np.concatenate([
min(pres_in.min(), pmin) * np.ones((1, lenout)),
pres_in
])
levmeshin = np.concatenate([
np.zeros((1, lenin), dtype=int) + nvertin + 1,
levmeshin
])
# Update output pressures to interfaces
pres_out = np.concatenate([
np.log(np.where(
np.exp(pres_out) - dpres_out / 2 <= 0,
max(1e-10, np.exp(pmin)),
np.exp(pres_out) - dpres_out / 2
)),
np.log(np.exp(pres_out[[-1]]) + dpres_out[[-1]] / 2)
])
# Now duplicates output pressures to sub-levels
levmeshout = np.arange(
0, nvertout + 1)[:, np.newaxis] * np.ones((1, lenout))
levmeshout_wgt = np.linspace(
0, nvertout, weights_nsubsteps * nvertout,
endpoint=False)[:, np.newaxis] * np.ones((1, lenout))
pres_out_wgt = np.empty_like(levmeshout_wgt)
for i in range(lenin):
pres_out_wgt[:, i] = np.log(np.interp(
levmeshout_wgt[:, i],
levmeshout[:, i], np.exp(pres_out[:, i])
))
pres_out = pres_out_wgt
# For linear interpolation return alpha high/low
levout = np.empty_like(pres_out)
for i in range(lenin):
levout[:, i] = np.interp(
pres_out[:, i],
pres_in[:, i], levmeshin[:, i]
)
# Getting interpolation coefficients
xlow = np.floor(levout).astype(int)
xhigh = xlow + 1
alphahigh = levout - xlow
alphalow = 1 - alphahigh
# Cropping ak pressure outside of model pressures
toolow = xlow <= 0
toohigh = xhigh > nvertin - 1
xlow[toolow] = 0
xhigh[toolow] = np.maximum(xhigh[toolow], 0)
xlow[toohigh] = np.minimum(xlow[toohigh], nvertin - 1)
xhigh[toohigh] = nvertin - 1
alphalow[toolow] = 1
alphahigh[toolow] = 0
alphalow[toohigh] = 0
alphahigh[toohigh] = 1
if cropstrato:
alphalow[toohigh] = 0
alphahigh[toohigh] = 0
# Scale coefficient by number of sub-step if weight interpolation
if vinterp_type == "weight":
mask = (alphalow > 0) | (alphahigh > 0)
alphahigh[:] = 0
alphalow[mask] = 1 / weights_nsubsteps
# Flip back outputs if pressure out has been flipped
if flip_out:
xlow = xlow[::-1]
xhigh = xhigh[::-1]
alphalow = alphalow[::-1]
alphahigh = alphahigh[::-1]
return xlow, xhigh, alphalow, alphahigh
