Bagaimana cara melakukan klip GIS sebenarnya dari lapisan poligon menggunakan lapisan poligon di R?

Saya ingin melakukan GIS Clip yang benar dalam R poligon tanah menggunakan serangkaian poligon batas tunggal, tetapi saya tidak dapat menemukan fungsi R untuk melakukannya dengan benar. Seharusnya berfungsi seperti clipfungsi di ESRI's ArcMap. Saya sudah mencoba overmetode ini dalam sppaket tetapi tampaknya tidak berfungsi untuk polys daripada polys.

Salah satu saran adalah menggunakan paket gIntersectionin rgeossebagai klip menggunakan kode berikut:

#------------------------------------
library(rgeos)
library(maptools)

#Read layers as SpatialPolygonsDataFrame (both the same Albers projection)
Soils_poly = readShapePoly("Soils_polygons")  #Note - Has 400 polygons
clipper_poly = readShapePoly("clipper_polygon")  #Note - Has 1 polygon

#Try gintersection as clip 
Clipped_polys = gIntersection(Clipper_Tile_poly, Soils_poly)

#-----------------------------------

Ini membutuhkan waktu 5 menit untuk berjalan (terlalu lambat) dan kesalahan dengan ini:

Kesalahan di RGEOSBinTopoFunc (spgeom1, spgeom2, byid, id, drop_not_poly, "rgeos_intersection"): TopologyException: tidak ada dirEdge keluar yang ditemukan di -721459.77681285271 2009506.5980877089

Saya juga mencoba kode ini untuk memeriksa tumpang tindih:

gIntersects(Clipper_Tile_poly, Soils_poly)

dan hasilnya BENAR. clipfungsi dalam ESRI ArcMap berfungsi dengan baik untuk data ini.

Adakah yang tahu fungsi R untuk melakukan klip video spasial dengan benar menggunakan spasial?

Petunjuk yang diberikan oleh @mdsummer tentang penggunaan byid=TRUEkarya secara akurat.

Lihat contoh yang dapat direproduksi, di bawah ini:

library(rgeos)
library(sp)

#Create SpatialPlygons objects
polygon1 <- readWKT("POLYGON((-190 -50, -200 -10, -110 20, -190 -50))") #polygon 1
polygon2 <- readWKT("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20))") #polygon 2

par(mfrow = c(1,2)) #in separate windows
plot(polygon1, main = "Polygon1") #window 1
plot(polygon2, main = "Polygon2") #window 2

polygon_set <- readWKT(paste("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20),",
                     "(-190 -50, -200 -10, -110 20, -190 -50))"))

par(mfrow = c(1,1)) #now, simultaneously
plot(polygon_set, main = "Polygon1 & Polygon2")

clip <- gIntersection(polygon1, polygon2, byid = TRUE, drop_lower_td = TRUE) #clip polygon 2 with polygon 1
plot(clip, col = "lightblue")

GT <- GridTopology(c(-175, -85), c(10, 10), c(36, 18))
gr <- as(as(SpatialGrid(GT), "SpatialPixels"), "SpatialPolygons")
plot(gr)

clip2 <- gIntersection(clip, gr, byid = TRUE, drop_lower_td = TRUE)
plot(clip2, col = "pink")

Anda juga dapat menggunakan paket raster raster::intersect(spdf1, spdf2). Ini memiliki keuntungan mempertahankan atribut jika Anda memiliki SpatialPolygonsDataFrame.

library(sp); library(rgeos)

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")))
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")))

# works but drop the attributes
#gIntersection(spdf1, spdf2, byid=T)

#better to keep attributes
inter1=raster::intersect(spdf1, spdf2)

plot(spdf1, col="red")
plot(spdf2, col="yellow", add=T)
plot(inter1,col="blue", add=T)

Berkat pertanyaan ini untuk menunjukkan itu dan untuk kode sampel.