Как я могу решить эту проблему трехмерной регулярной интерполяции сетки

Question

Я новый пользователь Python.У меня есть файл h5, который представляет собой снимок гравитационного потенциала при фиксированном красном смещении.Я прочитал файл h5 на python, и теперь я хочу написать код, который даст значение гравитационного потенциала для заданных значений (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]],..............................

Предположим, я хочу найти значение потенциала в точке (4.98376989e-09, -3.97621746e-09, 2.25554383e-09) поиспользуя функцию #RegularGridInterpolator.Как я могу это сделать?

Answer 1

Это достаточно интересный (и хитрый) вопрос, поэтому я решил, что заслуживает ответа, чтобы продемонстрировать пример скучной интерполяции с использованием данных из файла HDF5.Ниже есть 2 раздела кода.

1. Первый заполняет файл HDF5 определением сетки и mesh_data, используемыми при интерполяции.

2. Второй открывает файл HDF5 с шага 1и читает наборы данных 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()

Затем запустите этот код, чтобы прочитать файл HDF5 с помощью h5py и выполнить интерполяцию:

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 выше.Процесс чтения данных похож, только звонки разные.

Запустите этот код для чтения файла HDF5 с помощью Pytables и выполните интерполяцию:

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]