有没有比使用 np.where 更快的迭代一个非常大的 2D numpy 数组的方法？

Question

我有一个巨大的二维 numpy 数组，其中填充了整数值。我通过 gdal.GetRasterBand() 从 .tif-image 收集它们。图像的像素值表示唯一的集群标识号。因此，一个簇内的所有像素都具有相同的值。在我的脚本中，我想检查簇的像素是否超过特定阈值。如果集群大小大于阈值，我想保留集群并给它们一个像素值 1。如果集群的像素小于阈值，则该集群的所有像素都应获得值 0。

我的代码目前有效，但非常非常慢。而且因为我想改变阈值，所以需要永远。我将衷心感谢您的帮助。谢谢。

# Import GeoTIFF via GDAL and convert to NumpyArray
data = gdal.Open(image)
raster = data.GetRasterBand(1)
raster = raster.ReadAsArray()

# Different thresholds for iteration
thresh = [0,10,25,50,100,1000,2000]

for threshold in thresh:
        clusteredRaster = np.array(raster.copy(), dtype = int)

        for clump in np.unique(clusteredRaster): # Unique ids of the clusters in image

            if clusteredRaster[np.where(clusteredRaster == clump)].size >= threshold: 
                clusteredRaster[np.where(clusteredRaster == clump)] = int(1)

            else:
                clusteredRaster[np.where(clusteredRaster == clump)] = int(0)
'''

[ClusterImage][1]

In the image you can see the cluster image. Each color stands vor a specific clusternumber. I want to delete the small ones (under a specific size) and just keep the big ones.

  [1]: https://i.stack.imgur.com/miEKg.png

Answer 1

可以进行一些修改来提高性能，

clusteredRaster = np.array(raster.copy(), dtype = int)

可以替换为

clusteredRaster = raster.astype(int)

这本质上既是复制运算符又是转换运算符，因此此操作速度更快。

现在

clusteredRaster[np.where(clusteredRaster == clump)] = int(1)

您不需要调用 np.where，这样会更快

clusteredRaster[clusteredRaster == clump] = int(1)

这部分也完成了

clusteredRaster[np.where(clusteredRaster == clump)].size

您还可以删除 clusteredRaster==clump 的计算两次，如下所示：

for clump in np.unique(clusteredRaster): # Unique ids of the clusters in image
    indicies = clusteredRaster==clump

    if clusteredRaster[indicies].size >= threshold: 
        clusteredRaster[indicies] = int(1)

    else:
        clusteredRaster[indicies] = int(0)

我认为你的函数现在的运行速度是原来的两倍，但是如果你想运行更快，你必须使用更小的数据类型，比如 np.uint8 而不是普通的 int，前提是你的图像是用 RGB 编码的并且可以用 8 位整数表示（或者如果 8 位太低，则可能是 np.uint16？）这是从 python 端获得的最快速度。

有一些更快的方法，例如使用带有 openmp 的 C 模块来多线程处理您的工作 across 多个内核，这可以通过 numba or cython 之类的方法轻松完成，而不必担心编写 C 代码，但是如果你想获得有史以来最快的性能，还有很多阅读要做，比如使用哪个线程后端（TBB vs openmp）和一些 os 和设备相关的功能。

Answer 2

除了之外，您还可以在 for 循环之外计算唯一值、索引和计数。另外，您不需要每次都复制 raster - 您可以制作一个 np.uint8 的数组。

这给出了与您的原始实施相同的结果：

data = gdal.Open(image)
raster = data.GetRasterBand(1).ReadAsArray()

# Different thresholds for iteration
thresh = [0, 10, 25, 50, 100, 1000, 2000]

# determine the unique clumps and their frequencies outside of the for loops
clumps, counts = np.unique(raster, return_counts=True)

# only determine the indices once, rather than for each threshold
indices = np.asarray([raster==clump for clump in clumps])

for threshold in thresh:
    
    clustered_raster = np.zeros_like(raster, dtype=np.uint8)

    for clump_indices, clump_counts in zip(indices, counts):
    
        clustered_raster[clump_indices] = clump_counts >= threshold

Answer 3

根据您的有用回答，我得到了一个简单的解决方案！这个想法是找到每个阈值的唯一值和簇大小，并立即填写正确的值，从而避免循环。它将每次迭代的迭代时间从最初的 142 秒减少到 0.52 秒并重现相同的结果。

data = gdal.Open(image)
raster = data.GetRasterBand(1).ReadAsArray()

thresh = [0, 10, 25, 50, 100, 1000, 2000]   
for threshold in thresh:
    # Create new 0-raster with same dimensions as input raster
    clusteredRaster = np.zeros(raster.shape, dtype = uint8)
    
    # Get unique cluster IDs and count the size of the occurence
    clumps, counts = np.unique(raster, return_counts=True)

    # Get only the clumps which are bigger than the threshold
    biggerClumps = clumps[counts >= threshold]

    # fill in ones for the relevant cluster IDs
    clusteredRaster[np.isin(raster,biggerClumps)] = 1

有没有比使用 np.where 更快的迭代一个非常大的 2D numpy 数组的方法？

Is there a faster method for iterating over a very big 2D numpy array than using np.where?

python

numpy

cluster-computing

gdal