R: Unduh DEM besar, ubah proyeksi, dan sesuaikan ke skala yang lebih kecil

11

Ini adalah proses yang hanya membutuhkan beberapa detik dalam perangkat lunak GIS. Upaya saya untuk melakukannya di R menggunakan sejumlah besar memori kemudian gagal. Apakah ada yang salah dalam kode saya, atau ini hanya sesuatu yang R tidak bisa lakukan? Saya telah membaca R dapat bekerja di dalam Rumput, dapatkah saya menggunakan fungsi Rumput dari dalam R?

library(raster)

# I have many environmental rasters in this format
new_r <- raster(ncol=615, nrow=626, xmn=-156.2, xmx=-154.8, ymn=18.89, ymx=20.30)
res(new_r) <- 0.00225
projection(new_r) <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0"

R> new_r ### not too big with a few hundred cells per side
class       : RasterLayer 
dimensions  : 627, 622, 1  (nrow, ncol, nlayers)
ncell       : 389994 
resolution  : 0.00225, 0.00225  (x, y)
projection  : +proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0 
extent      : -156.2, -154.8, 18.89, 20.3  (xmin, xmax, ymin, ymax)
values      : none

# I get the DEM at much higher resolution (zipfile is 182Mb)
zipurl <- "ftp://soest.hawaii.edu/coastal/webftp/Hawaii/dem/Hawaii_DEM.zip"
DEMzip <- download.file(zipurl, destfile = "DEMzip")
unzip("DEMzip", exdir = "HIDEM")
HIDEM <- raster("HIDEM/hawaii_dem")

R> HIDEM ### 10m resolution, file is way too big
class       : RasterLayer 
dimensions  : 15067, 13136, 1  (nrow, ncol, nlayers)
ncell       : 197920112 
resolution  : 10, 10  (x, y)
projection  : +proj=utm +zone=5 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
extent      : 179066, 310426, 2093087, 2243757  (xmin, xmax, ymin, ymax)
values      : HIDEM/hawaii_dem 
min value   : 0 
max value   : 4200 

# the following line fails (after a long time)
new_HIDEM <- projectRaster(HIDEM, new_r)
J. Win.
sumber
Hanya ingin tahu, paket apa yang Anda gunakan?
djq
@celenius: paket ini disebutraster
J. Win.

Jawaban:

9

Dari pandangan saya pada sumbernya, rastersepertinya menebak apakah dataset cocok dengan memori, dan jika demikian, lakukan operasi dalam memori, jika tidak pada disk. Anda dapat memaksanya untuk melakukan perhitungan dengan menetapkan secara eksplisit chunksize(sel untuk diproses pada satu waktu) dan maxmemory(jumlah sel maksimum untuk dibaca ke dalam memori):

setOptions(chunksize = 1e+04, maxmemory = 1e+06)

Atau, Anda dapat melakukan transformasi dengan GDAL secara langsung:

gdalwarp -t_srs '+proj=utm +zone=5 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0' HIDEM/hawaii_dem hawaii_dem_utm.tif

Ini kemungkinan akan menjadi opsi tercepat, dan tidak memerlukan pengaturan lingkungan GIS secara eksplisit.

scw
sumber
Itu tidak berhasil, tetapi ini berhasil: setOptions(chunksize = 1e+04, maxmemory = 1e+06)Waktu delapan menit, jauh lebih sedikit daripada yang diperlukan untuk menginstal dan menggunakan GIS yang sebenarnya.
J. Win.
@J. Winchester: Saya telah memperbarui respons saya untuk memasukkan pengaturan Anda karena itulah pendekatan yang lebih baik. Penulis paket, kemungkinan akan tertarik mendengar kapan dan mengapa crash, dan mudah-mudahan memperbarui default untuk mencerminkan ini.
scw
1
itu ide yang baik untuk menambahkan kompresi (lossless) dan ubin (default ke 256x256) ke GeoTIFF dari gdalwarp jika target Anda dapat mengatasinya: -co COMPRESS = LZW -co TILED = YA
mdsumner
Saya memberi tahu Robert Hijmans tentang kasus ini. Pada DEM yang lebih kecil, pengaturan default hampir-optimal, jadi ini adalah misteri sejauh ini.
J. Win.
Bagus! Ini juga memungkinkan saya untuk mengekspor RasterLayer ke (3GB) netcdf dengan writeRaster.
David LeBauer
3

Anda juga dapat menggunakan paket spgrass6 untuk integrasi antara R dan rumput. Penulisnya adalah Roger Bivand (penulis sp)

Paket ini memiliki banyak fungsi untuk melengkapi menjalankan rumput di dalam R (atau sebaliknya) dan bertukar data antara R dan rumput

untuk informasi lebih lanjut: http://cran.r-project.org/web/packages/spgrass6/index.html

dickoa
sumber
1

HAI,

Ini adalah proses yang hanya membutuhkan beberapa detik dalam perangkat lunak GIS. Upaya saya untuk melakukannya di R menggunakan> jumlah memori yang besar kemudian gagal.

Anda menjawab pertanyaan Anda, lakukan itu dalam GRASS atau GDAL dan tinggalkan R untuk tugas-tugas lain.

Pablo
sumber
1
Untuk kelengkapan: Anda harus melihat paket spgrass untuk menjalankan rumput di R.
johanvdw
1
Dan opsi ketiga adalah menggunakan saga gis. Ada modul (RSAGA) yang menghubungkan saga dan R.
johanvdw
Fungsi R ini dirancang untuk menggunakan GDAL, tetapi tampaknya tidak menggunakannya dengan baik dalam kasus ini. Pertanyaan saya adalah "Bagaimana saya bisa menyelesaikan tugas ini dengan R", bukan "Perangkat lunak SIG apa yang tersedia yang dapat melakukan tugas ini."
J. Win.