How to perform a true GIS clip of polygons layer using a polygon layer in R?

Question

I would like to do a true GIS Clip in R of soils polygons using a series of single boundary polygons, but I cannot find an R function to properly do it. It should work just like the clip function in ESRI's ArcMap. I've tried the over method in sp package but it doesn't seem to work for polys over polys.

One suggestion was to use the gIntersection in rgeos package as a clip using the following code:

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

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

This takes 5 minutes to run (way too slow) and errors with this:

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

I also tried this code to check for overlap:

gIntersects(Clipper_Tile_poly, Soils_poly)

and the result was TRUE. clip function in ESRI ArcMap works fine for this data.

Anyone knows of an R function to properly do a clip on spatial polygons using spatial polygons?

Answer 1

The hint provided by @mdsummer of using byid=TRUE works accurately.

See the reproducible example, below:

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

---

Edit:

If one wants to intersect while keeping the attribute data, then one option is the raster::intersect function. See jeb's answer.

Answer 2

You can also use the raster package raster::intersect(spdf1, spdf2) . It has the advantage of maintaining the attributes in case you have a 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)

Thanks to this question to pointing that out and for sample code.

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

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