1

我正在尝试查看是否numpy.histogram2d会为我将 2 个数组中的数据交叉制表。我以前从未使用过此功能，并且遇到错误，我不知道如何解决。

import numpy as np
import random
zones = np.zeros((20,30), int)
values = np.zeros((20,30), int)
for i in range(20):
    for j in range(30):
        values[i,j] = random.randint(0,10)
zones[:8,:15] = 100
zones[8:,:15] = 101
zones[:8,15:] = 102
zones[8:,15:] = 103
np.histogram2d(zones,values)

此代码导致以下错误：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-18-53447df32000> in <module>()
----> 1 np.histogram2d(zones,values)

C:\Python27\ArcGISx6410.2\lib\site-packages\numpy\lib\twodim_base.pyc in histogram2d(x, y, bins, range, normed, weights)
    613         xedges = yedges = asarray(bins, float)
    614         bins = [xedges, yedges]
--> 615     hist, edges = histogramdd([x,y], bins, range, normed, weights)
    616     return hist, edges[0], edges[1]
    617 

C:\Python27\ArcGISx6410.2\lib\site-packages\numpy\lib\function_base.pyc in histogramdd(sample, bins, range, normed, weights)
    279         # Sample is a sequence of 1D arrays.
    280         sample = atleast_2d(sample).T
--> 281         N, D = sample.shape
    282 
    283     nbin = empty(D, int)

ValueError: too many values to unpack

这是我想要完成的事情：

我有 2 个数组。一个数组来自表示土地覆盖类的地理数据集（栅格）（例如，1=树、2=草、3=建筑物等）。另一个数组来自表示某种政治边界（例如地块、人口普查区、城镇等）的地理数据集（栅格）。我正在尝试获取一个表格，其中将每个唯一的政治边界区域（数组值表示唯一的 id）列为行，并将每个土地覆盖类的每个边界内的像素总数列为列。

4

3 回答 3

2

我假设values是土地覆盖，zones是政治边界。您可能想要使用np.bincount，这就像一个特殊的直方图，其中每个 bin 的间距和宽度正好是 1。

import numpy as np

zones = np.zeros((20,30), int)
zones[:8,:15] = 100
zones[8:,:15] = 101
zones[:8,15:] = 102
zones[8:,15:] = 103

values = np.random.randint(0,10,(20,30))  # no need for that loop

tab = np.array([np.bincount(values[zones==zone]) for zone in np.unique(zones)])

但是，如果您小心处理 bin 边缘，则可以使用直方图更简单地执行此操作：

np.histogram2d(zones.flatten(), values.flatten(), bins=[np.unique(zones).size, values.max()-values.min()+1])

其工作方式如下。最简单的示例是查看所有值，而不考虑区域：

np.bincount(values)

这为您提供了一行，其中包含每个值（0 到 10）的计数。下一步是查看区域。对于一个区域，您只有一排，它将是：

zone = 101         # the desired zone
mask = zone==zones # a mask that is True wherever your zones map matches the desired zone
np.bincount(values[mask])  # count the values where the mask is True

现在，我们只想为地图中的每个区域执行此操作。您可以使用以下方式获取区域地图中唯一值的列表

zs = np.unique(zones)

并使用列表理解循环遍历它，其中每个项目是上述行之一：

tab = np.array([np.bincount(values[zones==zone]) for zone in np.unique(zones)])

然后，您的表格如下所示：

print tab
# elements with cover = 
#  0  1  2  3  4  5  6  7  8  9     # in zone:
[[16 11 10 12 13 15 11  7 13 12]    # 100
 [13 23 15 16 24 16 24 21 15 13]    # 101
 [10 12 23 13 12 11 11  5 11 12]    # 102
 [19 25 20 12 16 19 13 18 22 16]]   # 103

最后，您可以在 matplotlib 中这样绘制：

import matplotlib.pyplot as plt

plt.hist2d(zones.flatten(), values.flatten(), bins=[np.unique(zones).size, values.max()-values.min()+1])
于 2013-10-07T22:05:56.530 回答
2

histogram2d 需要一维数组作为输入，并且您的区域和值是二维的。您可以使用 ravel 将它们线性化：

np.histogram2d(zones.ravel(), values.ravel())
于 2013-10-07T22:15:16.660 回答
1

如果效率不是问题，我认为这适用于您想要做的事情

from collections import Counter

c = Counter(zip(zones.flat[:], landcover_classes.flat[:]))

c 将包含键/值元组，其中键是（区域，土地覆盖类）的元组。如果你喜欢，你可以填充一个数组

for (i, j), count in c.items():
    my_table[i, j] = count

当然，这只有在 i 和 j 是从零开始的连续整数（即从 0 到 Ni 以及从 0 到 Nj）时才有效。

于 2013-10-07T23:35:16.220 回答