我正在使用Numpy / Scipy尝试将三维大气数据从一个垂直坐标插值到另一个。例如,我有温度和相对湿度的立方体数据,两者都在恒定、规则的压力表面上。我想将相对湿度插值到恒定温度表面上。
先前曾提出过类似问题,位于这里,但那个解决方案非常慢。在我的情况下,我的立方体(30x321x321)约有3M个点,并且该方法需要大约4分钟才能处理一组数据。
那篇文章已经是近5年前的了。也许新版本的Numpy / Scipy有更快的方法来处理这个问题?也许新的思路可以更好地解决这个问题?我乐意听取建议。
编辑:
“Slow”表示每组数据需要4分钟。我不知道还有什么其他的量化方式。
正在使用的代码……
先前曾提出过类似问题,位于这里,但那个解决方案非常慢。在我的情况下,我的立方体(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转储文件。
mask = grid[0, ...] > grid[-1, ...]; grid[:, mask] = grid[::-1, mask]
(未经测试)。 - user707650scipy.ndimage.map_coordinates
(或者也许是最近添加的scipy.interpolate.RegularGridInterpolator
,两者基本上都是做同样的事情)。我现在没有时间给出完整的答案,但如果你的原始数据在一个规则的网格上,就不需要更昂贵的方法了。 - Joe Kington