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:


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

  • Try gIntersection with byid=TRUE (I think it's a problem with dangling topologies for certain clips, and sometimes doing it this way helps), for speed check out gUnarySTRtreeQuery() or gBinarySTRtreeQuery() to identify intersecting bounding boxes of pairs of polygons and only intersect those pairs. There's no easy high-level wrappers for all of this afaik
    – mdsumner
    Apr 14, 2014 at 22:37

3 Answers 3


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

See the reproducible example, below:


#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

polygons side-by-side

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here


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

  • 1
    Works a dream - unbelievably fast. I'm using this to clip polylines. I wanted to make subsets of the shapefile for British rivers, because it's much faster to work with a smaller subset of data.
    – CJB
    Apr 6, 2016 at 9:12
  • 1
    @AndreSilva, thought I had, but guess I hadn't! Re: community, I wish the R GIS developer community was a bit larger TBH. I will probably have to fall back on ArcMap for some digitization, and other processes that are quick in the GUI. Sep 9, 2017 at 3:05
  • 1
    raster::intersect is faster and preserves the @data see gis.stackexchange.com/questions/173257/… (edit, or see @jeb below)
    – geneorama
    Mar 23, 2022 at 14:50

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)

enter image description here

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

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

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

### 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
writeOGR(luxembourg_clip, paste0(getwd(), "/Output_shapes"), "lux_clip", driver="ESRI Shapefile")

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

enter image description here

enter image description here

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.