scipy python中多维数组的最近邻中值和计数

byqmnocz  于 2022-11-10  发布在  Python
关注(0)|答案(1)|浏览(157)

从一个用浮点数填充的数组(在本例中为简单起见使用int)中,目标是创建两个新的np.数组(与data具有相同的形状):

  • 第一个np.数组,其中值为元素周围的最近邻居和元素本身,因此为9个值。
  • 第二个np.数组,包含相邻元素的计数和元素本身。没有值(当你超出边界时,没有值)或-99不应被计算在内。

我现在正在使用的示例代码:(注意,“数据”的真实的形状是360,802,填充有浮点数)

data= np.array([
    [  5,  8,  5,  3, -1,  7, -5,-99, -2,  5], 
    [ -2,  4,  1,-99,  1,  9, -1, -2,  5,  9],
    [  0,  4,  1,  0, -3,  2, -1,  1,  3, -1],
    [ -2,  9, -1,  3,  1,  3,  2,  0, -7,  2],
    [  4,-99,  4,  0, -9, -1, -1,-99,  1, -4],
    [  1, -4,  7,  1,  6,  6, -2, -9,  4,  9]
    ], np.float32)

distFmOriginX=1
distFmOriginY=1

medianArr       = np.empty([data.shape[0],data.shape[1]], dtype = float)
nrOfCountsArr   = np.zeros([data.shape[0],data.shape[1]], dtype = int)

for x in range(data.shape[0]):
 for y in range(data.shape[1]):
  cnt=0
  arr=[]
  if data[x][y]==-99:
   medianArr[x][y]=-99
   nrOfCountsArr[x][y]=-99
   continue
  xAlt = x-distFmOriginX 
  yAlt = y-distFmOriginY

  for i in range((2*distFmOriginX)+1):
   for j in range((2*distFmOriginY)+1):
    if (xAlt+i)>=data.shape[0] or (xAlt+i)<0:
     continue
    if (yAlt+j)>=data.shape[1] or (yAlt+j)<0:
     continue
    if data[xAlt+i][yAlt+j]==-99:
     continue
    arr.append(data[xAlt+i][yAlt+j])
    cnt+=1
  if cnt>0:
   medianArr[x][y]     = np.median(arr)
   nrOfCountsArr[x][y] = cnt

这导致了期望的结果:

medianArr:
 [[  4.5   4.5   4.    1.    3.    0.   -1.  -99.    5.    5. ]
 [  4.    4.    3.5 -99.    1.5  -1.    0.   -1.    2.    4. ]
 [  2.    1.    2.    1.    1.5   1.    1.    0.    1.    2.5]
 [  4.    2.5   2.    0.    0.   -1.    0.5   0.5   0.5   0. ]
 [  1.  -99.    2.    1.    1.    1.   -0.5 -99.    0.5   1.5]
 [  1.    4.    1.    2.5   0.5  -1.   -1.   -1.    1.    2.5]]
nrOfCountsArr:
 [[  4   6   5   5   5   6   5 -99   5   4]
 [  6   9   8 -99   8   9   8   8   8   6]
 [  6   9   8   8   8   9   9   9   9   6]
 [  5   8   8   9   9   9   8   8   8   6]
 [  5 -99   8   9   9   9   8 -99   8   6]
 [  3   5   5   6   6   6   5   5   5   4]]

现在,上面的工作解决方案与老式的多嵌套for循环是适当的。这需要太多的处理时间,我敢肯定有一个麻木/scipy的解决方案,我不知道。搜索了很多交付这个解决方案,但这是平均值,而不是中值:
https://gis.stackexchange.com/questions/254753/calculate-the-average-of-neighbor-pixels-for-raster-edge

import numpy as np
from scipy.signal import convolve2d

def average_convolve2d(data):
    kernel = np.ones((3, 3))

    array_sum = convolve2d(
        data, kernel, mode='same',
        boundary='fill', fillvalue=0)

    num_counts = convolve2d(
        np.ones(data.shape), kernel, mode='same',
        boundary='fill', fillvalue=0)

    return array_sum / num_counts

avg = average_convolve2d(data)

这个方法非常好,而且非常快(0.01s对10.5s)。我一直在想,我是否可以对中值(而不是本例中的平均值)做同样或类似的操作,并返回两个数组(中值和窗口内的数字计数)。

  • 注意:将来我希望能够指定在x和y维度上远离主元素的元素数量。例如:在x维度上向左和向右移动1步(distFmOriginX=1),在y维度上向上和向下移动2步(distFmOriginY=2),包括主元素,得到15个值。*

希望你们中有人能想出一个我还没想到的绝妙的解决方案。

14ifxucb

14ifxucb1#

请使用以下选项:

import numpy as np

data = np.array([
    [5, 8, 5, 3, -1, 7, -5, -99, -2, 5],
    [-2, 4, 1, -99, 1, 9, -1, -2, 5, 9],
    [0, 4, 1, 0, -3, 2, -1, 1, 3, -1],
    [-2, 9, -1, 3, 1, 3, 2, 0, -7, 2],
    [4, -99, 4, 0, -9, -1, -1, -99, 1, -4],
    [1, -4, 7, 1, 6, 6, -2, -9, 4, 9]
], np.float32)

# create padded array with np.nan around the edges

padded = np.pad(data, 1, 'constant', constant_values=np.nan)

# set -99 to np.nan to ignore in the median

padded[padded == -99] = np.nan

# create sliding window over the padded array

window = np.lib.stride_tricks.sliding_window_view(padded, (3, 3)).reshape((6, 10, -1))

# compute the median

median = np.nanmedian(window, axis=2)

# compute the count

count = np.count_nonzero(~np.isnan(window), axis=2)

# reset the elements to -99

median[data == -99] = -99
count[data == -99] = -99

print(median)
print("--")
print(count)

输出

[[  4.5   4.5   4.    1.    3.    0.   -1.  -99.    5.    5. ]
 [  4.    4.    3.5 -99.    1.5  -1.    0.   -1.    2.    4. ]
 [  2.    1.    2.    1.    1.5   1.    1.    0.    1.    2.5]
 [  4.    2.5   2.    0.    0.   -1.    0.5   0.5   0.5   0. ]
 [  1.  -99.    2.    1.    1.    1.   -0.5 -99.    0.5   1.5]
 [  1.    4.    1.    2.5   0.5  -1.   -1.   -1.    1.    2.5]]
--
[[  4   6   5   5   5   6   5 -99   5   4]
 [  6   9   8 -99   8   9   8   8   8   6]
 [  6   9   8   8   8   9   9   9   9   6]
 [  5   8   8   9   9   9   8   8   8   6]
 [  5 -99   8   9   9   9   8 -99   8   6]
 [  3   5   5   6   6   6   5   5   5   4]]

循序渐进

1.我们的想法是使用pad + sliding_window_view创建窗口数组,并将要忽略的值设置为nan
1.第二步是使用约简函数nanmediancount_nonzero来计算所需的结果。
1.最后设置回-99值以匹配精确的输出。

相关问题