Optimisation de ReadAsArray Python GDAL

J'utilise la méthode GDAL ReadAsArray pour travailler avec des données raster à l'aide de numpy (en particulier la reclassification). Comme mes rasters sont volumineux, je traite les tableaux en blocs, en répétant chaque bloc et en les traitant avec une méthode similaire à l' exemple GeoExamples .

J'examine maintenant la meilleure façon de définir la taille de ces blocs pour optimiser le temps nécessaire au traitement de l'ensemble du raster. Conscient des limites des tailles de tableau numpy et de l'utilisation de GDAL GetBlockSize pour utiliser la taille de bloc "naturelle" d'un raster, j'ai testé en utilisant quelques tailles de bloc différentes, constituées de multiples de la taille "naturelle", avec l'exemple de code ci-dessous:

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)

Ce qui produit le type de sortie suivant:

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

J'ai essayé de l'exécuter pour quelques rasters différents, avec différentes tailles et types de pixels, et semble obtenir des tendances similaires, où une augmentation de dix fois dans la dimension x ou y (dans certains cas, les deux) réduit de moitié le temps de traitement, ce qui bien que non significatif dans l'exemple ci-dessus, cela peut signifier un certain nombre de minutes pour mes plus grands rasters.

Donc ma question est, pourquoi ce comportement se produit-il?

Je m'attendais à utiliser moins de blocs pour améliorer le temps de traitement, mais les tests utilisant le moins ne sont pas les plus rapides. Aussi, pourquoi le test final prend-il autant de temps que celui qui le précède? Existe-t-il une sorte de préférence avec les rasters pour la lecture par ligne ou colonne, ou sous la forme du bloc en cours de lecture, la taille totale? Ce que j'espère obtenir de cela, ce sont les informations pour obtenir un algorithme de base qui sera en mesure de définir la taille de bloc d'un raster à une valeur optimale, en fonction de la taille de l'entrée.

Notez que mon entrée est un raster de grille ArcRIFO ESRI, qui a une taille de bloc "naturelle" de 256 x 16, et la taille totale de mon raster dans cet exemple était 41740 x 46280.

Avez-vous essayé d'utiliser une taille de bloc égale. Je traite des données raster qui sont de l'ordre de 200k x 200k pixels et assez clairsemées. De nombreuses analyses comparatives ont produit des blocs de 256 x 256 pixels comme les plus efficaces pour nos processus. Cela est lié au nombre de recherches de disque nécessaires pour récupérer un bloc. Si le bloc est trop volumineux, il est plus difficile de l'écrire sur le disque de manière contiguë, ce qui signifie plus de recherches. De même, s'il est trop petit, vous devrez effectuer plusieurs lectures pour traiter l'intégralité du raster. Cela permet également de garantir que la taille totale est une puissance de deux. 256x256 est d'ailleurs la taille de bloc géotiff par défaut dans gdal, alors peut-être qu'ils ont tiré la même conclusion

Je soupçonne que vous vous heurtez vraiment au cache de blocs de GDAL, et c'est un bouton qui va avoir un impact significatif sur votre courbe de vitesse de traitement.

Voir SettingConfigOptions , en particulier GDAL_CACHEMAX, pour plus de détails à ce sujet et examinez comment la modification de cette valeur en quelque chose de beaucoup plus grand interagit avec votre simulation.