Ekstrak Raster dari Raster menggunakan Polygon shapefile di R

13

Saya baru mengenal R dan menggunakan paket raster. Saya memiliki masalah dalam mengekstraksi poligon dari file raster yang ada. Jika saya gunakan

extract(raster, poly_shape)

fungsi pada raster itu selalu membuat daftar dengan data. Yang saya inginkan adalah mengekstrak file raster lain yang dapat saya muat dengan ArcGIS lagi. Setelah membaca sedikit lebih lanjut saya pikir fungsi krop adalah yang benar-benar saya butuhkan. Tetapi ketika saya mencoba menggunakan fungsi ini

crop(raster, poly_shape)

Saya mendapatkan Galat ini:

Error in .local(x, y, ...) : extents do not overlap
In addition: Warning message:
In intersect(extent(x), extent(y)) : Objects do not overlap

File raster dan poly_shape sama untuk kedua fungsi. Bisakah Anda memberi tahu saya apa yang salah di sini? Apakah benar bahwa fungsi memotong membuat raster lain dan bukan daftar?

EDIT : Fungsi Sejauh () tidak bekerja untuk saya. Saya tetap mengalami masalah yang sama. Tapi saya yakin 2 dataset tumpang tindih! Dengan

extract(raster, poly_shape)

Saya mendapatkan data yang benar dari itu. Sama seperti daftar dan bukan sebagai raster seperti saya ingin memilikinya. Saya baru saja memuat dataset di ArcGIS sebelumnya dan mereka sangat cocok jadi saya tidak memeriksa proyeksi. Sekarang saya mencoba

projection(raster) # "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
projection(poly_shape) # "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"

dan Anda dapat melihat bahwa proyeksi tidak sesuai. Fungsi ekstrak tampaknya dapat secara otomatis mengubah file dengan cara yang benar. Saya tahu itu karena saya melakukan yang berikut:

  • Saya memotong bagian yang pasti dari poligon yang saya ekstrak di R juga di ArcGIS
  • Saya menghitung jumlah semua nilai poligon R yang diekstraksi (daftar)
  • Saya menghitung jumlah semua sel raster yang saya potong di ArcGIS

2 memiliki hasil yang sama persis, jadi saya kira kesimpulannya adalah bahwa fungsi ekstrak tidak berfungsi dengan benar. Sekarang saya punya 2 pilihan, saya kira:

  1. Saya perlu cara untuk mengeluarkan Raster dari daftar yang diekstraksi lagi atau
  2. 2 set data (raster + poly_shape) harus menggunakan prjection yang sama dan fungsi crop seharusnya bekerja

Apa yang akan Anda lakukan di sini?

Lars
sumber
Bagaimana jika itu adalah raster 4 band rgbi? Band hilang sejauh ini ...
Doris
Selamat datang di GIS SE! Sebagai pengguna baru, pastikan untuk mengikuti tur singkat . Kemudian pertimbangkan untuk mengedit jawaban Anda untuk memberikan informasi dan referensi tambahan. Lihat Bagaimana cara menulis jawaban yang baik? untuk info lebih lanjut.
Andy
Jika Anda memiliki pertanyaan baru, silakan tanyakan dengan mengklik tombol Ajukan Pertanyaan . Sertakan tautan ke pertanyaan ini jika itu membantu menyediakan konteks. - Dari Ulasan
Erik

Jawaban:

19

Fungsi ekstrak berperilaku tepat seperti seharusnya. Anda bisa memaksa fungsi krop untuk menggunakan tingkat poligon dan kemudian menutupi objek untuk mengembalikan raster tepat yang mewakili area poligon. Jika Anda terus menerima kesalahan itu berarti bahwa data Anda, pada kenyataannya, tidak tumpang tindih.

Harap diingat bahwa R tidak melakukan proyeksi "on the fly" jadi, periksa proyeksi Anda. Anda dapat memeriksa apakah ekstensi Anda tumpang tindih menggunakan fungsi "Sejauh ()".

Berikut adalah contoh tanam menggunakan luas poligon kemudian menutupi raster yang dihasilkan menggunakan poligon "raster".

# Add required packages
require(raster)
require(rgeos)
require(sp)

# Create some data using meuse 
data(meuse)
  coordinates(meuse) <- ~x+y
    proj4string(meuse) <- CRS("+init=epsg:28992")
data(meuse.grid)
  coordinates(meuse.grid) = ~x+y
    proj4string(meuse.grid) <- CRS("+init=epsg:28992")
      gridded(meuse.grid) = TRUE    
        r <- raster(meuse.grid) 
          r[] <- runif(ncell(r))

