与多个2D np阵列相交以确定空间连续区域 [英] Intersect multiple 2D np arrays for determining spatially contiguous zones
问题描述
此问题的背景与过去的帖子.选择的答案提供了一个很好的解决方案.但是，我想补充一点细微差别.
The background to this question is related to a past post. The selected answer provides a nice solution. However, there is a further nuance I would like to add.
代替使用链接到帖子上方，这些数组突出显示了一个条件，即先前的解决方案集无法完全按照我的意愿进行处理.新的数组是:
Instead of using the arrays given in the linked-to post above, these arrays highlight a condition that the previous set of solutions don't handle quite the way I would like. The new arrays are:
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个(或更多)具有相同高程，坡度和长宽比的相邻(相邻或对角线)单元.下图中以红色突出显示的单元格对此进行了说明.
The previous solutions assigned unique zone numbers wherever elevation, slope, and aspect were the same, regardless of whether the unique zones were spatially contiguous. The image below attempts to show this result (see the upper of 2 arrays in lower left corner). The desired result would be something like the lower array of the lower left corner of the image. For example, note that each of the corner cells have the same elevation, slope, and aspect, but are not spatially contiguous and should therefore be assigned unique zone IDs of their own. "Spatially contiguous" is defined as 2 (or more) neighboring (adjacent or diagonal) cells that have the same elevation, slope, and aspect. The cells highlighted in red in the image below depict this.
在先前的解决方案中，没有考虑到所得区域中的空间连续性，其中最有意义的解决方案是@Stephen Rauch.他的解决方案是(其他):
Among the previous solutions that don't account for spatial continuity in the resulting zones, the one that made the most sense to me was by @Stephen Rauch. His solution was (among others):
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)
推荐答案
一旦计算出区域，便可以填充区域.我最初的想法是在唯一性计算之后执行此操作，但是鉴于您仅对连续性感兴趣，请使用原始数据并完全跳过第一个分区步骤.以下实现并非超级有效.我确定可以做一些优化工作:
You could floodfill the regions once they are computed. My original idea would have been to do this after the uniqueness computation, but given that you are only interested in contiguity, use the original data and skip the first zoning step entirely. The following implementation is not super efficient. I'm sure something could be done to optimize it:
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]]
这篇关于与多个2D np阵列相交以确定空间连续区域的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
