How can I solve this 3D regular grid interpolation problem

Question

I am a new python user. I have a h5 file, which is a snapshot of gravitational potential at a fixed redshift. I have read the h5 file in python and now I want to write a code which will give the value of the gravitational potential for given values of (x, y, z) by using trilinear interpolation. Can anyone of you please help me to do that? For your kind consideration, the code is given below:

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]],..............................

Suppose, I want to find the value of potential at point (4.98376989e-09, -3.97621746e-09, 2.25554383e-09) by using #RegularGridInterpolator function. How can I do that?

Answer 1

This is an interesting (and tricky) enough question that I decided it deserved an answer to demonstrate the scipy interpolate example using data from an HDF5 file. There are 2 code sections below.

1. The first populates a HDF5 file with the mesh definition and mesh_data used in the interpolation.

2. The second opens the HDF5 file from Step 1 and reads the x,y,z, mesh_data datasets as Numpy arrays used in the example.

Run this code to create the HDF5 file:

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()

Then, run this code to read the HDF5 file with h5py and do the interpolation:

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

Resulting output should look like this (which is identical to the scipy example):

[125.80469388 146.30069388]

For those that use the Pytables API to read HDF5 data, here is an alternative method for Step 2 above. The process to read the data is similar, only the calls are different.

Run this code to read the HDF5 file with Pytables and do the interpolation:

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

Resulting output should be identical to above (and to the scipy example):

[125.80469388 146.30069388]