Klip poligon dan simpan data?

13

Saya memiliki dua poligon ini:

library(sp); library(rgeos); library(maptools)

coords1 <- matrix(c(-1.841960, -1.823464, -1.838623, -1.841960, 55.663696,
                    55.659178, 55.650841, 55.663696), ncol=2)
coords2 <- matrix(c(-1.822606, -1.816790, -1.832712, -1.822606, 55.657887,
                    55.646806, 55.650679, 55.657887), ncol=2)
p1 <- Polygon(coords1)
p2 <- Polygon(coords2)
p1 <- Polygons(list(p1), ID = "p1")
p2 <- Polygons(list(p2), ID = "p2")
myPolys <- SpatialPolygons(list(p1, p2))
spdf1 = SpatialPolygonsDataFrame(myPolys, data.frame(variable1 = c(232,
                                                                   242), variable2 = c(235, 464), row.names = c("p1", "p2")))
proj4string(spdf1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf1, col="red")

coords1a <- matrix(c(-1.830219, -1.833753, -1.821154, -1.830219, 55.647353,
                     55.656629, 55.652122, 55.647353), ncol=2)
p1a <- Polygon(coords1a)
p1a <- Polygons(list(p1a), ID = "p1a")
myPolys1 <- SpatialPolygons(list(p1a))
spdf2 = SpatialPolygonsDataFrame(myPolys1, data.frame(variable1 = c(2),
                                                      variable2 = c(3), row.names = c("p1a")))
proj4string(spdf2) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf2, col="yellow", add=T)

masukkan deskripsi gambar di sini

Saya ingin memotong bagian spdf1yang berpotongan dengan spdf2. Namun, saya ingin spdf1tetap sebagai SpatialPolygonsDataFrame dan menyimpan informasi apa pun yang ada di dalamnya spdf1@data.

Saya sudah mencoba gDifference sebagai berikut, yang memotong bagian-bagian spdf1yang berpotongan dengan spdf2, tetapi kemudian dikonversi spdf1ke SpatialPolygons (yaitu membuang informasi yang terkandung dalam spdf1@data).

gDifference(spdf1, spdf2, byid=T)

Bagaimana saya bisa dipotong menjadi spdf1dengan spdf2tetapi mempertahankan data yang terdapat dalam spdf1@data? Saya telah memeriksa dan mencoba pertanyaan serupa ini tanpa cara overlay poligon pada SpatialPointsDataFrame dan menyimpan data SPDF?

luciano
sumber

Jawaban:

9

Pendekatan yang paling sederhana adalah

  library(raster)
  x <- spdf1 - spdf2

  # or, more formally
  y <- erase(spdf1,  spdf2)

Lihat? 'Paket raster' (bagian XIV) untuk lebih banyak fungsi yang berhubungan dengan overlay poligon. Fungsi-fungsi ini menggunakan fungsi-dasar rgeo di bawah tenda, di fungsi 'tingkat pengguna' (sebagai lawan dari 'tingkat pengembang').

Robert Hijmans
sumber
Apa yang Anda maksud dengan fungsi "di tingkat pengguna '(sebagai lawan dari' tingkat pengembang ')"?
luciano
rgeosmenyediakan operasi geometris tetapi tidak berurusan dengan atribut data. Jadi menggunakan fungsi-fungsi ini membutuhkan banyak pekerjaan untuk menjaga semuanya tetap bersama. Fungsi raster menyederhanakan ini dan berperilaku seperti fungsi serupa dalam perangkat lunak GIS,
Robert Hijmans
Ya, tetapi ini dapat menyebabkan Kesalahan dalam SpatialPolygonsDataFrame (part2, x @ data [match (row.names (part2),: row.names data dan ID Poligon tidak cocok
jebyrnes
Itu akan menjadi bug. Bisakah Anda memberikan contohnya?
Robert Hijmans
4

Solusinya adalah menambahkan kembali atribut setelah melakukan klip, saat mengkonversi dari SpatialPolygonske SpatialPolygonsDataFrame.

sp3 <- gDifference(spdf1, spdf2, byid = TRUE)
row.names(sp3) <- row.names(spdf1)

spdf3 <- SpatialPolygonsDataFrame(sp3, data = spdf1@data)

spdf3@data

   variable1 variable2
p1       232       235
p2       242       464

plot(spdf3, col="red")

masukkan deskripsi gambar di sini

Andre Silva
sumber
Ini terlihat seperti masalah yang saya miliki, hanya klip output dalam contoh khusus saya menghasilkan rownames dari spdf1 yang tidak ada (sebagai gsub sederhana untuk menghilangkan digit ke-2 berturut-turut. Nama harus mengurus pertandingan rownames, bukan? )
jebyrnes
Kesalahan dalam SpatialPolygonsDataFrame (klip, data = as.data.frame (all_spdfs_together @ data)): Ketidakcocokan panjang objek: klip memiliki 1718 objek Poligon, tetapi as.data.frame (all_spdfs_together @ data) memiliki 86 baris
jebyrnes
Tentu - Saya punya banyak poligon di lautan. Beberapa ditempatkan secara tidak benar di tanah, atau tumpang tindih dengan tanah. Saya ingin memotong itu jadi saya hanya memiliki daerah yang ada di lautan. Saya memiliki shapefile garis pantai untuk dibandingkan. Inilah kodenya - gist.github.com/jebyrnes/c2e8d2b6c82849dad3a813d952ab8bb0
jebyrnes
1
nevermind - the raster :: erase solution berfungsi (tidak pernah sebelumnya karena alasan yang aneh)
jebyrnes