我正在使用Numpy执行图像处理，特别是运行标准差拉伸。这读入X列，找到标准差，并执行百分比线性拉伸。然后迭代到下一个列“组”，并执行相同的操作。输入图像是一个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) )

以下是最新情况：在没有使用profile模块的情况下，因为我不想把代码的一小部分打包成函数，所以我使用了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创建了一个三维数组（有人愿意解释一下为什么吗？）而write array需要一个二维数组。总处理时间现在从不到2分钟到5分钟不等。你知道为什么会存在这个时间范围吗？
非常感谢每个发帖的人！下一步是真正进入Numpy，学习矢量化以进行额外的优化。