# Create a polygon
f <- gBuffer(meuse[10,], byid=FALSE, id=NULL, width=250, 
                         joinStyle="ROUND", quadsegs=10)   

# Plot full raster and polygon                       
plot(r)
  plot(f,add=T)

# Crop using extent, rasterize polygon and finally, create poly-raster
#          **** This is the code that you are after ****  
cr <- crop(r, extent(f), snap="out")                    
  fr <- rasterize(f, cr)   
    lr <- mask(x=cr, mask=fr)

# Plot results
plot(lr)
  plot(f,add=T)
Jeffrey Evans
sumber
4
ekstrak () tidak tampil pada proyeksi ulang fly, tapi tanaman () tidak. Itu mungkin menyebabkan beberapa kebingungan
mdsumner
Pemotongan @Jefferey () dan mask () hanya memotong raster sesuai dengan luasan segi empat dari poligon, ia tidak memotongnya dari dalam batas poligon. Adakah yang tahu perintah apa yang dapat memotong raster di dalam batas poligon yang diberikan?
csheth
1
@ Chintan Sheth, untuk mask untuk subset dalam poligon Anda harus memiliki raster yang mewakili nilai-nilai dalam poligon. Inilah sebabnya mengapa Anda merasterisasi subset poligon dan kemudian menutupinya. Langkah pemangkasan adalah mengurangi tingkat raster hingga sama dengan poligon himpunan bagian, yang di masa lalu, jika dilewati, akan menimbulkan kesalahan tingkat ketidakcocokan.
Jeffrey Evans
spTransformdari sppaket (yang kadang-kadang secara otomatis dimuat dengan paket R spasial lainnya) memungkinkan untuk memproyeksi ulang sehingga kedua file dalam proyeksi yang sama misalnya. good_poly=spTransform(spolygon, CRSobj=crs(raster_file))
user3386170
@ user3386170, Hah? Tidak yakin apa yang Anda maksud. Pertanyaan ini muncul pada saat paket raster baru saja menambahkan "on the fly projection" dalam beberapa fungsi. Fungsi-fungsi ini sebelumnya telah melemparkan kesalahan ketika proyeksi tidak cocok, tetapi posting ini dari 2014. Anda juga harus menyadari selalu memuat rgdal saat menggunakan sp :: spTransform () karena, ia menambahkan fungsionalitas tambahan, penting, di belakang layar.
Jeffrey Evans
8

Apa yang sebenarnya saya cari adalah mask()fungsinya.

mask(raster, poly_shape)

berfungsi tanpa kesalahan dan mengembalikan apa yang saya cari.

Lars
sumber
2
Proyeksi ulang data Anda sehingga berada di ruang proyeksi yang sama. Bahkan di ArcGIS, di mana on the fly proyeksi otomatis, praktik yang sangat buruk untuk melakukan analisis dalam proyeksi yang berbeda. Dengan data dalam proyeksi yang berbeda, data Anda tidak akan berbagi luasan umum, yang merupakan kesalahan yang Anda terima.
Jeffrey Evans
Gunakan projectExtent () untuk mendapatkan batas untuk memotong raster.
mdsumner
Agar sesuai dengan format T&J situs ini harus ditempatkan ke badan utama pertanyaan sebagai edit / pembaruan (dan kemudian tambahkan komentar ke jawaban ini dalam "balasan" agar mereka tahu ada lebih banyak untuk dilihat).
matt wilkie
@mattwilkie Maaf karena tidak pas dengan format tetapi teks saya terlalu panjang untuk dikirim sebagai komentar di sini. @JeffreyEvans saya mencoba berikut: projection(raster) = projection(poly_shape)dan sebaliknya projection(poly_shape) = projection(raster)tetapi kedua cara menghasilkan kesalahan yang sama: Error in .local(x, y, ...) : extents do not overlap In addition: Warning message: In intersect(extent(x), extent(y)) : Objects do not overlap. Apakah ada cara bagaimana saya bisa melihat proyeksi mana yang digunakan dengan cepat menggunakan fungsi ekstrak () (karena itu jelas bekerja)?
Lars
1
Apa yang sebenarnya saya cari adalah fungsi mask (). mask(raster, poly_shape)berfungsi tanpa kesalahan dan mengembalikan apa yang saya cari.
Lars
-1

Extent berfungsi dengan baik ... Saya pikir Xmin, Xmax, Ymin, dan Ymax sejauh Anda berbeda dari X dan Y raster Anda - yaitu mereka diatur berlawanan

ToNoY
sumber