Вычислить производную данных трехмерной регулярной сетки после интерполяции - PullRequest
Новая
       33

Вычислить производную данных трехмерной регулярной сетки после интерполяции

0 голосов
/ 04 февраля 2019

Мои данные - это обычная трехмерная сетка, и в качестве самого простого примера я привел здесь код простого файла h5 для простого куба, а для интерполяции я использую функцию RegularGridInterpolator.Но я хочу знать, как я могу изменить свой код, чтобы я мог получить производную от этой интерполированной функции?Для вашей доброй информации, я дал здесь мой код:

код для генерации файла h5: 

import numpy as np
import h5py

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

x = np.linspace(-1, 1, 2)
y = np.linspace(-1, 1, 2)
z = np.linspace(-1, 1, 2)
mesh_data = f(*np.meshgrid(x, y, z, indexing='ij', sparse=True))

h5file = h5py.File('cube.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()

Код для интерполяции: 

import numpy as np   
import h5py
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import RegularGridInterpolator
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
f = h5py.File('cube.h5', 'r')  #importing h5 file from gevolution output
list(f.keys())
dset = f[u'mesh_data']
dset.shape
dset.value.shape
dset[0:2,0:2,0:2]
x = np.linspace(-1, 1, 2)
y = np.linspace(-1, 1, 2)
z = np.linspace(-1, 1, 2)
def deriv_interp(x, y, z, dset):
    d_dset = np.gradient(dset)
    deriv_interpolators = [RegularGridInterpolator((x, y, z), d, method = 'linear') for d in d_dset]
    def interp(x,y,z):
        return np.dstack([d(x,y,z) for d in deriv_interpolators])
    return interp
pts = np.array([0.2, 0.9, 0.6]) 
deriv_interpolators(pts)

1 Ответ

0 голосов
/ 04 февраля 2019

Я не думаю, что RegularGridInterpolator имеет какой-либо полезный атрибут для поиска производной, поскольку он может выполнять только линейную интерполяцию или интерполяцию ближайших соседей, поэтому он не создает ничего похожего на производную под капотом.Вам понадобится какой-то полиномиальный интерполятор, который, как я вижу, не реализован в scipy.Конечно, вы всегда можете просто использовать np.gradient для данных в сетке, а затем создать еще один набор итерполяторов для каждого измерения. 

def deriv_interp(x, y, z, dset):
    d_dset = np.gradient(dset)
    deriv_interpolators = [RegularGridInterpolator((x, y, z), d, method = 'linear') for d in d_dset]
    def interp(coord):
        return np.dstack([d(coord) for d in deriv_interpolators])
    return interp

На самом деле это не проверялось, но этодолжна быть общая идея.

...