Intersect multiple 2D np arrays for determining spatially contiguous zones

Question

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

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.

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)

Answer 1

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)

The result is

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

https://ideone.com/9JuHjt