我不认为用 0 替换是这样做的正确方法,您正在将 covolved 值推向 0。这些缺失应该被视为“缺失”。因为它们代表了缺失的信息片段,没有理由仅仅假设它们可能是 0,它们根本不应该参与任何计算。
我尝试将缺失值设置为numpy.nan,然后进行卷积,结果表明内核之间的任何重叠和任何缺失都会在结果中给出nan,即使重叠是内核中的0,所以你会在结果中得到一个扩大的缺失洞。根据您的应用,这可能是您想要的结果。
但在某些情况下,您不希望仅仅因为 1 次缺失就丢弃这么多信息(也许 astropy 具有更好的实现:numpy.nans 被忽略(或替换为插值?)。
因此,使用astropy,您将执行以下操作:
from astropy.convolution import convolve
inarray=numpy.where(inarray.mask,numpy.nan,inarray) # masking still doesn't work, has to set to numpy.nan
result=convolve(inarray,kernel)
但是,您仍然无法控制可以容忍的缺失程度。为了实现这一点,我创建了一个函数,它使用scipy.ndimage.convolve() 进行初始卷积,但只要涉及缺失(numpy.nan),就需要手动重新计算值:
def convolve2d(slab,kernel,max_missing=0.5,verbose=True):
'''2D convolution with missings ignored
<slab>: 2d array. Input array to convolve. Can have numpy.nan or masked values.
<kernel>: 2d array, convolution kernel, must have sizes as odd numbers.
<max_missing>: float in (0,1), max percentage of missing in each convolution
window is tolerated before a missing is placed in the result.
Return <result>: 2d array, convolution result. Missings are represented as
numpy.nans if they are in <slab>, or masked if they are masked
in <slab>.
'''
from scipy.ndimage import convolve as sciconvolve
assert numpy.ndim(slab)==2, "<slab> needs to be 2D."
assert numpy.ndim(kernel)==2, "<kernel> needs to be 2D."
assert kernel.shape[0]%2==1 and kernel.shape[1]%2==1, "<kernel> shape needs to be an odd number."
assert max_missing > 0 and max_missing < 1, "<max_missing> needs to be a float in (0,1)."
#--------------Get mask for missings--------------
if not hasattr(slab,'mask') and numpy.any(numpy.isnan(slab))==False:
has_missing=False
slab2=slab.copy()
elif not hasattr(slab,'mask') and numpy.any(numpy.isnan(slab)):
has_missing=True
slabmask=numpy.where(numpy.isnan(slab),1,0)
slab2=slab.copy()
missing_as='nan'
elif (slab.mask.size==1 and slab.mask==False) or numpy.any(slab.mask)==False:
has_missing=False
slab2=slab.copy()
elif not (slab.mask.size==1 and slab.mask==False) and numpy.any(slab.mask):
has_missing=True
slabmask=numpy.where(slab.mask,1,0)
slab2=numpy.where(slabmask==1,numpy.nan,slab)
missing_as='mask'
else:
has_missing=False
slab2=slab.copy()
#--------------------No missing--------------------
if not has_missing:
result=sciconvolve(slab2,kernel,mode='constant',cval=0.)
else:
H,W=slab.shape
hh=int((kernel.shape[0]-1)/2) # half height
hw=int((kernel.shape[1]-1)/2) # half width
min_valid=(1-max_missing)*kernel.shape[0]*kernel.shape[1]
# dont forget to flip the kernel
kernel_flip=kernel[::-1,::-1]
result=sciconvolve(slab2,kernel,mode='constant',cval=0.)
slab2=numpy.where(slabmask==1,0,slab2)
#------------------Get nan holes------------------
miss_idx=zip(*numpy.where(slabmask==1))
if missing_as=='mask':
mask=numpy.zeros([H,W])
for yii,xii in miss_idx:
#-------Recompute at each new nan in result-------
hole_ys=range(max(0,yii-hh),min(H,yii+hh+1))
hole_xs=range(max(0,xii-hw),min(W,xii+hw+1))
for hi in hole_ys:
for hj in hole_xs:
hi1=max(0,hi-hh)
hi2=min(H,hi+hh+1)
hj1=max(0,hj-hw)
hj2=min(W,hj+hw+1)
slab_window=slab2[hi1:hi2,hj1:hj2]
mask_window=slabmask[hi1:hi2,hj1:hj2]
kernel_ij=kernel_flip[max(0,hh-hi):min(hh*2+1,hh+H-hi),
max(0,hw-hj):min(hw*2+1,hw+W-hj)]
kernel_ij=numpy.where(mask_window==1,0,kernel_ij)
#----Fill with missing if not enough valid data----
ksum=numpy.sum(kernel_ij)
if ksum<min_valid:
if missing_as=='nan':
result[hi,hj]=numpy.nan
elif missing_as=='mask':
result[hi,hj]=0.
mask[hi,hj]=True
else:
result[hi,hj]=numpy.sum(slab_window*kernel_ij)
if missing_as=='mask':
result=numpy.ma.array(result)
result.mask=mask
return result
下图展示了输出。左边是一个 30x30 的随机地图,有 3 个numpy.nan 洞,大小为:
- 1x1
- 3x3
- 5x5
右侧是卷积输出,由 5x5 内核(全为 1)和 50% 的容差水平 (max_missing=0.5)。
所以前 2 个较小的孔使用附近的值填充,而在最后一个中,因为缺失的数量 > 0.5x5x5 = 12.5、numpy.nans 被放置来表示缺失的信息。