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)