将多个二维 np 阵列相交以确定区域 [英] Intersect multiple 2D np arrays for determining zones
问题描述
使用这个可重现的小示例，到目前为止，我无法从 3 个数组生成一个新的整数数组，其中包含所有三个输入数组的唯一分组.
数组与地形属性相关:
将 numpy 导入为 npasp = np.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4)) #方面slp = np.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4)) #坡elv = np.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #海拔这个想法是使用 GIS 例程将地理等值线分解为 3 个不同的属性:
- 1-8 为方面(1=朝北，2=朝东北等)
- 坡度为 9-12(9=缓坡...12=最陡坡)
- 13-16 为海拔(13=最低海拔...16=最高海拔)
下面的小图试图描绘我所追求的结果(左下角显示的数组).请注意，图中给出的答案"只是一种可能的答案.我不关心结果数组中整数的最终排列，只要最终数组在每个行/列索引处包含一个整数，用于标识唯一分组.
例如，[0,1] 和 [0,2] 处的数组索引具有相同的纵横比、坡度和高程，因此在结果数组中接收相同的整数标识符.
网格中的每个位置都与一个元组相关联，该元组由来自asp、slp 和 elv.例如，左上角有元组 (8,9,13).我们想将这个元组映射到一个唯一标识这个元组的数字.
一种方法是将 (8,9,13) 视为 3D 数组的索引np.arange(9*13*17).reshape(9,13,17).选择了这个特定的阵列容纳 asp、slp 和 elv 中的最大值:
在[107]中:asp.max()+1出[107]:9在 [108] 中:slp.max()+1出[108]:13在 [110] 中:elv.max()+1出[110]:17现在我们可以将元组 (8,9,13) 映射到数字 1934:
在[113]中:x = np.arange(9*13*17).reshape(9,13,17)在 [114] 中:x[8,9,13]出[114]:1934如果我们对网格中的每个位置都这样做，那么我们会为每个位置获得一个唯一的编号.我们可以就此结束，让这些独特的数字作为标签.
或者，我们可以生成更小的整数标签(从 0 开始并增加 1)通过使用 np.unique与:return_inverse=True
uniqs, labels = np.unique(vals, return_inverse=True)标签 = 标签.reshape(vals.shape)例如，
将 numpy 导入为 npasp = np.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4)) #方面slp = np.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4)) #坡elv = np.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #海拔x = np.arange(9*13*17).reshape(9,13,17)vals = x[asp, slp, elv]uniqs, 标签 = np.unique(vals, return_inverse=True)标签 = 标签.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 用于索引中的项目])x = np.arange(M.prod()).reshape(M)vals = x[索引]uniqs, 标签 = np.unique(vals, return_inverse=True)标签 = 标签.reshape(vals.shape) 产生与上述相同的结果，但即使 asp、slp、elv 是浮点数和/或大整数也能工作.
最后，我们可以避免np.arange的生成:
x = np.arange(M.prod()).reshape(M)vals = x[索引]通过将 vals 计算为索引和步幅的乘积:
M = np.r_[1, M[:-1]]步幅 = M.cumprod()指数 = np.stack(指数，轴=-1)vals = (指数 * 步幅).sum(axis=-1)所以把它们放在一起:
将 numpy 导入为 npasp = np.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4)) #方面slp = np.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4)) #坡elv = np.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #海拔def find_labels(*arrs):指数 = [np.unique(arr, return_inverse=True)[1] for arr in arrs]M = np.array([item.max()+1 用于索引中的项目])M = np.r_[1, M[:-1]]步幅 = M.cumprod()指数 = np.stack(指数，轴=-1)vals = (指数 * 步幅).sum(axis=-1)uniqs, 标签 = np.unique(vals, return_inverse=True)标签 = 标签.reshape(arrs[0].shape)退货标签打印(find_labels(asp，slp，elv))# [[ 3 7 7 0]# [ 6 10 12 4]# [ 8 9 11 4]# [ 2 5 5 1]]Using this small reproducible example, I've so far been unable to generate a new integer array from 3 arrays that contains unique groupings across all three input arrays.
The arrays are related to topographic properties:
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
The idea is that the geographic contours are broken into 3 different properties using GIS routines:
- 1-8 for aspect (1=north facing, 2=northeast facing, etc.)
- 9-12 for slope (9=gentle slope...12=steepest slope)
- 13-16 for elevation (13=lowest elevations...16=highest elevations)
The small graphic below attempts to depict the kind of result I'm after (array shown in lower left). Note, the "answer" given in the graphic is but one possible answer. I'm not concerned about the final arrangement of integers in the resulting array so long as the final array contains an integer at each row/column index that identifies unique groupings.
For example, the array indexes at [0,1] and [0,2] have the same aspect, slope, and elevation and therefore receive the same integer identifier in the resulting array.
Does numpy have a built in routine for this kind of thing?
Each location in the grid is associated with a tuple composed of one value from
asp, slp and elv. For example, the upper left corner has tuple (8,9,13).
We would like to map this tuple to a number which uniquely identifies this tuple.
One way to do that would be to think of (8,9,13) as the index into the 3D array
np.arange(9*13*17).reshape(9,13,17). This particular array was chosen
to accommodate the largest values in asp, slp and 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
Now we can map the tuple (8,9,13) to the number 1934:
In [113]: x = np.arange(9*13*17).reshape(9,13,17)
In [114]: x[8,9,13]
Out[114]: 1934
If we do this for each location in the grid, then we get a unique number for each location. We could end right here, letting these unique numbers serve as labels.
Or, we can generate smaller integer labels (starting at 0 and increasing by 1)
by using np.unique with
return_inverse=True:
uniqs, labels = np.unique(vals, return_inverse=True)
labels = labels.reshape(vals.shape)
So, for example,
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)
yields
array([[11, 0, 0, 1],
[ 9, 12, 2, 3],
[10, 8, 5, 3],
[ 7, 6, 6, 4]])
The above method works fine as long as the values in asp, slp and elv are small integers. If the integers were too large, the product of their maximums could overflow the maximum allowable value one can pass to np.arange. Moreover, generating such a large array would be inefficient.
If the values were floats, then they could not be interpreted as indices into the 3D array x.
So to address these problems, use np.unique to convert the values in asp, slp and elv to unique integer labels first:
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)
which yields the same result as shown above, but works even if asp, slp, elv were floats and/or large integers.
Finally, we can avoid the generation of np.arange:
x = np.arange(M.prod()).reshape(M)
vals = x[indices]
by computing vals as a product of indices and strides:
M = np.r_[1, M[:-1]]
strides = M.cumprod()
indices = np.stack(indices, axis=-1)
vals = (indices * strides).sum(axis=-1)
So putting it all together:
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]]
这篇关于将多个二维 np 阵列相交以确定区域的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
