python - 使用 scipy.ndimage.map_cooperatives 对 ndarray 中的所有数组进行高效采样

Question

我有一个 3D 屏蔽数组堆栈。我想在相同的固定位置对堆栈中的所有数组进行采样。

stack.ma_stack.shape

(1461、390、327)

#Indices to be sampled
x = np.array([ 117.38670304,  119.1220485 ])
y = np.array([ 209.98120554,  210.37202372])

以下方法非常有效，但仅适用于整数索引:

x_int = np.rint(x).astype(int)
y_int = np.rint(y).astype(int)
samp = stack.ma_stack[:,y_int,x_int]
samp.shape

(1461, 2)

我正在尝试为浮点索引实现 scipy.ndimage.map_coordinates 插值采样，但我似乎无法弄清楚如何正确格式化坐标。

大多数示例使用 map_coordinates 对单个数组进行采样，以下适用于堆栈中的单个数组:

map_coord = np.array([[y,], [x,]])
samp = scipy.ndimage.map_coordinates(stack.ma_stack[0], map_coord, order=1)
samp.shape

(1, 2)

我可以循环遍历堆栈中的每个数组，但我知道有一个简单的索引技巧可以在一次调用中对整个堆栈进行采样。我阅读了有关 mgrid 的内容，并做了一些实验，但找不到正确的解决方案(我仍在学习高级索引)。我知道有人会立即知道答案。谢谢。

相关说明:有人知道如何对屏蔽数组执行此操作，而不用 fill_value 或 np.nan 替换丢失的数据吗？ ndimage 插值不能很好地处理屏蔽数组:

https://github.com/scipy/scipy/issues/1682

Answer 1

必须有一种方法让它自动广播...同时，您可以使用 np.arange(...) 强制广播以从每个二维数组中获取一个点在堆栈中:

map_coords = np.broadcast_arrays(np.arange(stack.ma_stack.shape[0])[:, None], y, x)
samp = ndimage.map_coordinates(stack.ma_stack, map_coords, order=1)

虽然效率很低，因为“广播”是提前完成的(并且可能复制所有数据)，但它仍然比循环快一点:

In [88]: a = np.random.rand(1461, 390, 327)

In [89]: x = np.array([ 117.38670304,  119.1220485 ])

In [90]: y = np.array([ 209.98120554,  210.37202372])

In [107]: %%timeit
   .....: map_coord = np.array([[y,], [x,]])
   .....: np.concatenate([ndimage.map_coordinates(ai, map_coord, order=1) for ai in a])
   .....: 
10 loops, best of 3: 33.1 ms per loop

In [108]: %%timeit
   .....: map_coords = np.broadcast_arrays(np.arange(a.shape[0])[:, None], y, x)
   .....: ndimage.map_coordinates(a, map_coords, order=1)
   .....: 
100 loops, best of 3: 4.67 ms per loop

In [109]: samp_OP = np.concatenate([ndimage.map_coordinates(ai, map_coord, order=1) for ai in a])

In [110]: samp_chan = ndimage.map_coordinates(a, map_coords, order=1)

In [111]: np.allclose(samp_chan, samp_OP)
Out[111]: True