如何解决这个三维规则网格插值问题

Question

我是一个新的python用户。我有一个h5文件，它是在一个固定的红移下引力势的快照。我已经读过python中的h5文件，现在我想写一个代码，通过三线性插值给出给定值（x，y，z）的重力势值。你们谁能帮我吗？为便于贵方参考，代码如下：

In [1]: import numpy as np

In [2]: import h5py

In [3]: from scipy.interpolate import RegularGridInterpolator

In [4]: f = h5py.File('my.h5', 'r')

In [5]: list(f.keys())
Out[5]: [u'data']

In [6]: data = f[u'data']

In [7]: data.shape
Out[7]: (64, 64, 64)

In [8]: data.dtype
Out[8]: dtype(('<f8', (3,)))

In [9]: data[0:63, 0:63, 0:63]
Out[9]: 
array([[[[ 7.44284016e-09, -3.69665900e-09,  8.75937447e-10],
         [ 8.00073078e-09, -2.62747161e-09,  9.82415717e-11],
         [ 7.81088465e-09, -2.03862452e-09, -4.00492778e-10],
         ...,
         [ 4.98376989e-09, -3.97621746e-09,  2.25554383e-09],
         [ 5.54899844e-09, -4.09876187e-09,  2.01146743e-09],
         [ 6.03652599e-09, -4.03159468e-09,  1.47328647e-09]],..............................

假设，我想用#RegularGridInterpolator函数求点（4.98376989e-09，-3.97621746e-09，2.25554383e-09）的电位值。我怎么能做到呢？在

Answer 1

这是一个很有趣（也很棘手）的问题，因此我决定用HDF5文件中的数据演示scipy interpolate示例是值得的。下面有两个代码部分。在

1. 第一种方法使用网格定义和 插值中使用的网格数据。

2. 第二个打开步骤1中的HDF5文件，并将x,y,z, mesh_data数据集作为示例中使用的Numpy数组读取。

运行此代码创建HDF5文件：

import numpy as np
import h5py

def f(x,y,z):
   return 2 * x**3 + 3 * y**2 - z

x = np.linspace(1, 4, 11)
y = np.linspace(4, 7, 22)
z = np.linspace(7, 9, 33)
mesh_data = f(*np.meshgrid(x, y, z, indexing='ij', sparse=True))

h5file = h5py.File('interpolate_test.h5', 'w')
h5file.create_dataset('/x', data=x)
h5file.create_dataset('/y', data=y)
h5file.create_dataset('/z', data=z)
h5file.create_dataset('/mesh_data', data=mesh_data)

h5file.close()

然后，运行此代码以读取带有h5py的HDF5文件并执行插值：

^{pr2}$

结果输出应该如下所示（与scipy示例相同）：

[125.80469388 146.30069388]

对于那些使用Pytables API读取HDF5数据的用户，下面是上面第2步的另一种方法。读取数据的过程相似，只是调用不同。在

运行此代码以读取带有Pytables的HDF5文件并执行插值：

import numpy as np
from scipy.interpolate import RegularGridInterpolator
import tables

h5file = tables.open_file('interpolate_test.h5')

x = h5file.root.x.read()
y = h5file.root.y.read()
z = h5file.root.z.read()
mesh_data = h5file.root.mesh_data.read()

my_interpolating_function = RegularGridInterpolator((x, y, z), mesh_data)

pts = np.array([[2.1, 6.2, 8.3], [3.3, 5.2, 7.1]])
print (my_interpolating_function(pts))

结果输出应与上述（以及scipy示例）相同：

[125.80469388 146.30069388]