交叉多个 2D np 数组以确定区域
Intersect multiple 2D np arrays for determining zones
使用这个可重现的小示例，到目前为止，我无法从 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=最高海拔）
下面的小图试图描述我想要的那种结果（数组显示在左下角）。请注意，图中给出的 "answer" 只是一种可能的答案。我不关心结果数组中整数的最终排列，只要最终数组在每个 row/column 索引处包含一个标识唯一分组的整数。
例如，[0,1] 和 [0,2] 处的数组索引具有相同的坡向、坡度和高程，因此在结果数组中接收相同的整数标识符。
numpy是否有针对此类事情的内置例程？
这似乎与标记图像中的独特区域类似。这是我为此编写的函数，但您首先需要将 3 个数组连接到 1 个 3D 数组。
def labelPix(pix):
height, width, _ = pix.shape
pixRows = numpy.reshape(pix, (height * width, 3))
unique, counts = numpy.unique(pixRows, return_counts = True, axis = 0)
unique = [list(elem) for elem in unique]
labeledPix = numpy.zeros((height, width), dtype = int)
offset = 0
for index, zoneArray in enumerate(unique):
index += offset
zone = list(zoneArray)
zoneArea = (pix == zone).all(-1)
elementsArray, numElements = scipy.ndimage.label(zoneArea)
elementsArray[elementsArray!=0] += offset
labeledPix[elementsArray!=0] = elementsArray[elementsArray!=0]
offset += numElements
return labeledPix
这将标记唯一的 3 值组合，同时还将单独的标签分配给具有相同 3 值组合但彼此不联系的区域。
asp = numpy.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4)) #aspect
slp = numpy.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4)) #slope
elv = numpy.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #elevation
pix = numpy.zeros((4,4,3))
pix[:,:,0] = asp
pix[:,:,1] = slp
pix[:,:,2] = elv
print(labelPix(pix))
returns:
[[ 0 1 1 2]
[10 12 3 4]
[11 9 6 4]
[ 8 7 7 5]]
这可以使用 numpy.unique() 然后像这样的映射来完成：
代码：
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)
测试代码：
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))
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)
print(combined_unique)
结果：
[[12 1 1 2]
[10 13 3 4]
[11 9 6 4]
[ 8 7 7 5]]
网格中的每个位置都与一个元组相关联，该元组由一个值组成
asp、slp 和 elv。比如左上角有元组(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
如果我们对网格中的每个位置都这样做，那么我们会为每个位置获得一个唯一的编号。
我们可以在这里结束，让这些唯一的数字作为标签。
或者，我们可以生成更小的整数标签（从 0 开始，递增 1）
通过使用 np.unique
return_inverse=True:
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、slp和elv中的值是小整数，上述方法就可以正常工作。如果整数太大，它们的最大值的乘积可能会溢出可以传递给 np.arange 的最大允许值。此外，生成如此大的数组效率很低。
如果值是浮点数，则它们不能被解释为 3D 数组的索引 x。
所以为了解决这些问题，首先使用np.unique将asp、slp和elv中的值转换为唯一的整数标签：
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 是浮点数 and/or 大整数也有效。
终于可以避免生成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]]
这是使用 itertools.groupby 的简单 Python 技术。它要求输入是一维列表，但这不应该是一个主要问题。策略是将列表压缩在一起，连同索引号，然后对结果列进行排序。然后我们将相同的列组合在一起，在比较列时忽略索引号。然后我们收集每个组的索引号，并使用它们构建最终输出列表。
from itertools import groupby
def show(label, seq):
print(label, ' '.join(['{:2}'.format(u) for u in seq]))
asp = [8, 1, 1, 2, 7, 8, 2, 3, 7, 6, 4, 3, 6, 5, 5, 4]
slp = [9, 10, 10, 9, 9, 12, 12, 9, 10, 11, 11, 9, 9, 9, 9, 9]
elv = [13, 14, 14, 13, 14, 15, 16, 14, 14, 15, 16, 14, 13, 14, 14, 13]
size = len(asp)
a = sorted(zip(asp, slp, elv, range(size)))
groups = sorted([u[-1] for u in g] for _, g in groupby(a, key=lambda t:t[:-1]))
final = [0] * size
for i, g in enumerate(groups, 1):
for j in g:
final[j] = i
show('asp', asp)
show('slp', slp)
show('elv', elv)
show('out', final)
输出
asp 8 1 1 2 7 8 2 3 7 6 4 3 6 5 5 4
slp 9 10 10 9 9 12 12 9 10 11 11 9 9 9 9 9
elv 13 14 14 13 14 15 16 14 14 15 16 14 13 14 14 13
out 1 2 2 3 4 5 6 7 8 9 10 7 11 12 12 13
没有必要做第二个排序，我们可以只使用一个普通的列表 comp
groups = [[u[-1] for u in g] for _, g in groupby(a, key=lambda t:t[:-1])]
或生成器表达式
groups = ([u[-1] for u in g] for _, g in groupby(a, key=lambda t:t[:-1]))
我这样做只是为了让我的输出与问题中的输出匹配。
这是使用基于字典的查找解决此问题的一种方法。
from collections import defaultdict
import itertools
group_dict = defaultdict(list)
idx_count = 0
for a, s, e in np.nditer((asp, slp, elv)):
asp_tuple = (a.tolist(), s.tolist(), e.tolist())
if asp_tuple not in group_dict:
group_dict[asp_tuple] = [idx_count+1]
idx_count += 1
else:
group_dict[asp_tuple].append(group_dict[asp_tuple][-1])
list1d = list(itertools.chain(*list(group_dict.values())))
np.array(list1d).reshape(4, 4)
# result
array([[ 1, 2, 2, 3],
[ 4, 5, 6, 7],
[ 7, 8, 9, 10],
[11, 12, 12, 13]])
使用这个可重现的小示例，到目前为止，我无法从 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=最高海拔）
下面的小图试图描述我想要的那种结果（数组显示在左下角）。请注意，图中给出的 "answer" 只是一种可能的答案。我不关心结果数组中整数的最终排列，只要最终数组在每个 row/column 索引处包含一个标识唯一分组的整数。
例如，[0,1] 和 [0,2] 处的数组索引具有相同的坡向、坡度和高程，因此在结果数组中接收相同的整数标识符。
numpy是否有针对此类事情的内置例程？
这似乎与标记图像中的独特区域类似。这是我为此编写的函数，但您首先需要将 3 个数组连接到 1 个 3D 数组。
def labelPix(pix):
height, width, _ = pix.shape
pixRows = numpy.reshape(pix, (height * width, 3))
unique, counts = numpy.unique(pixRows, return_counts = True, axis = 0)
unique = [list(elem) for elem in unique]
labeledPix = numpy.zeros((height, width), dtype = int)
offset = 0
for index, zoneArray in enumerate(unique):
index += offset
zone = list(zoneArray)
zoneArea = (pix == zone).all(-1)
elementsArray, numElements = scipy.ndimage.label(zoneArea)
elementsArray[elementsArray!=0] += offset
labeledPix[elementsArray!=0] = elementsArray[elementsArray!=0]
offset += numElements
return labeledPix
这将标记唯一的 3 值组合，同时还将单独的标签分配给具有相同 3 值组合但彼此不联系的区域。
asp = numpy.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4)) #aspect
slp = numpy.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4)) #slope
elv = numpy.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #elevation
pix = numpy.zeros((4,4,3))
pix[:,:,0] = asp
pix[:,:,1] = slp
pix[:,:,2] = elv
print(labelPix(pix))
returns:
[[ 0 1 1 2]
[10 12 3 4]
[11 9 6 4]
[ 8 7 7 5]]
这可以使用 numpy.unique() 然后像这样的映射来完成：
代码：
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)
测试代码：
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))
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)
print(combined_unique)
结果：
[[12 1 1 2]
[10 13 3 4]
[11 9 6 4]
[ 8 7 7 5]]
网格中的每个位置都与一个元组相关联，该元组由一个值组成
asp、slp 和 elv。比如左上角有元组(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
如果我们对网格中的每个位置都这样做，那么我们会为每个位置获得一个唯一的编号。 我们可以在这里结束，让这些唯一的数字作为标签。
或者，我们可以生成更小的整数标签（从 0 开始，递增 1）
通过使用 np.unique
return_inverse=True:
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、slp和elv中的值是小整数，上述方法就可以正常工作。如果整数太大，它们的最大值的乘积可能会溢出可以传递给 np.arange 的最大允许值。此外，生成如此大的数组效率很低。
如果值是浮点数，则它们不能被解释为 3D 数组的索引 x。
所以为了解决这些问题，首先使用np.unique将asp、slp和elv中的值转换为唯一的整数标签：
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 是浮点数 and/or 大整数也有效。
终于可以避免生成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]]
这是使用 itertools.groupby 的简单 Python 技术。它要求输入是一维列表，但这不应该是一个主要问题。策略是将列表压缩在一起，连同索引号，然后对结果列进行排序。然后我们将相同的列组合在一起，在比较列时忽略索引号。然后我们收集每个组的索引号，并使用它们构建最终输出列表。
from itertools import groupby
def show(label, seq):
print(label, ' '.join(['{:2}'.format(u) for u in seq]))
asp = [8, 1, 1, 2, 7, 8, 2, 3, 7, 6, 4, 3, 6, 5, 5, 4]
slp = [9, 10, 10, 9, 9, 12, 12, 9, 10, 11, 11, 9, 9, 9, 9, 9]
elv = [13, 14, 14, 13, 14, 15, 16, 14, 14, 15, 16, 14, 13, 14, 14, 13]
size = len(asp)
a = sorted(zip(asp, slp, elv, range(size)))
groups = sorted([u[-1] for u in g] for _, g in groupby(a, key=lambda t:t[:-1]))
final = [0] * size
for i, g in enumerate(groups, 1):
for j in g:
final[j] = i
show('asp', asp)
show('slp', slp)
show('elv', elv)
show('out', final)
输出
asp 8 1 1 2 7 8 2 3 7 6 4 3 6 5 5 4
slp 9 10 10 9 9 12 12 9 10 11 11 9 9 9 9 9
elv 13 14 14 13 14 15 16 14 14 15 16 14 13 14 14 13
out 1 2 2 3 4 5 6 7 8 9 10 7 11 12 12 13
没有必要做第二个排序，我们可以只使用一个普通的列表 comp
groups = [[u[-1] for u in g] for _, g in groupby(a, key=lambda t:t[:-1])]
或生成器表达式
groups = ([u[-1] for u in g] for _, g in groupby(a, key=lambda t:t[:-1]))
我这样做只是为了让我的输出与问题中的输出匹配。
这是使用基于字典的查找解决此问题的一种方法。
from collections import defaultdict
import itertools
group_dict = defaultdict(list)
idx_count = 0
for a, s, e in np.nditer((asp, slp, elv)):
asp_tuple = (a.tolist(), s.tolist(), e.tolist())
if asp_tuple not in group_dict:
group_dict[asp_tuple] = [idx_count+1]
idx_count += 1
else:
group_dict[asp_tuple].append(group_dict[asp_tuple][-1])
list1d = list(itertools.chain(*list(group_dict.values())))
np.array(list1d).reshape(4, 4)
# result
array([[ 1, 2, 2, 3],
[ 4, 5, 6, 7],
[ 7, 8, 9, 10],
[11, 12, 12, 13]])