加速迭代 Numpy 数组答案

【问题标题】：Speeding up iterating over Numpy Arrays加速迭代 Numpy 数组
【发布时间】：2011-10-03 00:24:46
【问题描述】：

我正在使用 Numpy 执行图像处理，特别是运行标准偏差拉伸。这会读入 X 列，找到 Std。并执行百分比线性拉伸。然后它迭代到下一个“组”列并执行相同的操作。输入图像是一个 1GB 的 32 位单波段栅格，需要相当长的时间来处理（小时）。下面是代码。

我意识到我有 3 个嵌套的 for 循环，这可能是瓶颈发生的地方。如果我在“框”中处理图像，也就是说加载一个 [500,500] 的数组并迭代图像处理时间非常短。不幸的是，相机错误要求我迭代极长的条带 (52,000 x 4) (y,x) 以避免条带。

任何关于加快速度的建议将不胜感激：

def box(dataset, outdataset, sampleSize, n):

    quiet = 0
    sample = sampleSize
    #iterate over all of the bands
    for j in xrange(1, dataset.RasterCount + 1): #1 based counter

        band = dataset.GetRasterBand(j)
        NDV = band.GetNoDataValue()

        print "Processing band: " + str(j)       

        #define the interval at which blocks are created
        intervalY = int(band.YSize/1)    
        intervalX = int(band.XSize/2000) #to be changed to sampleSize when working

        #iterate through the rows
        scanBlockCounter = 0

        for i in xrange(0,band.YSize,intervalY):

            #If the next i is going to fail due to the edge of the image/array
            if i + (intervalY*2) < band.YSize:
                numberRows = intervalY
            else:
                numberRows = band.YSize - i

            for h in xrange(0,band.XSize, intervalX):

                if h + (intervalX*2) < band.XSize:
                    numberColumns = intervalX
                else:
                    numberColumns = band.XSize - h

                scanBlock = band.ReadAsArray(h,i,numberColumns, numberRows).astype(numpy.float)

                standardDeviation = numpy.std(scanBlock)
                mean = numpy.mean(scanBlock)

                newMin = mean - (standardDeviation * n)
                newMax = mean + (standardDeviation * n)

                outputBlock = ((scanBlock - newMin)/(newMax-newMin))*255
                outRaster = outdataset.GetRasterBand(j).WriteArray(outputBlock,h,i)#array, xOffset, yOffset


                scanBlockCounter = scanBlockCounter + 1
                #print str(scanBlockCounter) + ": " + str(scanBlock.shape) + str(h)+ ", " + str(intervalX)
                if numberColumns == band.XSize - h:
                    break

                #update progress line
                if not quiet:
                    gdal.TermProgress_nocb( (float(h+1) / band.YSize) )

这里有一个更新：在不使用配置文件模块的情况下，由于我不想开始将一小部分代码包装到函数中，因此我混合使用了 print 和 exit 语句来大致了解哪些行花费的时间最多。幸运的是（我确实知道我是多么幸运）一条线拖累了一切。

    outRaster = outdataset.GetRasterBand(j).WriteArray(outputBlock,h,i)#array, xOffset, yOffset

在打开输出文件并写出数组时，GDAL 似乎效率很低。考虑到这一点，我决定将修改后的数组“outBlock”添加到 python 列表中，然后写出块。这是我更改的部分：

刚刚修改了 outputBlock ...

         #Add the array to a list (tuple)
            outputArrayList.append(outputBlock)

            #Check the interval counter and if it is "time" write out the array
            if len(outputArrayList) >= (intervalX * writeSize) or finisher == 1:

                #Convert the tuple to a numpy array.  Here we horizontally stack the tuple of arrays.
                stacked = numpy.hstack(outputArrayList)

                #Write out the array
                outRaster = outdataset.GetRasterBand(j).WriteArray(stacked,xOffset,i)#array, xOffset, yOffset
                xOffset = xOffset + (intervalX*(intervalX * writeSize))

                #Cleanup to conserve memory
                outputArrayList = list()
                stacked = None
                finisher=0

Finisher 只是一个处理边缘的标志。花了一些时间来弄清楚如何从列表中构建一个数组。在那，使用 numpy.array 正在创建一个 3-d 数组（有人愿意解释为什么吗？）并且写入数组需要一个 2d 数组。总处理时间现在从不到 2 分钟到 5 分钟不等。知道为什么可能存在时间范围吗？

非常感谢所有发帖的人！下一步是真正进入 Numpy 并了解向量化以进行额外优化。

【问题讨论】：

您是否进行了分析以找出您的热点在哪里？我可以想象您受到文件 IO 限制并且应该以更大的块从磁盘中提取数据的情况。同样，您可能内存不足，应注意创建不必要的副本。您甚至可能会受到计算限制，应该考虑更好的算法。
你能解释一下'band'是什么类型的对象吗？我同意 matt 的观点——您需要分析您的代码以确定哪些部分会减慢您的速度。
你有多少内存？ 1GB 并不是那么大的阵列。您应该能够将整个内容加载到现代机器上的内存中。像（我认为？）你想要做的对比拉伸可以就地完成。（例如，data -= whatever 和 data /= whatever 将在整个数组上逐元素操作而不进行复制）。
@Luke band 是 GDALRasterBandShadow 类型的 SwigObject。基本上是 RGB 图像中的 R，尽管这些是单波段黑白。 @matt 瓶颈发生在“outraster =”@Joe Kington 写回磁盘我有 4GB 并且可以加载整个图像。这只是一个测试，因为该程序将处理 50GB 以上的全局马赛克。这就是为什么我试图避免读取整个数组。
如果瓶颈真的是 I/O，那么你几乎别无选择；优化其余代码无济于事。

