矩形网格上的Python 4D线性插值 [英] Python 4D linear interpolation on a rectangular grid

问题描述

我需要在4个维度(纬度，经度，海拔和时间)上线性插值温度数据.
点的数量相当多(360x720x50x8)，我需要一种快速的方法来计算数据范围内在空间和时间上任何一点的温度.

I need to interpolate temperature data linearly in 4 dimensions (latitude, longitude, altitude and time).
The number of points is fairly high (360x720x50x8) and I need a fast method of computing the temperature at any point in space and time within the data bounds.

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

I have tried using scipy.interpolate.LinearNDInterpolator but using Qhull for triangulation is inefficient on a rectangular grid and takes hours to complete.

通过阅读 SciPy票证，该解决方案似乎正在实现新的解决方案内插器使用标准interp1d来计算更多的数据点，然后对新数据集使用最近邻居"方法.

By reading this SciPy ticket, the solution seemed to be implementing a new nd interpolator using the standard interp1d to calculate a higher number of data points, and then use a "nearest neighbor" approach with the new dataset.

但是，这又花费了很长时间(几分钟).

This, however, takes a long time again (minutes).

是否有一种快速的方法可以在4维的矩形网格上插值数据而无需花费几分钟的时间来完成?

Is there a quick way of interpolating data on a rectangular grid in 4 dimensions without it taking minutes to accomplish?

我考虑过使用interp1d 4次而没有来计算更高的点密度，但是让用户使用坐标来调用它，但是我无法理解如何做到这一点.

I thought of using interp1d 4 times without calculating a higher density of points, but leaving it for the user to call with the coordinates, but I can't get my head around how to do this.

否则，可以根据自己的需要编写自己的4D插值器吗?

Otherwise would writing my own 4D interpolator specific to my needs be an option here?

这是我一直在测试的代码:

Here's the code I've been using to test this:

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

非常感谢您的帮助！

推荐答案

在您链接的同一张票证中，有一个示例性实施例，称其为张量积插值，显示了正确的方法将对interp1d的递归调用嵌套.如果为interp1d选择默认的kind='linear'参数，则等效于四线性插值.

In the same ticket you have linked, there is an example implementation of what they call tensor product interpolation, showing the proper way to nest recursive calls to interp1d. This is equivalent to quadrilinear interpolation if you choose the default kind='linear' parameter for your interp1d's.

虽然这可能足够好，但这不是线性插值，并且插值函数中将存在更高阶的项，如

While this may be good enough, this is not linear interpolation, and there will be higher order terms in the interpolation function, as this image from the wikipedia entry on bilinear interpolation shows:

这对于您所追求的东西可能已经足够好了，但是在某些应用中，首选三角剖分，真正分段线性插值的应用程序.如果您确实需要此功能，可以通过一种简单的方法来解决qhull的缓慢性.

This may very well be good enough for what you are after, but there are applications where a triangulated, really piecewise linear, interpoaltion is preferred. If you really need this, there is an easy way of working around the slowness of qhull.

设置好LinearNDInterpolator后，需要两个步骤来得出给定点的插值:

Once LinearNDInterpolator has been setup, there are two steps to coming up with an interpolated value for a given point:

1. 弄清楚该点在哪个三角形内(在您的情况下为4D超四面体)，并且
2. 使用相对于顶点的点的重心坐标进行权重.

1. figure out inside which triangle (4D hypertetrahedron in your case) the point is, and
2. interpolate using the barycentric coordinates of the point relative to the vertices as weights.

您可能不想弄乱重心坐标，因此最好将其留给LinearNDInterpolator.但是您确实知道有关三角剖分的一些知识.通常，由于您有规则的网格，因此在每个超立方体中，三角剖分将是相同的.因此，要内插一个值，您可以首先确定您的点在哪个子多维数据集中，用该多维数据集的16个顶点构建一个LinearNDInterpolator，并使用它来内插您的值:

You probably do not want to mess with barycentric coordinates, so better leave that to LinearNDInterpolator. But you do know some things about the triangulation. Mostly that, because you have a regular grid, within each hypercube the triangulation is going to be the same. So to interpolate a single value, you could first determine in which subcube your point is, build a LinearNDInterpolator with the 16 vertices of that cube, and use it to interpolate your value:

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，尽管它可能比对整个对象进行三角剖分要快，但仍然非常慢.

This cannot work on vectorized data, since that would require storing a LinearNDInterpolator for every possible subcube, and even though it probably would be faster than triangulating the whole thing, it would still be very slow.

这篇关于矩形网格上的Python 4D线性插值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！