由于这是大气数据，我猜想您的网格没有均匀的间距;然而，如果你的网格是直线的（比如每个垂直列都有相同的z坐标），那么你有一些选择。
例如，如果您只需要线性插值（比如简单的可视化），您可以执行以下操作：

# Find nearest grid point
idx = grid[:,0,0].searchsorted(value)
upper = grid[idx,0,0]
lower = grid[idx - 1, 0, 0]
s = (value - lower) / (upper - lower)
result = (1-s) * data[idx - 1, :, :] + s * data[idx, :, :]

(You当然，需要添加value是否超出范围的检查）。对于您的网格大小，这将是非常快的（如在一秒的微小分数）
如果需要的话，你可以很容易地修改上面的函数来执行三次插值;挑战在于为不均匀的垂直间隔选择正确的权重。
使用www.example.com_coordinates的问题在于，尽管它提供了更高阶的插值并可以处理任意采样点，但它确实假设输入数据是均匀分布的。它仍然会产生平滑的结果，但不是可靠的近似值。scipy.ndimage.map_coordinates is that, although it provides higher order interpolation and can handle arbitrary sample points, it does assume that the input data be uniformly spaced. It will still produce smooth results, but it won't be a reliable approximation.
如果您的坐标网格不是直线的，因此给定索引的z值会随不同的x和y索引而变化，那么您现在使用的方法可能是最好的方法，而无需对特定问题进行大量分析。
更新：
一个巧妙的技巧（再次假设每列具有相同的、不一定规则的坐标）是使用interp1d来提取权重，如下所示：

NZ = grid.shape[0]
zs = grid[:,0,0]
ident = np.identity(NZ)
weight_func = interp1d(zs, ident, 'cubic')

每个网格只需要执行上述操作一次;你甚至可以重用weight_func，只要垂直坐标不变。
当需要插值时，weight_func(value)将为您提供权重，您可以使用该权重计算(x_idx, y_idx)处的单个插值：

weights = weight_func(value)
interp_val = np.dot(data[:, x_idx, y_idx), weights)

如果要计算插值的整个平面，可以使用np.inner，但由于z坐标在前，因此需要执行以下操作：

result = np.inner(data.T, weights).T

同样，计算应该实际上是即时的。

这是一个很老的问题，但现在最好的方法是使用MetPy的interpolate_1d函数：
https://unidata.github.io/MetPy/latest/api/generated/metpy.interpolate.interpolate_1d.html

在1、2和3维的规则网格上存在Numba加速插值的新实现：https://github.com/dbstein/fast_interp
用法如下：

from fast_interp import interp2d
import numpy as np

nx = 50
ny = 37
xv, xh = np.linspace(0, 1,       nx, endpoint=True,  retstep=True)
yv, yh = np.linspace(0, 2*np.pi, ny, endpoint=False, retstep=True)
x, y = np.meshgrid(xv, yv, indexing='ij')

test_function = lambda x, y: np.exp(x)*np.exp(np.sin(y))
f = test_function(x, y)
test_x = -xh/2.0
test_y = 271.43
fa = test_function(test_x, test_y)

interpolater = interp2d([0,0], [1,2*np.pi], [xh,yh], f, k=5, p=[False,True], e=[1,0])
fe = interpolater(test_x, test_y)