标签： python for-loop numpy gdal

【解决方案1】：

按要求回答。

如果您受 IO 限制，您应该分块您的读/写。尝试将约 500 MB 的数据转储到 ndarray，处理所有数据，将其写出，然后抓取下一个约 500 MB。确保重复使用 ndarray。

【讨论】：

【解决方案2】：

加快对numpy 数据的操作的一种方法是使用vectorize。本质上，vectorize 接受一个函数f 并创建一个新函数g 将f 映射到一个数组a 上。 g 然后被称为：g(a)。

>>> sqrt_vec = numpy.vectorize(lambda x: x ** 0.5)
>>> sqrt_vec(numpy.arange(10))
array([ 0.        ,  1.        ,  1.41421356,  1.73205081,  2.        ,
        2.23606798,  2.44948974,  2.64575131,  2.82842712,  3.        ])

如果没有您正在使用的数据，我无法确定这是否会有所帮助，但也许您可以将上述内容重写为一组可以是vectorized 的函数。也许在这种情况下，您可以将一组索引向量化为ReadAsArray(h,i,numberColumns, numberRows)。以下是潜在好处的示例：

>>> print setup1
import numpy
sqrt_vec = numpy.vectorize(lambda x: x ** 0.5)
>>> print setup2
import numpy
def sqrt_vec(a):
    r = numpy.zeros(len(a))
    for i in xrange(len(a)):
        r[i] = a[i] ** 0.5
    return r
>>> timeit.timeit(stmt='a = sqrt_vec(numpy.arange(1000000))', setup=setup1, number=1)
0.30318188667297363
>>> timeit.timeit(stmt='a = sqrt_vec(numpy.arange(1000000))', setup=setup2, number=1)
4.5400981903076172

15 倍加速！另请注意，numpy 切片可以优雅地处理ndarrays 的边缘：

>>> a = numpy.arange(25).reshape((5, 5))
>>> a[3:7, 3:7]
array([[18, 19],
       [23, 24]])

因此，如果您可以将您的 ReadAsArray 数据转换为 ndarray，您就不必进行任何边缘检查恶作剧了。

关于您关于重塑的问题 - 重塑根本不会改变数据。它只是改变了numpy 索引数据的“步幅”。当您调用reshape 方法时，返回的值是数据的新视图；数据根本不会被复制或更改，具有旧步幅信息的旧视图也不会。

>>> a = numpy.arange(25)
>>> b = a.reshape((5, 5))
>>> a
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24])
>>> b
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])
>>> a[5]
5
>>> b[1][0]
5
>>> a[5] = 4792
>>> b[1][0]
4792
>>> a.strides
(8,)
>>> b.strides
(40, 8)

【讨论】：

感谢您的信息。我现在正在阅读矢量化文档并致力于修改我的代码。 readAsArray 确实返回了一个 ndarray。我遇到的问题是数组的大小可能大于 3gb（这些是巨大的图像），我需要在摄取它们之前对它们进行分段。我相信我看到的瓶颈实际上是逐块写出新数组。我正在尝试写入元组并减少写入次数。结合矢量化，这可能会奏效！
我搜索的另一个问题。为了很好地处理边缘，我需要重塑吗？我犹豫不决，因为我正在处理卫星图像，并试图在拉伸（基本上）反照率时利用空间自相关。
@Jzl5325，请参阅我关于重塑的编辑。简而言之，重塑根本不会改变数据。
来自 numpy 文档：“提供 vectorize 函数主要是为了方便，而不是为了性能。实现本质上是一个 for 循环。” numpy.org/doc/stable/reference/generated/numpy.vectorize.html 换句话说，它不依赖于优化的、预编译的 C 函数，因此没有运行时优势
@Ender，提供vectorize 主要是为了方便这一事实并不能得出没有运行时优势。请再看一下我包含的测试。有一个数量级的加速。当然，如果您可以使用内置的 numpy 操作，您可能会获得三个数量级的加速。但如果你不能，vectorize 有时仍然可以提供一点帮助。

【解决方案3】：

没有试图完全理解你在做什么，我注意到你没有使用任何numpy slices 或array broadcasting，这两者都可以加速你的代码，或者至少，让它更多可读。如果这些与您的问题无关，我深表歉意。

【讨论】：

我曾考虑过索引。从文档看来，我需要将整个数组加载到内存中才能做到这一点？使用 band.ReadAsArray - GDAL 函数。
好吧。尝试重新考虑您的问题并将其写为矩阵运算。如果你能做到这一点，那么你可以很容易地用 numpy 将它应用到 python 中（你的矩阵将是一个二维数组等等，......）。如果您遇到内存问题，请尝试在子矩阵或其他矩阵表示中进行处理。抱歉，除了这几个提示，我无能为力。我无法真正进入你的 sn-p ;)