相交多个 2D np 阵列以确定空间上连续的区域

Question

这个问题的背景与a past post有关。选定的答案提供了一个很好的解决方案。但是，我还想补充一点细微差别。

这些数组没有使用上面linked-to post 中给出的数组，而是突出显示上一组解决方案无法按照我想要的方式处理的情况。新的数组是：

import numpy as np

asp = np.array([2,1,8,2,5,4,6,4,3,3,3,3,2,1,8,2]).reshape((4,4))  #aspect
slp = np.array([11,11,12,11,12,11,12,10,9,9,9,10,11,10,11,11]).reshape((4,4))  #slope
elv = np.array([16,15,15,16,17,15,15,15,15,15,14,14,16,15,16,16]).reshape((4,4)) #elevation

以前的解决方案在高程、坡度和坡向相同的地方分配了唯一的区域编号，无论这些唯一区域在空间上是否连续。下图试图显示此结果（请参见左下角 2 个数组中的上一个）。所需的结果将类似于图像左下角的下部数组。例如，请注意每个角单元具有相同的高程、坡度和坡向，但在空间上并不连续，因此应为其分配唯一的区域 ID。 “空间连续”被定义为具有相同高程、坡度和坡向的 2 个（或更多）相邻（相邻或对角线）像元。下图中以红色突出显示的单元格描述了这一点。

在之前没有考虑结果区域中空间连续性的解决方案中，对我来说最有意义的是@Stephen Rauch。他的解决方案是（除其他外）：

combined = 10000 * asp + 100 * slp + elv
unique = dict(((v, i + 1) for i, v in enumerate(np.unique(combined))))
combined_unique = np.vectorize(unique.get)(combined)

Answer 1

您可以在计算区域后对区域进行填充。我最初的想法是在唯一性计算之后执行此操作，但鉴于您只对连续性感兴趣，请使用原始数据并完全跳过第一个分区步骤。以下实现不是超级有效。我确信可以做一些事情来优化它：

from collections import deque
import numpy as np

asp = np.array([[ 2, 1, 8, 2],
                [ 5, 4, 6, 4],
                [ 3, 3, 3, 3],
                [ 2, 1, 8, 2]])  #aspect
slp = np.array([[11,11,12,11],
                [12,11,12,10],
                [ 9, 9, 9,10],
                [11,10,11,11]])  #slope
elv = np.array([[16,15,15,16],
                [17,15,15,15],
                [15,15,14,14],
                [16,15,16,16]]) #elevation

criteria = np.stack((asp, slp, elv))

def floodfill(n, start):
    def check(r, c):
        return result[r, c] == 0 and np.array_equal(crit, criteria[:, r, c])

    stack = deque([start])
    crit = criteria [(slice(None), *start)]
    while stack:
        r, c = stack.popleft()
        if result[r, c]: continue
        result[r, c] = n
        if r > 0 and check(r - 1, c):
            stack.append((r - 1, c))
        if r < result.shape[0] - 1 and check(r + 1, c):
            stack.append((r + 1, c))
        if c > 0 and check(r, c - 1):
            stack.append((r, c - 1))
        if c < result.shape[1] - 1 and check(r, c + 1):
            stack.append((r, c + 1))


result = np.zeros_like(asp)
count = 1
it = np.nditer(result, flags=['multi_index'])
for x in it:
    if x == 0:
        floodfill(count, it.multi_index)
        count += 1

print(result)

结果是

[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9  9 10 11]
 [12 13 14 15]]

https://ideone.com/9JuHjt