相交多个2D np阵列以确定区域

Question

使用这个可重复的小例子,我到目前为止还无法从3个数组生成一个新的整数数组,这些数组包含所有三个输入数组的唯一分组.

这些数组与地形属性有关:

import numpy as np
asp = np.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4))  #aspect
slp = np.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4))  #slope
elv = np.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #elevation

我们的想法是使用GIS例程将地理轮廓分解为3个不同的属性:

• 方位1-8(1 =朝北,2 =东北朝向等)
• 坡度为9-12(9 =缓坡... 12 =最陡斜坡)
• 海拔13-16(13 =最低海拔... 16 =最高海拔)

下面的小图试图描绘我所追求的结果类型(数组显示在左下方).请注意,图中给出的"答案"只是一个可能的答案.只要最终数组在每个行/列索引处包含一个标识唯一分组的整数,我就不关心结果数组中整数的最终排列.

例如,[0,1]和[0,2]处的数组索引具有相同的方面,斜率和高程,因此在结果数组中接收相同的整数标识符.

是否numpy的有一个内置的例程这种事情？

Answer 1

网格中的每个位置都与一个由 asp、slp和中的一个值组成的元组相关联elv。例如，左上角有 tuple (8,9,13)。我们希望将此元组映射到唯一标识该元组的数字。

一种方法是将其视为(8,9,13)3D 数组的索引 np.arange(9*13*17).reshape(9,13,17)。asp选择这个特定数组是为了容纳、slp和中的最大值elv：

In [107]: asp.max()+1
Out[107]: 9

In [108]: slp.max()+1
Out[108]: 13

In [110]: elv.max()+1
Out[110]: 17

现在我们可以将元组 (8,9,13) 映射到数字 1934：

In [113]: x = np.arange(9*13*17).reshape(9,13,17)

In [114]: x[8,9,13]
Out[114]: 1934

如果我们对网格中的每个位置执行此操作，那么我们会为每个位置获得一个唯一的编号。我们可以在这里结束，让这些唯一的数字作为标签。

np.unique或者，我们可以使用with return_inverse=True生成更小的整数标签（从 0 开始并增加 1）：

uniqs, labels = np.unique(vals, return_inverse=True)
labels = labels.reshape(vals.shape)

所以，举例来说，

import numpy as np

asp = np.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4))  #aspect
slp = np.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4))  #slope
elv = np.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #elevation

x = np.arange(9*13*17).reshape(9,13,17)
vals = x[asp, slp, elv]
uniqs, labels = np.unique(vals, return_inverse=True)
labels = labels.reshape(vals.shape)

产量

array([[11,  0,  0,  1],
       [ 9, 12,  2,  3],
       [10,  8,  5,  3],
       [ 7,  6,  6,  4]])

只要 、和中的值是小整数asp，上述方法就可以正常工作。如果整数太大，它们最大值的乘积可能会溢出可以传递给 的最大允许值。而且，生成如此大的数组效率很低。如果值是浮点数，则它们不能被解释为 3D 数组的索引。slpelvnp.arangex

因此，要解决这些问题，请首先使用将,和np.unique中的值转换为唯一的整数标签：aspslpelv

indices = [ np.unique(arr, return_inverse=True)[1].reshape(arr.shape) for arr in [asp, slp, elv] ]
M = np.array([item.max()+1 for item in indices])
x = np.arange(M.prod()).reshape(M)
vals = x[indices]
uniqs, labels = np.unique(vals, return_inverse=True)
labels = labels.reshape(vals.shape)

它产生与上面所示相同的结果，但即使asp, slp,elv是浮点数和/或大整数也能工作。

最后，我们可以避免生成np.arange：

x = np.arange(M.prod()).reshape(M)
vals = x[indices]

通过计算vals索引和步幅的乘积：

M = np.r_[1, M[:-1]]
strides = M.cumprod()
indices = np.stack(indices, axis=-1)
vals = (indices * strides).sum(axis=-1)

所以把它们放在一起：

import numpy as np

asp = np.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4))  #aspect
slp = np.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4))  #slope
elv = np.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #elevation

def find_labels(*arrs):
    indices = [np.unique(arr, return_inverse=True)[1] for arr in arrs]
    M = np.array([item.max()+1 for item in indices])
    M = np.r_[1, M[:-1]]
    strides = M.cumprod()
    indices = np.stack(indices, axis=-1)
    vals = (indices * strides).sum(axis=-1)
    uniqs, labels = np.unique(vals, return_inverse=True)
    labels = labels.reshape(arrs[0].shape)
    return labels

print(find_labels(asp, slp, elv))

# [[ 3  7  7  0]
#  [ 6 10 12  4]
#  [ 8  9 11  4]
#  [ 2  5  5  1]]