Suite

Optimiser Python GDAL ReadAsArray


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 par blocs, en itérant à travers chaque bloc et en traitant avec une méthode similaire à l'exemple GeoExamples.

Je cherche maintenant à définir au mieux la taille de ces blocs pour optimiser le temps de traitement de l'ensemble du raster. Étant 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, composées de multiples de la taille "naturelle", avec l'exemple de code ci-dessous :

import timeit try: import gdal except: from osgeo import gdal # Fonction permettant de lire le raster sous forme de tableaux pour la taille de bloc choisie. def read_raster(x_block_size, y_block_size): raster = "chemin vers un grand raster" ds = gdal.Open(raster) band = ds.GetRasterBand(1) xsize = band.XSize ysize = band.YSize blocks = 0 pour y dans 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) # Fonction pour exécuter le test et imprimer le temps nécessaire pour terminer. 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 "	{:. 2f}s
".format(t.timeit(1)) raster = "chemin vers un grand raster" ds = gdal.Open(raster) band = ds.GetRasterBand(1) # Obtenir la taille de bloc "naturelle" et le total taille raster XY. block_sizes = band.GetBlockSize() x_block_size = block_sizes[0] y_block_size = block_sizes[1] xsize = band.XSize ysize = band.YSize band = Aucun ds = Aucun # Tests avec différentes tailles de blocs. 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_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 blocs taille 256 x 16 : 9,12s 4930 blocs taille 2560 x 160 : 5,32s 58 blocs taille 25600 x 1600 : 5,72s 49181 blocs taille 2560 x 16 : 4,22s 5786 blocs taille 25600 x 16 : 5,67s 47560 blocs taille 256 x 160 : 4,21s 4756 blocs taille 256 x 1600 : 5,62s 2893 blocs taille 41740 x 16 : 5,85s 164 blocs taille 256 x 46280 : 5,97s 46280 blocs taille 41740 x 1 : 5,00s 41740 blocs taille 1 x 46280 : 800,24s

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

Ma question est donc : 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 ? Y a-t-il une sorte de préférence avec les rasters pour la lecture par ligne ou colonne, ou dans la forme du bloc en cours de lecture, la taille totale ? Ce que j'espère en tirer, ce sont les informations pour obtenir un algorithme de base qui sera capable 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 ESRI ArcINFO, qui a une taille de bloc "naturelle" de 256 x 16, et la taille totale de mon raster dans cet exemple était de 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 éparses. De nombreuses analyses comparatives ont révélé que les blocs de 256x256 pixels sont les plus efficaces pour nos processus. Tout cela a à voir avec le 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 de nombreuses lectures pour traiter l'ensemble du raster. Cela permet également de s'assurer que la taille totale est une puissance de deux. 256x256 est d'ailleurs la taille de bloc geotiff 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 particulierGDAL_CACHEMAX, pour plus de détails à ce sujet et étudiez comment le changement de cette valeur en quelque chose de beaucoup plus important interagit avec votre simulation.


Voir la vidéo: Read and write raster files with GDAL in Python (Octobre 2021).