在两个numpy数组中用唯一的值组合标记区域？

Question

我有两个带标签的2dnumpy数组a和{}形状相同。我想用类似于两个数组的GIS geometric union来重新标记数组b，这样在数组a和b中具有唯一值组合的单元格被分配新的唯一id:

我不关心输出中区域的具体编号，只要值都是唯一的。我在下面附上了示例数组和期望的输出：我的实际数据集要大得多，两个数组都有整数标签，范围从“1”到“200000”。到目前为止，我已经尝试过将数组id连接起来以形成惟一的值组合，但理想情况下，我希望以1、2、3…等形式输出一组简单的新id

import numpy as np
import matplotlib.pyplot as plt

# Example labelled arrays a and b
input_a = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0],
                    [0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0],
                    [0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0],
                    [0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0],
                    [0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0],
                    [0, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 0],
                    [0, 0, 3, 3, 3, 3, 2, 2, 2, 2, 0, 0],
                    [0, 0, 3, 3, 3, 3, 2, 2, 2, 2, 0, 0],
                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

input_b = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [0, 0, 1, 1, 1, 3, 3, 3, 3, 3, 0, 0],
                    [0, 0, 1, 1, 1, 3, 3, 3, 3, 3, 0, 0],
                    [0, 0, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
                    [0, 0, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
                    [0, 0, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
                    [0, 0, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

# Plot inputs
plt.imshow(input_a, cmap="spectral", interpolation='nearest')
plt.imshow(input_b, cmap="spectral", interpolation='nearest')

# Desired output, union of a and b
output = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                   [0, 0, 1, 1, 1, 2, 3, 3, 3, 3, 0, 0],
                   [0, 0, 1, 1, 1, 2, 3, 3, 3, 3, 0, 0],
                   [0, 0, 1, 1, 1, 4, 7, 7, 7, 7, 0, 0],
                   [0, 0, 5, 5, 5, 6, 7, 7, 7, 7, 0, 0],
                   [0, 0, 5, 5, 5, 6, 7, 7, 7, 7, 0, 0],
                   [0, 0, 5, 5, 5, 6, 7, 7, 7, 7, 0, 0],
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

# Plot desired output
plt.imshow(output, cmap="spectral", interpolation='nearest')

Answer 1

这里有一种方法可以在概念上用集合并集来实现，但不是用地理信息系统的几何联合来实现，因为这是在我回答之后提到的。在

列出所有可能唯一的2元组值的列表，其中一个来自a，另一个来自b。将列表中的每个元组映射到其中的索引。使用该映射创建联合数组。在

例如，假设a和b是数组，每个数组都包含范围（4）中的值，为了简单起见，假设它们具有相同的形状。然后：

v = range(4)
from itertools import permutations
p = list(permutations(v,2))
m = {}
for i,x in enumerate(p):
    m[x] = i
union = np.empty_like(a)
for i,x in np.ndenumerate(a):
    union[i] = m[(x,b[i])]

为了演示，用

^{pr2}$

制作：

a = array([[3, 0, 3],
           [1, 3, 2],
           [0, 0, 3]])

b = array([[1, 3, 1],
           [0, 0, 1],
           [2, 3, 0]])

m = {(0, 1): 0,
     (0, 2): 1,
     (0, 3): 2,
     (1, 0): 3,
     (1, 2): 4,
     (1, 3): 5,
     (2, 0): 6,
     (2, 1): 7,
     (2, 3): 8,
     (3, 0): 9,
     (3, 1): 10,
     (3, 2): 11}

union = array([[10,  2, 10],
               [ 3,  9,  7],
               [ 1,  2,  9]])

在这种情况下，联合应该大于或等于其组合的属性反映在数值的增加，而不是元素数量的增加。在

Answer 2

如果我正确地理解了情况，那么您正在寻找来自a和{}的唯一对。因此，来自a的1和来自b的1将在输出中有一个唯一的标记；1来自a和b的3将在输出中具有另一个唯一的标记。再看看问题中的期望输出，这里似乎还有一个附加的条件情况，如果b为零，则输出也为零，而不考虑唯一对。在

下面的实现试图解决所有这些问题-

c = a*(b.max()+1) + b
c[b==0] = 0
_,idx = np.unique(c,return_inverse= True)
out = idx.reshape(b.shape)

样本运行-

^{pr2}$

样地-

# Plot inputs
plt.figure()                                                    
plt.imshow(a, cmap="spectral", interpolation='nearest')
plt.figure() 
plt.imshow(b, cmap="spectral", interpolation='nearest')

# Plot output
plt.figure()
plt.imshow(out, cmap="spectral", interpolation='nearest')