Optimizando Python GDAL ReadAsArray

Estoy usando el método GDAL ReadAsArray para trabajar con datos ráster usando numpy (específicamente reclasificación). Como mis rásteres son grandes, proceso las matrices en bloques, iterando en cada bloque y procesando con un método similar al ejemplo de GeoExamples .

Ahora estoy buscando la mejor manera de establecer el tamaño de estos bloques para optimizar el tiempo necesario para procesar todo el ráster. Conociendo las limitaciones con los tamaños de matriz numpy, y el uso de GDAL GetBlockSize para usar el tamaño de bloque "natural" de un ráster, he realizado pruebas usando algunos tamaños de bloque diferentes, formados por múltiplos del tamaño "natural", con el siguiente código de ejemplo:

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)

Lo que produce el siguiente tipo de salida:

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

He intentado ejecutar esto para algunos rásteres diferentes, con diferentes tamaños y tipos de píxeles, y parece que están obteniendo tendencias similares, donde un aumento de diez veces en la dimensión x o y (en algunos casos, ambos) reduce a la mitad el tiempo de procesamiento, que aunque no es tan significativo en el ejemplo anterior, puede significar una cantidad de minutos para mis rásters más grandes.

Entonces mi pregunta es, ¿por qué ocurre este comportamiento?

Esperaba usar menos bloques para mejorar el tiempo de procesamiento, pero las pruebas que usan menos no son las más rápidas. Además, ¿por qué la prueba final tarda mucho más que la anterior? ¿Hay algún tipo de preferencia con los rásteres para leer por fila o columna, o en la forma del bloque que se lee, el tamaño total? Lo que espero obtener de esto es la información para reunir un algoritmo básico que pueda establecer el tamaño de bloque de un ráster en un valor óptimo, dependiendo del tamaño de la entrada.

Tenga en cuenta que mi entrada es un ráster de cuadrícula ArcRIFO de ESRI, que tiene un tamaño de bloque "natural" de 256 x 16, y el tamaño total de mi ráster en este ejemplo fue 41740 x 46280.

¿Has intentado usar un tamaño de bloque igual? Trato con datos ráster que son del orden de 200k x 200k píxeles y bastante dispersos. Una gran cantidad de benchmarking ha generado bloques de 256x256 píxeles como los más eficientes para nuestros procesos. Todo esto tiene que ver con la cantidad de búsquedas de disco que se requieren para recuperar un bloque. Si el bloque es demasiado grande, entonces es más difícil escribirlo en el disco contiguamente, lo que significa más búsquedas. Del mismo modo, si es demasiado pequeño, deberá realizar muchas lecturas para procesar todo el ráster. También ayuda a garantizar que el tamaño total sea una potencia de dos. Por cierto, 256x256 es el tamaño de bloque geotiff predeterminado en gdal, por lo que quizás sacaron la misma conclusión

Mi sospecha es que realmente estás chocando con el bloque de caché de GDAL, y esa es una perilla que tendrá un impacto significativo en tu curva de velocidad de procesamiento.

Vea SettingConfigOptions , específicamente GDAL_CACHEMAX, para más detalles sobre esto e investigue cómo cambiar ese valor a algo significativamente mayor interactúa con su simulación.