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

Question

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

pythonnumpyscipyinterpolation

4

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

先前曾提出过类似问题，位于这里，但那个解决方案非常慢。在我的情况下，我的立方体（30x321x321）约有3M个点，并且该方法需要大约4分钟才能处理一组数据。

那篇文章已经是近5年前的了。也许新版本的Numpy / Scipy有更快的方法来处理这个问题？也许新的思路可以更好地解决这个问题？我乐意听取建议。

编辑：

“Slow”表示每组数据需要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转储文件。

- wedgef5

您IP地址为143.198.54.68，由于运营成本限制，当前对于免费用户的使用频率限制为每个IP每72小时10次对话，如需解除限制，请点击左下角设置图标按钮（手机用户先点击左上角菜单按钮）。 - user559633

你可以尝试将列翻转移出双重循环，使用numpy的掩码。类似这样：mask = grid[0, ...] > grid[-1, ...]; grid[:, mask] = grid[::-1, mask]（未经测试）。 - user707650

1

请查看numba。 - elyase

听起来你想要使用 scipy.ndimage.map_coordinates（或者也许是最近添加的 scipy.interpolate.RegularGridInterpolator，两者基本上都是做同样的事情）。我现在没有时间给出完整的答案，但如果你的原始数据在一个规则的网格上，就不需要更昂贵的方法了。 - Joe Kington

3个回答

1

这是一个比较老的问题，但如今最好的方法是使用 MetPy 的 interpolate_1d 函数:

https://unidata.github.io/MetPy/latest/api/generated/metpy.interpolate.interpolate_1d.html

- Fabzi

0

在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)

- Ales

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Gretchen · Accepted Answer

考虑到这是大气数据，我猜测你的网格并不具有均匀间距；但是，如果你的网格是规则的（每个垂直列具有相同的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, :, :]

（当然，您需要添加对值超出范围的检查）。对于您的大小网格，这将非常快（小数秒级）。

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

使用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

再次强调，计算应该是非常快的。