如何解决这个 3D 规则网格插值问题

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 插值示例。下面有 2 个代码部分。

1. 第一个使用网格定义填充 HDF5 文件，然后 插值中使用的 mesh_data。

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 文件并进行插值：

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

h5file = h5py.File('interpolate_test.h5')

x = h5file['x'][:]
y = h5file['y'][:]
z = h5file['z'][:]
mesh_data = h5file['mesh_data'][:,:,:]

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]

对于那些使用 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]