矩形网格上的Python 4D线性插值

Question

我需要在4个维度(纬度,经度,高度和时间)中线性插值温度数据.
点数相当高(360x720x50x8),我需要一种快速的方法来计算数据范围内空间和时间任意点的温度.

我尝试过使用scipy.interpolate.LinearNDInterpolator但是使用Qhull进行三角测量在矩形网格上效率低下并且需要数小时才能完成.

通过阅读此SciPy票证,解决方案似乎是使用标准实现新的nd插值器interp1d来计算更多数据点,然后使用"最近邻居"方法和新数据集.

然而,这需要很长时间(分钟).

有没有一种快速的方法可以在4维的矩形网格上插入数据而无需花费几分钟才能完成？

我想过使用interp1d4次而不计算更高密度的点,但留给用户用坐标调用,但我无法理解如何做到这一点.

否则我会根据自己的需要编写自己的4D内插器吗？

这是我用来测试这个的代码:

使用scipy.interpolate.LinearNDInterpolator:

import numpy as np
from scipy.interpolate import LinearNDInterpolator

lats = np.arange(-90,90.5,0.5)
lons = np.arange(-180,180,0.5)
alts = np.arange(1,1000,21.717)
time = np.arange(8)
data = np.random.rand(len(lats)*len(lons)*len(alts)*len(time)).reshape((len(lats),len(lons),len(alts),len(time)))

coords = np.zeros((len(lats),len(lons),len(alts),len(time),4))
coords[...,0] = lats.reshape((len(lats),1,1,1))
coords[...,1] = lons.reshape((1,len(lons),1,1))
coords[...,2] = alts.reshape((1,1,len(alts),1))
coords[...,3] = time.reshape((1,1,1,len(time)))
coords = coords.reshape((data.size,4))

interpolatedData = LinearNDInterpolator(coords,data)

使用scipy.interpolate.interp1d:

import numpy as np
from scipy.interpolate import LinearNDInterpolator

lats = np.arange(-90,90.5,0.5)
lons = np.arange(-180,180,0.5)
alts = np.arange(1,1000,21.717)
time = np.arange(8)
data = np.random.rand(len(lats)*len(lons)*len(alts)*len(time)).reshape((len(lats),len(lons),len(alts),len(time)))

interpolatedData = np.array([None, None, None, None])
interpolatedData[0] = interp1d(lats,data,axis=0)
interpolatedData[1] = interp1d(lons,data,axis=1)
interpolatedData[2] = interp1d(alts,data,axis=2)
interpolatedData[3] = interp1d(time,data,axis=3)

非常感谢您的帮助!

Answer 1

在您链接的同一张票中,有一个他们称之为张量积插值的示例实现,显示了嵌套递归调用的正确方法interp1d.如果您选择默认kind='linear'参数,这相当于quadrilinear插值interp1d.

虽然这可能足够好,但这不是线性插值,并且插值函数中会有更高阶的项,因为来自维基百科条目的双线性插值的图像显示:

对于你所追求的东西来说,这可能已经足够好了,但是有些应用程序需要三角形,真正的分段线性插值.如果你真的需要这个,有一种简单的方法可以解决qhull的缓慢问题.

一旦LinearNDInterpolator已设置,有两个步骤来想出给定点的插补值:

1. 弄清楚哪个三角形(在你的情况下是4D高三面体)的点是,和
2. 使用相对于顶点的点的重心坐标作为权重进行插值.

你可能不想搞乱重心坐标,所以最好留下它LinearNDInterpolator.但你确实知道有关三角测量的一些事情.大多数情况下,因为你有一个规则的网格,在每个超立方体内,三角测量将是相同的.因此,要插入单个值,您可以首先确定您的点位于哪个子多维数据LinearNDInterpolator集中,使用该多维数据集的16个顶点构建a ,并使用它来插值:

from itertools import product

def interpolator(coords, data, point) :
    dims = len(point)
    indices = []
    sub_coords = []
    for j in xrange(dims) :
        idx = np.digitize([point[j]], coords[j])[0]
        indices += [[idx - 1, idx]]
        sub_coords += [coords[j][indices[-1]]]
    indices = np.array([j for j in product(*indices)])
    sub_coords = np.array([j for j in product(*sub_coords)])
    sub_data = data[list(np.swapaxes(indices, 0, 1))]
    li = LinearNDInterpolator(sub_coords, sub_data)
    return li([point])[0]

>>> point = np.array([12.3,-4.2, 500.5, 2.5])
>>> interpolator((lats, lons, alts, time), data, point)
0.386082399091

这不适用于矢量化数据,因为这需要LinearNDInterpolator为每个可能的子立方体存储一个,即使它可能比对整个事物进行三角测量更快,但它仍然会非常慢.

Answer 2

scipy.ndimage.map_coordinates 是一个很好的快速插值器，用于统一网格（所有框的大小相同）。有关清晰的描述，请参阅SO 上的multivariate-spline-interpolation-in-python-scipy。

对于非均匀矩形网格，一个简单的包装器 Intergrid将非均匀网格映射/缩放到均匀网格，然后执行 map_coordinates。在像您这样的 4d 测试用例中，每次查询大约需要 1 微秒：

Intergrid: 1000000 points in a (361, 720, 47, 8) grid took 652 msec