每次您拨打以下电话时，都会发生几件事情
scipy.interpolate.griddata：
- 首先，调用来
sp.spatial.qhull.Delaunay
对不规则的网格坐标进行三角剖分。 - 然后，对于新网格中的每个点，都将搜索三角剖分，以找到将其放置在哪个三角形中（实际上，在哪个单形中，在您的3D情况下哪个在四面体中）。
- 计算每个新网格点相对于封闭单形顶点的重心坐标。
- 使用重心坐标和该函数在封闭单形顶点处的值，为该网格点计算一个插值。
对于所有插值，前三个步骤都是相同的，因此，如果您可以为每个新的网格点存储封闭单形顶点的索引和插值权重，则可以最大程度地减少计算量。遗憾的是，要直接使用可用功能并不容易，尽管确实有可能：
import scipy.interpolate as spintimport scipy.spatial.qhull as qhullimport itertoolsdef interp_weights(xyz, uvw): tri = qhull.Delaunay(xyz) simplex = tri.find_simplex(uvw) vertices = np.take(tri.simplices, simplex, axis=0) temp = np.take(tri.transform, simplex, axis=0) delta = uvw - temp[:, d] bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta) return vertices, np.hstack((bary, 1 - bary.sum(axis=1, keepdims=True)))def interpolate(values, vtx, wts): return np.einsum('nj,nj->n', np.take(values, vtx), wts)该函数
interp_weights对上面列出的前三个步骤进行计算。然后，该函数
interpolate使用这些计算值非常快速地执行步骤4：
m, n, d = 3.5e4, 3e3, 3# make sure no new grid point is extrapolatedbounding_cube = np.array(list(itertools.product([0, 1], repeat=d)))xyz = np.vstack((bounding_cube, np.random.rand(m - len(bounding_cube), d)))f = np.random.rand(m)g = np.random.rand(m)uvw = np.random.rand(n, d)In [2]: vtx, wts = interp_weights(xyz, uvw)In [3]: np.allclose(interpolate(f, vtx, wts), spint.griddata(xyz, f, uvw))Out[3]: TrueIn [4]: %timeit spint.griddata(xyz, f, uvw)1 loops, best of 3: 2.81 s per loopIn [5]: %timeit interp_weights(xyz, uvw)1 loops, best of 3: 2.79 s per loopIn [6]: %timeit interpolate(f, vtx, wts)10000 loops, best of 3: 66.4 us per loopIn [7]: %timeit interpolate(g, vtx, wts)10000 loops, best of 3: 67 us per loop
所以首先，它的作用与相同
griddata，这很好。其次，设置插值，即计算
vtx和
wts与调用大致相同
griddata。但是第三，您现在几乎可以立即在同一网格上插值不同的值。
在
griddata此未想到的唯一事情是分配
fill_value给必须外推的点。您可以通过检查至少一个权重为负的点来做到这一点，例如：
def interpolate(values, vtx, wts, fill_value=np.nan): ret = np.einsum('nj,nj->n', np.take(values, vtx), wts) ret[np.any(wts < 0, axis=1)] = fill_value return ret
