Secara efisien membaca dan mengklasifikasi ulang banyak raster di R?

8

Saya telah ditugaskan untuk membuat analisis kesesuaian kondisi gelombang di Teluk Meksiko. Saya memiliki 2 ribu file raster yang masing-masing sekitar 8 MB (2438 kolom, 1749 baris, ukuran sel 1 km). Parameter dalam file raster adalah periode gelombang dan saya ingin mengklasifikasi ulang semua raster sehingga jika 4<= wave period <=9 then make cell = 1, sel lain = 0. Kemudian jumlah semua raster menjadi raster terakhir dan bagi dengan jumlah total raster untuk memberikan persentase total pengamatan yang sesuai dan hasil ekspor ke dalam beberapa format yang kompatibel dengan ESRI ... mungkin sesuatu yang dapat mendukung pelampung jika perlu. Saya tidak bekerja sama sekali dengan Python atau R, tetapi setelah mencari online sepertinya masuk akal untuk melakukan proses ini dalam salah satu bahasa tersebut. Sejauh ini saya telah menemukan beberapa kode di R, tetapi saya bingung bagaimana cara membuatnya.

library(rgdal)
raster_data <- list.files(path=getwd())    #promt user for dir containing raster files
num_files <- length(raster_data)
for (i in raster_data) {                   #read in rasters
   my_data <- readGDAL(raster_data[i])

Pada titik ini saya bingung apakah saya juga harus mereklasifikasi dan mulai menjumlahkan data dalam loop ini atau tidak. Dugaan saya adalah ya karena kalau tidak saya pikir saya akan kehabisan memori di komputer saya, tetapi tidak yakin. Saya juga tidak yakin tentang cara mereklasifikasi data.

Dalam riset online saya pikir saya akan menggunakan reclass(my_data, c(-Inf,3,0, 4,9,1, 10,Inf,0)), tetapi apakah itu terlihat benar?

Dan untuk menjumlahkan akan saya gunakan sum(stack(my_data))dan entah bagaimana output itu. Juga ... jika ini lebih efisien dilakukan atau ditulis dalam Python saya akan terbuka untuk itu juga. Saya benar-benar pemula dalam hal pemrograman.

Nigel
sumber
Cukup gunakan python-gdal. Ini akan jauh lebih efisien dalam kasus Anda.
SS_Rebelious
Terima kasih, Rebelious. Hanya ingin tahu mengapa python-gdal lebih efisien dalam situasi ini? Apakah mungkin untuk melihat beberapa langkah dalam kode yang perlu dilakukan? Mencoba mencari cara terbaik untuk melakukan ini dengan memanfaatkan sedikit memori, dan cpu, mungkin ... membingungkan untuk mengetahui cara menulis kode sehingga akan membaca dalam data, memproses, mengeluarkannya dari memori dan lalu pindah ke raster berikutnya.
Nigel
Saya tidak bisa memberi tahu Anda persisnya mengapa, tetapi penyebab umumnya adalah bahwa R dirancang untuk tujuan lain dan diketahui berkinerja lambat dengan siklus. Sebagai contoh kode, jika tidak ada yang memberikannya, saya akan membagikan satu dengan Anda dalam waktu sekitar 10 jam ketika saya akan mendapatkan akses ke mesin tempat skrip yang sesuai disimpan.
SS_Rebelious

Jawaban:

8

Ini adalah cara ringkas untuk melakukannya di R --- di sini tanpa file perantara:

library(raster)
raster_data <- list.files(path=getwd())    #promt user for dir containing raster files
s <- stack(raster_data)
f <- function(x) { rowSums(x >= 4 & x <= 9) }
x <- calc(s, f, progress='text', filename='output.tif')
Robert Hijmans
sumber
1
+1 Ini bagus untuk masalah kecil, tetapi mari kita lakukan perhitungannya untuk ini: 2438 kolom kali 1749 baris kali 8 byte / nilai kali 2 ribu kisi = 63,6 GB, yang semuanya Rharus disimpan dalam RAM untuk dibuat s. (Kemungkinan RAM dua kali lebih banyak diperlukan karena smungkin tidak akan diganti raster_data.) Saya harap Anda memiliki banyak RAM! Solusi Anda dapat dibuat praktis dengan memecah 2000 kisi menjadi grup yang lebih kecil, melakukan perhitungan untuk setiap grup, dan kemudian menggabungkan perhitungan tersebut.
whuber
2
@whuber: 's' adalah objek kecil, hanya sekelompok petunjuk ke file. Fungsi calc, seperti fungsi lain dalam paket raster, tidak akan memuat semua nilai ke dalam memori; itu akan memprosesnya dalam potongan. Artinya, memecah menjadi kelompok seperti yang Anda sarankan dilakukan secara otomatis, di belakang layar. Ukuran chunk dapat dioptimalkan untuk jumlah RAM yang tersedia melalui rasterOptions ().
Robert Hijmans
1
Terima kasih telah menjelaskannya! Saya berasumsi, tanpa memeriksa, itu stackdan calcbekerja seperti kebanyakan Rfungsi lainnya dengan terlebih dahulu memuat semua data ke dalam RAM.
whuber
+1 Menyukai kesimpulan R vs contoh Python yang disediakan ...
SlowLearner
5

Seperti yang saya perhatikan dalam komentar, umumnya Anda harus menghindari penggunaan R untuk tujuan non-statistik karena masalah kinerja dalam aspek-aspek tertentu (bekerja dengan siklus adalah contoh). Berikut adalah contoh kode untuk Anda di Pyhton (terima kasih untuk artikel ini ) untuk klasifikasi ulang satu file dengan satu band. Anda akan dapat memodifikasinya dengan mudah untuk pemrosesan batch jika Anda tahu cara mendapatkan semua file dari direktori . Perhatikan bahwa raster direpresentasikan sebagai array, jadi Anda dapat menggunakan metode array untuk meningkatkan kinerja bila berlaku. Untuk bekerja dengan array di Python lihat dokumentasi Numpy .

UPD: kode yang saya posting pada awalnya adalah versi terpotong dari filter khusus yang diperlukan per pemrosesan piksel. Tetapi untuk pertanyaan ini, penggunaan Numpy akan meningkatkan kinerja (lihat kode).

from osgeo import gdal
import sys
import numpy

gdalData = gdal.Open("path_to_file")
if gdalData is None:
  sys.exit("ERROR: can't open raster")

#print "Driver short name", gdalData.GetDriver().ShortName
#print "Driver long name", gdalData.GetDriver().LongName
#print "Raster size", gdalData.RasterXSize, "x", gdalData.RasterYSize
#print "Number of bands", gdalData.RasterCount
#print "Projection", gdalData.GetProjection()
#print "Geo transform", gdalData.GetGeoTransform()


raster = gdalData.ReadAsArray()
xsize = gdalData.RasterXSize
ysize = gdalData.RasterYSize

#print xsize, 'x', ysize

## per pixel processing
#for col in range(xsize):
#  for row in range(ysize):
#    # if pixel value is 16 - change it to 7
#    if raster[row, col] == 16:
#      raster[row, col] = 7

#reclassify raster values equal 16 to 7 using Numpy
temp = numpy.equal(raster, 16)
numpy.putmask(raster, temp, 7)

# write results to file (but at first check if we are able to write this format)
format = "GTiff"
driver = gdal.GetDriverByName(format)
metadata = driver.GetMetadata()
if metadata.has_key(gdal.DCAP_CREATE) and metadata[gdal.DCAP_CREATE] == "YES":
  pass
else:
  print "Driver %s does not support Create() method." % format
  sys.exit(1)
if metadata.has_key(gdal.DCAP_CREATECOPY) and metadata[gdal.DCAP_CREATECOPY] == "YES":
  pass
else:
  print "Driver %s does not support CreateCopy() method." % format
  sys.exit(1)

# we already have the raster with exact parameters that we need
# so we use CreateCopy() method instead of Create() to save our time
outData = driver.CreateCopy("path_to_new_file", gdalData, 0)
outData.GetRasterBand(1).WriteArray(raster)
SS_Rebelious
sumber
4
"R lambat ketika melakukan siklus (loop)" sering disalahgunakan sebagai alasan untuk menghindari R. Ya, jika Anda memutar sel-sel raster dalam R itu akan lambat, tetapi paket raster bekerja pada seluruh raster sekaligus , dan memiliki banyak kode C dan berjalan dengan kecepatan C. Untuk raster dengan ukuran seperti itu, sebagian besar pekerjaan akan berada pada kecepatan C, overhead looping tidak signifikan.
Spacedman
@ Spacedman, ya 'raster' adalah paket yang berguna (dan saya menyukainya), tetapi saya tidak pernah puas dengan kinerjanya bahkan ketika loop tidak terlibat.
SS_Rebelious
2
Oke, bandingkan waktu yang dibutuhkan dalam R dengan waktu yang dibutuhkan dalam Python. Bisakah Anda tidak beroperasi pada array numpy keseluruhan daripada perulangan?
Spacedman
@Spacedman, saya baru saja memperbarui jawabannya.
SS_Rebelious
Terima kasih banyak untuk kalian berdua. Saya akan mencoba bermain-main dengan kode Python yang Anda berikan dan beberapa R dan melihat apa yang bisa saya capai. Saya akan memperbarui dengan hasil atau masalah.
Nigel
2
  1. Jangan gunakan readGDAL. Bunyinya menjadi objek Spasial * yang mungkin bukan ide yang baik ..

  2. Gunakan rasterpaket. Itu dapat membaca hal-hal GDAL menjadi objek Raster. Ini adalah hal yang baik. r = raster("/path/to/rasterfile.tif")akan membacanya r.

  3. Klasifikasi Anda kemudian t = r > 4 & r <= 9

  4. Pertanyaan besar adalah apakah akan menampilkan ini ke file raster baru dan kemudian melakukan langkah ringkasan di loop lain. Jika Anda memiliki penyimpanan, saya akan menulisnya ke file baru hanya karena jika loop Anda gagal (karena salah satu dari 2000 file tersebut adalah sampah) Anda harus mulai lagi. Jadi gunakan writeRasteruntuk membuat raster yang di-threshold jika Anda memutuskan untuk melakukannya.

  5. Maka loop Anda hanyalah sesuatu seperti

ini.

count = raster(0,size of your raster)
for(i in 1:number of rasters){
  r = raster(path to binary raster file 'i')
  count = count + r
}
return(count)

Manajemen memori R mungkin mengenai Anda di sini - saat Anda melakukannya, count=count+rR mungkin akan membuat salinan baru count. Jadi itu berpotensi 2.000 salinannya - tetapi pengumpulan sampah akan masuk dan membersihkan, dan ini di sini bahwa perulangan R dulu sangat buruk.

Dalam hal pengaturan waktu, pada PC 4yo saya, op ambang membutuhkan sekitar 1,5 detik pada raster sebesar angka acak antara nol dan dua puluh. Times 2000 = Anda mengerjakan matematika. Seperti biasa, buat set tes kecil untuk mengembangkan kode, lalu luncurkan pada data Anda yang sebenarnya dan nikmati kopi.

Spacedman
sumber
Saya ingin tahu tentang apa yang Anda maksud dengan "op ambang." Pada sistem saya (menjalankan R2.15.0), membaca grid 1,6 megapiksel (dalam format floating point ESRI asli) membutuhkan waktu 0,19 detik dan melakukan lima operasi loop - dua perbandingan, konjungsi, tambahan, dan toko - mengambil yang lain 0,09 detik, untuk 0,28 detik per iterasi. Ketika sebaliknya loop dilakukan dengan caching 100 grid sekaligus ke dalam array dan menggunakan rowSumsuntuk melakukan total, penggunaan RAM tentu saja naik: tetapi, sama jelasnya, tidak ada pengumpulan sampah. Waktu turun hanya ke 0,26 detik per iterasi.
whuber