Optimieren von Python GDAL ReadAsArray

Ich verwende die GDAL ReadAsArray-Methode, um mit Rasterdaten unter Verwendung von numpy (insbesondere Neuklassifizierung) zu arbeiten. Da meine Raster groß sind, verarbeite ich die Arrays in Blöcken, durchlaufe jeden Block und verarbeite sie mit einer ähnlichen Methode wie im GeoExamples- Beispiel.

Ich schaue mir jetzt an, wie die Größe dieser Blöcke am besten eingestellt werden kann, um die Zeit zu optimieren, die für die Verarbeitung des gesamten Rasters benötigt wird. Da ich mir der Einschränkungen bei numpy-Array-Größen und der Verwendung von GDAL GetBlockSize zur Verwendung der "natürlichen" Blockgröße eines Rasters bewusst bin, habe ich einige verschiedene Blockgrößen getestet, die aus Vielfachen der "natürlichen" Größe bestehen. mit dem folgenden Beispielcode:

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)

Welches erzeugt die folgende Art von Ausgabe:

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

Ich habe versucht, dies für einige verschiedene Raster mit unterschiedlichen Größen und Pixeltypen auszuführen, und es scheint ähnliche Trends zu geben, bei denen eine zehnfache Vergrößerung der x- oder y-Dimension (in einigen Fällen beides) die Verarbeitungszeit halbiert Obwohl im obigen Beispiel nicht so wichtig, kann dies für meine größten Raster eine Anzahl von Minuten bedeuten.

Meine Frage ist also, warum dieses Verhalten auftritt.

Ich hatte erwartet, weniger Blöcke zu verwenden, um die Verarbeitungszeit zu verbessern, aber die Tests mit den wenigsten sind nicht die schnellsten. Warum dauert der endgültige Test so viel länger als der vorhergehende? Gibt es eine Präferenz bei Rastern für das Lesen nach Zeilen oder Spalten oder in der Form des zu lesenden Blocks die Gesamtgröße? Was ich davon erhoffe, sind die Informationen, um einen grundlegenden Algorithmus zusammenzustellen, mit dem die Blockgröße eines Rasters abhängig von der Größe der Eingabe auf einen optimalen Wert eingestellt werden kann.

Beachten Sie, dass meine Eingabe ein ESRI ArcINFO-Raster ist, das eine "natürliche" Blockgröße von 256 x 16 hat, und die Gesamtgröße meines Rasters in diesem Beispiel 41740 x 46280 betrug.

Haben Sie versucht, eine gleiche Blockgröße zu verwenden? Ich beschäftige mich mit Rasterdaten, die in der Größenordnung von 200.000 x 200.000 Pixel liegen und recht spärlich sind. Viele Benchmarking-Tests haben Blöcke mit 256 x 256 Pixel ergeben, die für unsere Prozesse am effizientesten sind. Dies hängt alles damit zusammen, wie viele Festplattensuchen erforderlich sind, um einen Block abzurufen. Wenn der Block zu groß ist, ist es schwieriger, ihn zusammenhängend auf die Festplatte zu schreiben, was mehr Suchvorgänge bedeutet. Wenn es zu klein ist, müssen Sie ebenfalls viele Lesevorgänge ausführen, um das gesamte Raster zu verarbeiten. Es hilft auch sicherzustellen, dass die Gesamtgröße eine Zweierpotenz ist. 256x256 ist im Übrigen die Standard-Geotiff-Blockgröße in GDAL, daher haben sie möglicherweise die gleiche Schlussfolgerung gezogen

Mein Verdacht ist, dass Sie wirklich auf den Block-Cache von GDAL stoßen, und das ist ein Regler, der einen erheblichen Einfluss auf Ihre Verarbeitungsgeschwindigkeitskurve haben wird.

Siehe SettingConfigOptions , speziell GDAL_CACHEMAX, für weitere Einzelheiten zu diesen und untersuchen , wie sich dieser Wert auf etwas wesentlich größer wirkt mit Ihrer Simulation ändern.