python - Numpy/Scipy 中大气数据的快速 3D 插值

Question

我正在尝试使用 Numpy/Scipy 将 3D 大气数据从一个垂直坐标插入到另一个垂直坐标。例如，我有温度和相对湿度的立方体，它们都位于恒定、规则的压力表面上。我想将相对湿度内插到恒温表面。

我试图解决的确切问题之前已经被问过 here但是，那里的解决方案非常慢。就我而言，我的立方体 (30x321x321) 中有大约 3M 个点，该方法需要大约 4 分钟来处理一组数据。

该帖子已有近 5 年历史。较新版本的 Numpy/Scipy 是否有处理速度更快的方法？也许看待问题的新眼光有更好的方法？我乐于接受建议。

编辑:

一组数据立方体慢 = 4 分钟。我不确定我还能如何量化它。

正在使用的代码...

def interpLevel(grid,value,data,interp='linear'):
    """
    Interpolate 3d data to a common z coordinate.

    Can be used to calculate the wind/pv/whatsoever values for a common
    potential temperature / pressure level.

    grid : numpy.ndarray
       The grid. For example the potential temperature values for the whole 3d
       grid.

    value : float
       The common value in the grid, to which the data shall be interpolated.
       For example, 350.0

    data : numpy.ndarray
       The data which shall be interpolated. For example, the PV values for
       the whole 3d grid.

    kind : str
       This indicates which kind of interpolation will be done. It is directly
       passed on to scipy.interpolate.interp1d().

    returns : numpy.ndarray
       A 2d array containing the *data* values at *value*.

    """
    ret = np.zeros_like(data[0,:,:])
    for yIdx in xrange(grid.shape[1]):
        for xIdx in xrange(grid.shape[2]):
            # check if we need to flip the column
            if grid[0,yIdx,xIdx] > grid[-1,yIdx,xIdx]:
                ind = -1
            else:
                ind = 1
            f = interpolate.interp1d(grid[::ind,yIdx,xIdx], \
                data[::ind,yIdx,xIdx], \
                kind=interp)
            ret[yIdx,xIdx] = f(value)
    return ret

编辑 2:我可以分享示例数据的 npy 转储，如果有人有足够的兴趣看看我正在使用什么。

Answer 1

因为这是大气数据，我想你的网格没有统一的间距；但是，如果您的网格是直线的(这样每个垂直列都具有相同的 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, :, :]

(当然，您需要添加检查 value 是否超出范围)。对于像您这样大小的网格，这将非常快(如几分之一秒)

如果需要，您可以很容易地修改上面的内容以执行三次插值；挑战在于为不均匀的垂直间距选择正确的权重。

使用 scipy.ndimage.map_coordinates 的问题在于，尽管它提供了更高阶的插值并且可以处理任意样本点，但它确实假定输入数据是均匀间隔的。它仍会产生平滑的结果，但不是可靠的近似值。

如果您的坐标网格不是直线的，因此给定索引的 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

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