How Can I Solve This 3d Regular Grid Interpolation Problem

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

Solution 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]