¿Cómo realizar un verdadero recorte SIG de capa de polígonos utilizando una capa de polígonos en R?

Question

Me gustaría hacer un verdadero clip SIG en R de polígonos de suelos utilizando una serie de polígonos de límites simples, pero no puedo encontrar una función R para hacerlo correctamente. Debería funcionar igual que el clip en ArcMap de ESRI. He probado la función over método en sp pero no parece funcionar para polis sobre polis.

Una sugerencia fue utilizar el gIntersection en rgeos como un clip utilizando el siguiente código:

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

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

Tarda 5 minutos en ejecutarse (demasiado lento) y da error con esto:

Error en RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_not_poly, "rgeos_intersection") : TopologyException: no outgoing dirEdge found at -721459.77681285271 2009506.5980877089

También he probado este código para comprobar si hay solapamiento:

gIntersects(Clipper_Tile_poly, Soils_poly)

y el resultado fue VERDADERO. clip de ESRI ArcMap funciona correctamente para estos datos.

¿Alguien conoce una función de R para hacer correctamente un recorte en polígonos espaciales utilizando polígonos espaciales?

Answer 1

La sugerencia de @mdsummer de utilizar byid=TRUE funciona con precisión.

Véase el ejemplo reproducible a continuación:

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

---

Editar :

Si se desea realizar la intersección manteniendo los datos de atributos, una opción es la función raster::intersect. Véase respuesta de jeb .

Answer 2

También puede utilizar el paquete raster raster::intersect(spdf1, spdf2) . Tiene la ventaja de mantener los atributos en caso de tener 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)

Gracias a esta pregunta para señalarlo y para el código de ejemplo.

Answer 3

library(raster)
# library(rgdal) # if require
# library(rgeos) # if require

### load a shape from raster package
luxembourg <- shapefile(system.file("external/lux.shp", package="raster"))
plot(luxembourg) 

### draw a simple shape
another_shape <- as(extent(6, 6.4, 49.75, 50), 'SpatialPolygons')
plot(another_shape, add=T, border="red", lwd=3)

### crop (SpatialPolygon) luxembourg with another_shape
luxembourg_clip <- crop(luxembourg, another_shape) 

### See the outputs
plot(luxembourg_clip, add=T, col=4) ### plot on 
plot(luxembourg_clip, add=F, col=3) ### just plot

### Write on disk
library(rgdal)
writeOGR(luxembourg_clip, paste0(getwd(), "/Output_shapes"), "lux_clip", driver="ESRI Shapefile")

Fuente: https://rdrr.io/cran/raster/man/crop.html