我的数据是一个普通的3D网格和作为一个最简单的例子我这里给出一个简单的立方体一个简单的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)
我不认为RegularGridInterpolator
有求导,因为它只能做线性或最近邻插值,所以它不是真正创造aything般衍生的引擎盖下任何有用的属性。你会需要某种正d多项式内插这似乎并不在scipy
据我可以看到实现的。当然,你总是可以只用np.gradient
你的栅格数据,然后让另一组iterpolators每个维度的。
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
还没有真正经检查发现,但应该是一般的想法。