Otimizando o Python GDAL ReadAsArray

Estou usando o método GDAL ReadAsArray para trabalhar com dados rasterizados usando numpy (especificamente reclassificação). Como minhas rasters são grandes, eu processo as matrizes em blocos, repetindo cada bloco e processando com um método semelhante ao exemplo de GeoExamples .

Agora, estou analisando a melhor forma de definir o tamanho desses blocos para otimizar o tempo necessário para processar toda a varredura. Consciente das limitações dos tamanhos de matriz numpy e do uso do GDAL GetBlockSize para usar o tamanho de bloco "natural" de uma varredura, tenho testes usando alguns tamanhos de bloco diferentes, compostos por múltiplos do tamanho "natural", com o código de exemplo abaixo:

import timeit
try:
    import gdal
except:
    from osgeo import gdal

# Function to read the raster as arrays for the chosen block size.
def read_raster(x_block_size, y_block_size):
    raster = "path to large raster"
    ds = gdal.Open(raster)
    band = ds.GetRasterBand(1)
    xsize = band.XSize
    ysize = band.YSize
    blocks = 0
    for y in xrange(0, ysize, y_block_size):
        if y + y_block_size < ysize:
            rows = y_block_size
        else:
            rows = ysize - y
        for x in xrange(0, xsize, x_block_size):
            if x + x_block_size < xsize:
                cols = x_block_size
            else:
                cols = xsize - x
            array = band.ReadAsArray(x, y, cols, rows)
            del array
            blocks += 1
    band = None
    ds = None
    print "{0} blocks size {1} x {2}:".format(blocks, x_block_size, y_block_size)

# Function to run the test and print the time taken to complete.
def timer(x_block_size, y_block_size):
    t = timeit.Timer("read_raster({0}, {1})".format(x_block_size, y_block_size),
                     setup="from __main__ import read_raster")
    print "\t{:.2f}s\n".format(t.timeit(1))

raster = "path to large raster"
ds = gdal.Open(raster)
band = ds.GetRasterBand(1)

# Get "natural" block size, and total raster XY size. 
block_sizes = band.GetBlockSize()
x_block_size = block_sizes[0]
y_block_size = block_sizes[1]
xsize = band.XSize
ysize = band.YSize
band = None
ds = None

# Tests with different block sizes.
timer(x_block_size, y_block_size)
timer(x_block_size*10, y_block_size*10)
timer(x_block_size*100, y_block_size*100)
timer(x_block_size*10, y_block_size)
timer(x_block_size*100, y_block_size)
timer(x_block_size, y_block_size*10)
timer(x_block_size, y_block_size*100)
timer(xsize, y_block_size)
timer(x_block_size, ysize)
timer(xsize, 1)
timer(1, ysize)

Que produz o seguinte tipo de saída:

474452 blocks size 256 x 16:
        9.12s

4930 blocks size 2560 x 160:
        5.32s

58 blocks size 25600 x 1600:
        5.72s

49181 blocks size 2560 x 16:
        4.22s

5786 blocks size 25600 x 16:
        5.67s

47560 blocks size 256 x 160:
        4.21s

4756 blocks size 256 x 1600:
        5.62s

2893 blocks size 41740 x 16:
        5.85s

164 blocks size 256 x 46280:
        5.97s

46280 blocks size 41740 x 1:
        5.00s

41740 blocks size 1 x 46280:
        800.24s

Tentei executá-lo em algumas classes diferentes, com tamanhos e tipos de pixels diferentes, e parece estar obtendo tendências semelhantes, nas quais um aumento de dez vezes na dimensão x ou y (em alguns casos, ambos) reduz pela metade o tempo de processamento, o que embora não seja tão significativo no exemplo acima, pode significar alguns minutos para meus maiores rasters.

Então, minha pergunta é: por que esse comportamento está ocorrendo?

Eu esperava usar menos blocos para melhorar o tempo de processamento, mas os testes que usam menos não são os mais rápidos. Além disso, por que o teste final leva muito mais tempo do que o anterior? Existe algum tipo de preferência com rasters para ler por linha ou coluna, ou na forma do bloco que está sendo lido, o tamanho total? O que espero obter disso é a informação para reunir um algoritmo básico capaz de definir o tamanho do bloco de uma varredura para um valor ideal, dependendo do tamanho da entrada.

Observe que minha entrada é uma varredura de grade ESRI ArcINFO, que possui um tamanho de bloco "natural" de 256 x 16, e o tamanho total da minha varredura neste exemplo foi 41740 x 46280.

Você já tentou usar um tamanho de bloco igual. Eu lido com dados raster, da ordem de 200k x 200k pixels e bastante esparsos. Muitos testes de referência produziram blocos de 256x256 pixels como os mais eficientes para nossos processos. Isso tem tudo a ver com quantas pesquisas de disco são necessárias para recuperar um bloco. Se o bloco for muito grande, será mais difícil gravá-lo no disco de forma contígua, o que significa mais pesquisas. Da mesma forma, se for muito pequeno, você precisará fazer muitas leituras para processar toda a varredura. Também ajuda a garantir que o tamanho total seja uma potência de dois. Aliás, 256x256 é o tamanho padrão do bloco geotiff em gdal, então talvez eles tenham chegado à mesma conclusão

Suspeito que você esteja realmente colidindo com o cache de blocos do GDAL, e esse é um botão que terá um impacto significativo na sua curva de velocidade de processamento.

Consulte SettingConfigOptions , especificamente GDAL_CACHEMAX, para obter mais detalhes sobre isso e investigar como alterar esse valor para algo significativamente maior interage com sua simulação.