Пересечение нескольких 2D массивов NP для определения пространственно смежных зон

Question

История вопроса связана с прошлым сообщением .Выбранный ответ обеспечивает хорошее решение.Однако есть еще один нюанс, который я хотел бы добавить.

Вместо использования массивов, указанных в сообщении , связанного с выше, эти массивы подчеркивают условие, что предыдущий набор решений обрабатывает не совсем так, как мне хотелось бы.Новые массивы:

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 массивов в левом нижнем углу).Желаемый результат будет что-то вроде нижнего массива левого нижнего угла изображения.Например, обратите внимание, что каждая из угловых ячеек имеет одинаковую высоту, уклон и аспект, но не является пространственно смежной, и поэтому ей должны быть назначены собственные уникальные идентификаторы зоны.«Пространственно смежный» определяется как 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