查找连续的未屏蔽值 [英] Find consecutive unmasked values

问题描述

我有一个很大的3维(时间，经度，纬度)输入数组.大多数条目都被屏蔽.我需要找到那些掩码为False的条目的时间长于特定数量的连续时间步长(在此将其称为threshold).结果应该是与输入蒙版具有相同形状的蒙版.

I have a large 3 dimensional (time, longitude, latitude) input array. Most of the entries are masked. I need to find those entries where the mask is False for longer than a specific number of consecutive time steps (which I call threshold here). The result should be a mask with the same shape as the input mask.

以下是一些伪代码，希望可以使我的意思更清楚:

Here is some pseudo-code to hopefully make clearer what I mean:

new_mask = find_consecutive(mask, threshold=3)
mask[:, i_lon, i_lat]
# [1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0]
new_mask[:, i_lon, i_lat]
# [1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]

我不确定到目前为止我的方法是否有意义.它在性能方面做得很好，并为我提供了一个带标签的数组，并提供了有关所需标签的知识.我只是想不出一种有效的方法来再次将labels转换为蒙版.

I'm not sure if my approach so far make sense. It does good performance-wise and gives me a labeled array and knowledge about which labels I want. I just couldn't figure out an efficient way to transform labels into a mask again.

from scipy.ndimage import measurements
structure = np.zeros((3, 3, 3))
structure[:, 1, 1] = 1
labels, nr_labels = measurements.label(1 - mask, structure=structure)
_, counts = np.unique(labels, return_counts=True)
labels_selected = [i_count for i_count, count in enumerate(counts)
                   if count >= threshold]

推荐答案

这是 binary closing operation in image-processing .要解决此问题，您可以从scipy模块寻求帮助，特别是-ONES和长度threshold的适当一维内核后，rel ="nofollow"> scipy.ndimage.morphology.binary_closing .另外，Scipy的binary closing函数仅给我们提供了封闭的遮罩.因此，要获得所需的输出，我们需要使用输入掩码OR对其进行输出.因此，实现看起来像这样-

That's a classical case of binary closing operation in image-processing. To solve it you can take help from scipy module, specifically - scipy.ndimage.morphology.binary_closing after we feed an appropriate 1D kernel of all ONES and of length threshold. Also, Scipy's binary closing function gives us the closed mask only. So, to get the desired output, we need to OR it with the input mask. Thus, the implementation would look something like this -

from scipy.ndimage import binary_closing
out = mask | binary_closing(mask, structure=np.ones(threshold))

关于二进制关闭的NumPy版本！

现在，关闭操作基本上是 image-dilation和image-erosion ，因此我们可以使用可信赖的卷积运算来模拟该行为，并且在NumPy中将其作为

Now, a closing operation is basically image-dilation and image-erosion, so we can simulate that behiviour using the trusty convolution operation and we do have that here in NumPy as np.convolve. Similar to scipy's binary closing operation, we need the same kernel here as well and we would use it both for dilation and erosion. The implementation would be -

def numpy_binary_closing(mask,threshold):

    # Define kernel
    K = np.ones(threshold)

    # Perform dilation and threshold at 1
    dil = np.convolve(mask,K,mode='same')>=1

    # Perform erosion on the dilated mask array and threshold at given threshold
    dil_erd = np.convolve(dil,K,mode='same')>= threshold
    return dil_erd

样品运行-

In [133]: mask
Out[133]: 
array([ True, False, False, False, False,  True,  True, False, False,
        True, False], dtype=bool)

In [134]: threshold = 3

In [135]: binary_closing(mask, structure=np.ones(threshold))
Out[135]: 
array([False, False, False, False, False,  True,  True,  True,  True,
        True, False], dtype=bool)

In [136]: numpy_binary_closing(mask,threshold)
Out[136]: 
array([False, False, False, False, False,  True,  True,  True,  True,
        True, False], dtype=bool)

In [137]: mask | binary_closing(mask, structure=np.ones(threshold))
Out[137]: 
array([ True, False, False, False, False,  True,  True,  True,  True,
        True, False], dtype=bool)

In [138]: mask| numpy_binary_closing(mask,threshold)
Out[138]: 
array([ True, False, False, False, False,  True,  True,  True,  True,
        True, False], dtype=bool)

运行时测试(Scipy与Numpy！)

案例1:均匀稀疏

In [163]: mask = np.random.rand(10000) > 0.5

In [164]: threshold = 3

In [165]: %timeit binary_closing(mask, structure=np.ones(threshold))
1000 loops, best of 3: 582 µs per loop

In [166]: %timeit numpy_binary_closing(mask,threshold)
10000 loops, best of 3: 178 µs per loop

In [167]: out1 = binary_closing(mask, structure=np.ones(threshold))

In [168]: out2 = numpy_binary_closing(mask,threshold)

In [169]: np.allclose(out1,out2) # Verify outputs
Out[169]: True

案例2:稀疏且阈值更大

Case #2 : More sparse and bigger threshold

In [176]: mask = np.random.rand(10000) > 0.8

In [177]: threshold = 11

In [178]: %timeit binary_closing(mask, structure=np.ones(threshold))
1000 loops, best of 3: 823 µs per loop

In [179]: %timeit numpy_binary_closing(mask,threshold)
1000 loops, best of 3: 331 µs per loop

In [180]: out1 = binary_closing(mask, structure=np.ones(threshold))

In [181]: out2 = numpy_binary_closing(mask,threshold)

In [182]: np.allclose(out1,out2) # Verify outputs
Out[182]: True

优胜者是 Numpy ，而且赢利很大！

Winner is Numpy and by a big margin!

边界条件

如果1s足够接近Boemdries，似乎边界也需要关闭.要解决这些情况，可以在输入布尔数组的开头和结尾分别填充一个1，使用发布的代码，然后在末尾取消选择第一个和最后一个元素.因此，使用scipy的binary_closing方法的完整实现将是-

It seems the boundaries need the closing too, if the 1s are close enough to the bounadries. To solve those cases, you can pad one 1 each at the start and end of the input boolean array, use the posted code and then at the end de-select the first and last element. Thus, the complete implementation using scipy's binary_closing approach would be -

mask_ext = np.pad(mask,1,'constant',constant_values=(1))
out = mask_ext | binary_closing(mask_ext, structure=np.ones(threshold))
out = out[1:-1]

样品运行-

In [369]: mask
Out[369]: 
array([False, False,  True, False, False, False, False,  True,  True,
       False, False,  True, False], dtype=bool)

In [370]: threshold = 3

In [371]: mask_ext = np.pad(mask,1,'constant',constant_values=(1))
     ...: out = mask_ext | binary_closing(mask_ext, structure=np.ones(threshold))
     ...: out = out[1:-1]
     ...: 

In [372]: out
Out[372]: 
array([ True,  True,  True, False, False, False, False,  True,  True,
        True,  True,  True,  True], dtype=bool)

这篇关于查找连续的未屏蔽值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！