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

Question

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

#1 de Andre Silva (19 votos)
#2 de jeb (1 votos)

15

Me gustaría hacer un verdadero Clip GIS en R de polígonos de suelos utilizando una serie de polígonos de un solo contorno, pero no puedo encontrar una función R para hacerlo correctamente. Debería funcionar igual que la función CLIP en ArcMap de ESRI. He probado el método Over en el paquete SP, pero no parece funcionar para polys sobre polys. Una sugerencia fue utilizar gIntersection in rgeos package como Clip usando 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)

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

Esto tarda 5 minutos en ejecutarse (demasiado lento) y errores con esto:

Error en RGEOSBinTopoFunc (spgeom1, spgeom2, byid, id, drop_not_poly, "rgeos_intersection"): TopologyException: no se ha encontrado un dirEdge saliente en -721459.77681285271 2009506.5980877089

También probé este código para verificar la superposición:

gIntersects(Clipper_Tile_poly, Soils_poly)

y el resultado fue VERDADERO. La función CLIP en ESRI ArcMap funciona bien para estos datos.

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

polygon r layers clip

pregunta user29199 14.04.2014 - 23:52

2 respuestas

Lea otras preguntas en las etiquetas polygon r layers clip

¿Cómo dividir un polígono en varios polígonos de igual ancho en ArcMap 10? Diseño de base de datos: topología de carreteras en SIG

score 19 · Answer 1

La sugerencia proporcionada por @mdsummer sobre el uso de byid=TRUE funciona con precisión.

Vea el ejemplo reproducible, abajo:

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)#clippolygon2withpolygon1plot(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")

score 1 · Answer 2

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

Graciasa esta pregunta por señalarlo y por el código de muestra.