numpy 3D meshgrid 仅作为视图答案

【问题标题】：numpy 3D meshgrid only as a viewnumpy 3D meshgrid 仅作为视图
【发布时间】：2018-12-28 00:43:48
【问题描述】：

我有一个由间距 xi,yi,zi 定义的立方网格：

xi,yi,zi = [linspace(ox,ox+s*d,s) for ox,s,d in zip(origin,size,delta)]

我还在该网格上设置了一组标量值W。 W.shape() == size。我想使用scipy's linear interpolation，它需要作为输入：

类scipy.interpolate.LinearNDInterpolator(points, values):

参数：

points : ndarray 浮点数，形状 (npoints, ndims) 数据点坐标。

values ：浮点数或复数的ndarray，形状(npoints, ...) 数据值。

如何从xi,yi,zi 创建一组假的points（通过魔法广播）？现在我正在创建一个中间数组来提供给插值函数 - 有没有更好的方法？

相关问题：Numpy meshgrid in 3D。这篇文章中的答案实际上创建了网格 - 我只想将其模拟为另一个函数的输入（首选纯 numpy 解决方案）。

【问题讨论】：

您想要 Nxdim 数组但欺骗 numpy 使其实际上不分配完整数组？这是不可能的。您必须使用为常规网格设计的工具，但我想这对于 scipy 中的更高维度不存在。
@Sebastian，您可以从较小的数组中模拟较大的数组。例如，如果x.shape,y.shape = (n,m)，您可以使用[X,Y] = np.meshgrid(x,y); S=X+Y 而不是S=x+y[:,np.newaxis]，创建f(x,y) = x+y not 的广播数组。有关详细信息，请参阅链接的问题。
是的，但您似乎想要模拟更多元素，然后实际上是在 xi,zi,yi 数组中。这在理论上是可能的（使用 stride_tricks）。但是，由于生成的数组应该是二维的，因此无法在这种情况下构造它，即使不太可能 scipy 也不会在以后创建一个副本。
@Sebastian xi,yi,zi -> (n,m,l) 的形状表明我想要一个形状为(n*m*l,3) 的对象points。这正是meshgrid 对二维xi,yi -> (n,m) 到(n*m,2) 所做的。在某些时候，scipy 无论如何都会为插值创建一个内部副本，但我试图在这一步避免不必要的副本。
@Hooked，这不是meshgrid 所做的。你的标题有点误导。我已经给出了一个答案，它可以按照您的要求制作网格，但我认为您不能保持这种方式并将其变成某种形状 (n*m*l,3)。

标签： python numpy scipy array-broadcasting

【解决方案1】：

>>> xi, yi, zi = [np.arange(3) for i in range(3)]
>>> xx, yy, zz = np.broadcast_arrays(xi,yi[:,np.newaxis],zi[:,np.newaxis,np.newaxis])
>>> xx.shape
(3, 3, 3)
>>> xx.strides
(0, 0, 8)

您可以看到它没有创建新副本，因为前两个维度的步幅为 0。

我也写了一个 n 维版本：

def ndmesh(*args):
   args = map(np.asarray,args)
   return np.broadcast_arrays(*[x[(slice(None),)+(None,)*i] for i, x in enumerate(args)])

【讨论】：

避免 (mnl,3) 点数组的任何额外临时数组（不是视图）的唯一方法是使用 xx.flat（而不是 ravel 等）插入到未初始化的点数组中。
拨打ndmesh时；论据是什么？
要广播的数组。在我给出的示例中，它是ndmesh(xi,yi,zi)。

【解决方案2】：

您可以按照其他答案中解释的类似方式构造必要的points 数组：

xx, yy, zz = np.broadcast_arrays(xi[:,None,None], yi[None,:,None], zi[None,None,:])
points = (xx.ravel(), yy.ravel(), zz.ravel())
ip = LinearNDInterpolator(points, data.ravel())

但是，如果您有一个规则网格，那么使用LinearNDInterpolator 很可能不是最佳选择，因为它是为分散数据插值而设计的。它构建了数据点的 Delaunay 三角剖分，但在这种情况下，原始数据已经具有非常规则的结构，可以更有效地利用。

由于您的网格是矩形的，您可以将插值构建为三个一维插值的张量积。 Scipy 没有这个内置的（到目前为止），但它很容易做到，见这个线程：http://mail.scipy.org/pipermail/scipy-user/2012-June/032314.html （使用例如 interp1d 而不是 pchip 来获得一维插值）

【讨论】：

【解决方案3】：

我不相信有任何方法可以将缺少完整副本的东西传递给 LinearNDInterpolator（因为也没有用于三个维度的常规网格的函数）。所以唯一避免创建完整数组的地方是在创建这个点数组期间，我不知道你现在是怎么做的，所以在这方面它可能已经很有效了，但我想它可能不值得麻烦避免这个。

其他然后 np.mgrid+reshape 也许这样的东西可能是一个选项（对于 n 维也不难写）：

# Create broadcastest versions of xi, yi and zi
# np.broadcast_arrays does not allocate the full arrays
xi, yi, zi = np.broadcast_arrays(xi[:,None,None], yi[:,None,None], zi[:,None,None])

# then you could use .flat to fill a point array:
points = np.empty((xi.size, 3), dtype=xi.dtype)
points[:,0] = xi.flat
points[:,1] = yi.flat
points[:,2] = zi.flat

与.repeat函数相反，这里创建的临时数组不比原来的xi等数组大。

【讨论】：