在Python numpy掩码数组中用最近邻填充缺失值？

Question

我正在用Python处理一个二维的Numpy掩码数组。

我需要把被掩盖的区域的数据值改成离它最近的未被掩盖的值。

注意：如果有多个离它最近的未被掩盖的值，随便取一个就行（选哪个都可以，只要编码起来简单就好…）

比如：

import numpy
import numpy.ma as ma

a = numpy.arange(100).reshape(10,10)
fill_value=-99
a[2:4,3:8] = fill_value
a[8,8] = fill_value
a = ma.masked_array(a,a==fill_value)

>>> a  [[0 1 2 3 4 5 6 7 8 9]
  [10 11 12 13 14 15 16 17 18 19]
  [20 21 22 -- -- -- -- -- 28 29]
  [30 31 32 -- -- -- -- -- 38 39]
  [40 41 42 43 44 45 46 47 48 49]
  [50 51 52 53 54 55 56 57 58 59]
  [60 61 62 63 64 65 66 67 68 69]
  [70 71 72 73 74 75 76 77 78 79]
  [80 81 82 83 84 85 86 87 -- 89]
  [90 91 92 93 94 95 96 97 98 99]],

• 我希望它看起来像这样：

>>> a.data
 [[0 1 2 3 4 5 6 7 8 9]
 [10 11 12 13 14 15 16 17 18 19]
 [20 21 22 ? 14 15 16 ? 28 29]
 [30 31 32 ? 44 45 46 ? 38 39]
 [40 41 42 43 44 45 46 47 48 49]
 [50 51 52 53 54 55 56 57 58 59]
 [60 61 62 63 64 65 66 67 68 69]
 [70 71 72 73 74 75 76 77 78 79]
 [80 81 82 83 84 85 86 87 ? 89]
 [90 91 92 93 94 95 96 97 98 99]],

注意：“？”可以取任何一个相邻的未被掩盖的值。

有没有什么高效的方法来做到这一点？

谢谢你的帮助。

Answer 1

对于更复杂的情况，你可以使用scipy.spatial这个库：

from scipy.spatial import KDTree
x,y=np.mgrid[0:a.shape[0],0:a.shape[1]]

xygood = np.array((x[~a.mask],y[~a.mask])).T
xybad = np.array((x[a.mask],y[a.mask])).T

a[a.mask] = a[~a.mask][KDTree(xygood).query(xybad)[1]]

print a
  [[0 1 2 3 4 5 6 7 8 9]
  [10 11 12 13 14 15 16 17 18 19]
  [20 21 22 13 14 15 16 17 28 29]
  [30 31 32 32 44 45 46 38 38 39]
  [40 41 42 43 44 45 46 47 48 49]
  [50 51 52 53 54 55 56 57 58 59]
  [60 61 62 63 64 65 66 67 68 69]
  [70 71 72 73 74 75 76 77 78 79]
  [80 81 82 83 84 85 86 87 78 89]
  [90 91 92 93 94 95 96 97 98 99]]

Answer 2

我通常使用距离变换，这个方法是Juh_在这个问题中明智地建议的。

虽然这个方法并不直接适用于带掩码的数组，但我觉得把它转化过来应该不会太难，而且它的效率很高，我在处理大约100百万像素的图像时没有遇到任何问题。

这里复制相关的方法供参考：

import numpy as np
from scipy import ndimage as nd

def fill(data, invalid=None):
    """
    Replace the value of invalid 'data' cells (indicated by 'invalid') 
    by the value of the nearest valid data cell

    Input:
        data:    numpy array of any dimension
        invalid: a binary array of same shape as 'data'. True cells set where data
                 value should be replaced.
                 If None (default), use: invalid  = np.isnan(data)

    Output: 
        Return a filled array. 
    """
    #import numpy as np
    #import scipy.ndimage as nd

    if invalid is None: invalid = np.isnan(data)

    ind = nd.distance_transform_edt(invalid, return_distances=False, return_indices=True)
    return data[tuple(ind)]

Answer 3

你可以使用 np.roll 来制作 a 的移动副本，然后利用布尔逻辑在掩码上找出需要填充的位置：

import numpy as np
import numpy.ma as ma

a = np.arange(100).reshape(10,10)
fill_value=-99
a[2:4,3:8] = fill_value
a[8,8] = fill_value
a = ma.masked_array(a,a==fill_value)
print(a)

# [[0 1 2 3 4 5 6 7 8 9]
#  [10 11 12 13 14 15 16 17 18 19]
#  [20 21 22 -- -- -- -- -- 28 29]
#  [30 31 32 -- -- -- -- -- 38 39]
#  [40 41 42 43 44 45 46 47 48 49]
#  [50 51 52 53 54 55 56 57 58 59]
#  [60 61 62 63 64 65 66 67 68 69]
#  [70 71 72 73 74 75 76 77 78 79]
#  [80 81 82 83 84 85 86 87 -- 89]
#  [90 91 92 93 94 95 96 97 98 99]]

for shift in (-1,1):
    for axis in (0,1):        
        a_shifted=np.roll(a,shift=shift,axis=axis)
        idx=~a_shifted.mask * a.mask
        a[idx]=a_shifted[idx]

print(a)

# [[0 1 2 3 4 5 6 7 8 9]
#  [10 11 12 13 14 15 16 17 18 19]
#  [20 21 22 13 14 15 16 28 28 29]
#  [30 31 32 43 44 45 46 47 38 39]
#  [40 41 42 43 44 45 46 47 48 49]
#  [50 51 52 53 54 55 56 57 58 59]
#  [60 61 62 63 64 65 66 67 68 69]
#  [70 71 72 73 74 75 76 77 78 79]
#  [80 81 82 83 84 85 86 87 98 89]
#  [90 91 92 93 94 95 96 97 98 99]]

---

如果你想使用更多的最近邻居，你可以尝试这样做：

neighbors=((0,1),(0,-1),(1,0),(-1,0),(1,1),(-1,1),(1,-1),(-1,-1),
           (0,2),(0,-2),(2,0),(-2,0))

注意，neighbors 中元素的顺序很重要。你可能希望用最近的邻居来填补缺失的值，而不是随便找个邻居。虽然可能有更聪明的方法来生成邻居序列，但我现在想不到。

a_copy=a.copy()
for hor_shift,vert_shift in neighbors:
    if not np.any(a.mask): break
    a_shifted=np.roll(a_copy,shift=hor_shift,axis=1)
    a_shifted=np.roll(a_shifted,shift=vert_shift,axis=0)
    idx=~a_shifted.mask*a.mask
    a[idx]=a_shifted[idx]

要注意的是，np.roll 会把底部的边缘滚动到顶部，所以如果顶部有缺失值，可能会用最底部的值来填补。如果这造成了问题，我需要再想想怎么解决。一个明显但不太聪明的办法是使用 if 语句，给边缘提供一个不同的可接受邻居序列……