Mengoptimalkan Python GDAL ReadAsArray

9

Saya menggunakan metode ReadalArray GDAL untuk bekerja dengan data raster menggunakan numpy (khusus klasifikasi ulang). Karena raster saya besar, saya memproses array dalam blok, iterasi setiap blok dan pemrosesan dengan metode yang mirip dengan contoh GeoExamples .

Saya sekarang mencari cara terbaik untuk mengatur ukuran blok ini untuk mengoptimalkan waktu yang dibutuhkan untuk memproses seluruh raster. Menyadari keterbatasan dengan ukuran array numpy, dan penggunaan GDAL GetBlockSize untuk menggunakan ukuran blok "natural" dari sebuah raster, saya telah menguji menggunakan beberapa ukuran blok yang berbeda, terdiri dari beberapa ukuran "natural", dengan contoh kode di bawah ini:

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)

Yang menghasilkan semacam output berikut:

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

Saya telah mencoba menjalankan ini untuk beberapa raster yang berbeda, dengan ukuran dan tipe piksel yang berbeda, dan tampaknya mendapatkan tren yang sama, di mana peningkatan sepuluh kali lipat dalam dimensi x atau y (dalam beberapa kasus, keduanya) membagi dua waktu pemrosesan, yang walaupun tidak begitu penting dalam contoh di atas, bisa berarti beberapa menit untuk raster terbesar saya.

Jadi pertanyaan saya adalah, mengapa perilaku ini terjadi?

Saya memang berharap menggunakan blok lebih sedikit untuk meningkatkan waktu pemrosesan, tetapi tes menggunakan paling tidak bukan yang tercepat. Juga, mengapa ujian akhir memakan waktu lebih lama dari yang sebelumnya? Apakah ada semacam preferensi dengan raster untuk dibaca oleh baris atau kolom, atau dalam bentuk blok yang dibaca, ukuran total? Apa yang saya harapkan dari ini adalah informasi untuk mendapatkan algoritma dasar bersama yang akan dapat mengatur ukuran blok raster ke nilai optimal, tergantung pada ukuran input.

Perhatikan bahwa input saya adalah raster kisi ESRI ArcINFO, yang memiliki ukuran blok "alami" 256 x 16, dan ukuran total raster saya dalam contoh ini adalah 41740 x 46280.

ssast
sumber
Tes luar biasa! banyak membantu! Bersulang!
Isaque Daniel

Jawaban:

4

Sudahkah Anda mencoba menggunakan ukuran blok yang sama. Saya berurusan dengan data raster yang berada di urutan 200k x 200k piksel dan cukup jarang. Banyak pembandingan yang menghasilkan blok 256x256 piksel sebagai yang paling efisien untuk proses kami. Ini semua berkaitan dengan berapa banyak disk yang diperlukan untuk mengambil blok. Jika blok terlalu besar maka lebih sulit untuk menulisnya ke disk secara bersebelahan, artinya mencari lebih banyak. Demikian juga, jika terlalu kecil, Anda harus melakukan banyak pembacaan untuk memproses seluruh raster. Ini juga membantu memastikan ukuran total adalah kekuatan dua. 256x256 adalah kebetulan ukuran blok geotiff default di gdal, jadi mungkin mereka menarik kesimpulan yang sama

James
sumber
Blok 256 x 256 sedikit lebih cepat daripada kebanyakan tes lainnya (dan sama dengan 2560 x 16 dan 41740 x 1), tetapi hanya sekitar 5%. Namun, dengan mengonversi raster saya ke format geotiff, itu adalah opsi tercepat setidaknya 20%, jadi untuk tiffs setidaknya yang terlihat menjadi pilihan ukuran blok yang bagus. Gdal saya memiliki ukuran default geotiff pada 128 x 128.
ssast
1
Ya, jika Anda memiliki pilihan format, geotiff adalah pilihan terbaik - sejauh ini waktu pengembangan yang paling banyak telah dimasukkan ke driver ini. Juga bereksperimen dengan kompresi, dan jika data Anda jarang (banyak nilai nol), Anda harus melihat menggunakan opsi pembuatan SPARSE_OK dan melewati membaca / menulis blok nol
James
Baik untuk mengetahui referensi di masa mendatang, meskipun saya terjebak dengan membaca grid ESRI ArcINFO untuk contoh yang diberikan.
ssast
Juga, untuk perbedaan antara dua contoh terakhir, Anda akan ingin membaca tentang urutan utama baris vs urutan utama kolom. Sekali lagi, tergantung pada berapa banyak disk yang diperlukan untuk membangun blok yang diminta.
James
1

Kecurigaan saya adalah Anda benar-benar bertemu dengan cache block GDAL, dan itu adalah tombol yang akan memiliki dampak signifikan pada kurva kecepatan pemrosesan Anda.

Lihat SettingConfigOptions , secara khusus GDAL_CACHEMAX, untuk detail lebih lanjut tentang ini dan menyelidiki bagaimana mengubah nilai itu menjadi sesuatu yang secara signifikan lebih besar berinteraksi dengan simulasi Anda.

Howard Butler
sumber
Setelah mengatur cache ke 1 GB dengan gdal.SetCacheMax (1000000000), ada penurunan beberapa detik di sebagian besar tes, kecuali yang terakhir 1 x ukuran satu, yang mempercepat hingga 40-an. Mengurangi cache ke 1 mb sebenarnya mempercepat semua tes kecuali dua yang terakhir. Saya pikir ini karena saya menggunakan ukuran blok alami untuk sebagian besar tes, dan tidak memiliki blok yang tumpang tindih, sehingga cache tidak diperlukan kecuali untuk dua tes terakhir.
ssast