Comment réaliser un vrai clip SIG de couche de polygones en utilisant une couche de polygones dans R?

Je voudrais faire un vrai clip SIG dans R de polygones de sols en utilisant une série de polygones à frontière unique, mais je ne trouve pas de fonction R pour le faire correctement. Cela devrait fonctionner exactement comme la clipfonction d'ArcMap d'ESRI. J'ai essayé la overméthode dans le sppaquet mais cela ne semble pas fonctionner pour les polys sur polys.

Une suggestion a été d'utiliser le package gIntersectionin rgeoscomme un clip en utilisant le code suivant:

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

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

Cela prend 5 minutes à exécuter (beaucoup trop lent) et des erreurs avec ceci:

Erreur dans RGEOSBinTopoFunc (spgeom1, spgeom2, byid, id, drop_not_poly, "rgeos_intersection"): TopologyException: aucun dirEdge sortant trouvé à -721459.77681285271 2009506.5980877089

J'ai également essayé ce code pour vérifier le chevauchement:

gIntersects(Clipper_Tile_poly, Soils_poly)

et le résultat était VRAI. clipdans ESRI ArcMap fonctionne très bien pour ces données.

Quelqu'un connaît une fonction R pour faire correctement un clip sur des polygones spatiaux en utilisant des polygones spatiaux?

L'astuce fournie par @mdsummer d'utiliser les byid=TRUEœuvres avec précision.

Voir l'exemple reproductible ci-dessous:

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

Vous pouvez également utiliser le package raster raster::intersect(spdf1, spdf2). Il a l'avantage de conserver les attributs au cas où vous auriez un 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)

Merci à cette question de le signaler et pour un exemple de code.