Optimizando Python GDAL ReadAsArray

9

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.

ssast
fuente
Prueba increíble! ayuda mucho! ¡salud!
Isaque Daniel

Respuestas:

4

¿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

James
fuente
Los bloques de 256 x 256 fueron un poco más rápidos que la mayoría de las otras pruebas (e iguales a 2560 x 16 y 41740 x 1), pero solo en aproximadamente un 5%. Sin embargo, al convertir mi ráster a formato geotiff, fue la opción más rápida en al menos un 20%, por lo que para tiffs al menos eso parece ser una buena opción de tamaño de bloque. Sin embargo, mi gdal tenía el tamaño de bloque geotiff predeterminado en 128 x 128.
ssast
1
Sí, si tiene la opción de formatos, geotiff es la mejor opción, con mucho, el mayor tiempo de desarrollo se ha dedicado a este controlador. También experimente con la compresión, y si sus datos son escasos (muchos valores nulos), debería considerar usar la opción de creación SPARSE_OK y omitir la lectura / escritura de bloques nulos
James
Es bueno saberlo para referencia futura, aunque estoy atascado con la lectura de cuadrículas ArcRIFO de ESRI para el ejemplo dado.
ssast
Además, por la diferencia entre los dos ejemplos finales, querrá leer sobre el orden principal de la fila frente al orden principal de la columna. Nuevamente, se trata de cuántas búsquedas de disco se requieren para construir el bloque solicitado.
James
1

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.

Howard Butler
fuente
Después de configurar el caché a 1 GB con gdal.SetCacheMax (1000000000), hubo una disminución de algunos segundos en la mayoría de las pruebas, excepto el último 1 x tamaño grande, que aceleró hasta 40 segundos. Disminuir el caché a 1 mb en realidad acelera todas las pruebas, excepto las dos últimas. Creo que esto se debe a que estoy usando el tamaño de bloque natural para la mayoría de las pruebas, y no tengo bloques superpuestos, por lo que la caché no es necesaria, excepto para las dos pruebas finales.
ssast