python Numpy/Scipy环境下大气数据的快速三维插值

nlejzf6q  于 2023-02-07  发布在  Python
关注(0)|答案(3)|浏览(238)

我尝试使用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转储。

qij5mzcb

qij5mzcb1#

由于这是大气数据,我猜想您的网格没有均匀的间距;然而,如果你的网格是直线的(比如每个垂直列都有相同的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

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

t0ybt7op

t0ybt7op2#

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

wgx48brx

wgx48brx3#

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

相关问题