问题描述
我需要在 4 个维度(纬度、经度、高度和时间)上线性插入温度数据.
点数相当多 (360x720x50x8),我需要一种快速计算数据范围内时空任意点温度的方法.
我尝试过使用 scipy.interpolate.LinearNDInterpolator
,但使用 Qhull 进行三角剖分在矩形网格上效率低下,需要数小时才能完成.
通过阅读这个 作为权重进行插值.
您可能不想弄乱重心坐标,因此最好将其留给 LinearNDInterpolator
.但是你确实知道一些关于三角剖分的事情.大多数情况下,因为你有一个规则的网格,在每个超立方体内,三角剖分将是相同的.因此,要插入单个值,您可以首先确定您的点在哪个子立方体中,使用该立方体的 16 个顶点构建一个 LinearNDInterpolator
,并使用它来插入您的值:
from itertools 导入产品定义插值器(坐标,数据,点):暗淡 = len(点)指数 = []sub_coords = []对于 xrange(dims) 中的 j :idx = np.digitize([point[j]], coords[j])[0]索引 += [[idx - 1, idx]]sub_coords += [coords[j][indices[-1]]]指数 = 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)返回 li([点])[0]>>>点 = np.array([12.3,-4.2, 500.5, 2.5])>>>插值器((纬度,经度,小数,时间),数据,点)0.386082399091
这不适用于矢量化数据,因为这需要为每个可能的子多维数据集存储一个 LinearNDInterpolator
,即使它可能比对整个事物进行三角剖分更快,它仍然会非常慢.
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.
I have tried using scipy.interpolate.LinearNDInterpolator
but using Qhull for triangulation is inefficient on a rectangular grid and takes hours to complete.
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).
Is there a quick way of interpolating data on a rectangular grid in 4 dimensions without it taking minutes to accomplish?
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.
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:
Using 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)
Using 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)
Thank you very much for your help!
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:
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.
Once LinearNDInterpolator
has been setup, there are two steps to coming up with an interpolated value for a given point:
- figure out inside which triangle (4D hypertetrahedron in your case) the point is, and
- interpolate using the barycentric coordinates of the point relative to the vertices as weights.
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
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 线性插值的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!