I have these two polygons:

library(sp); library(rgeos); library(maptools)

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")))
proj4string(spdf1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
plot(spdf1, col="red")

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")))
proj4string(spdf2) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
plot(spdf2, col="yellow", add=T)

enter image description here

I want to cut out parts of spdf1 that are intersected by spdf2. However, I want spdf1 to remain as a SpatialPolygonsDataFrame and to retain any information contained within spdf1@data.

I¹ve tried gDifference as follows, which cuts out parts of spdf1 that are intersected by spdf2, but then converts spdf1 to SpatialPolygons (i.e. discarding information contained in spdf1@data).

gDifference(spdf1, spdf2, byid=T)

How can I cut into spdf1 with spdf2 but retain data contained in spdf1@data? I have checked and tried this similar question without how to overlay a polygon over SpatialPointsDataFrame and preserving the SPDF data?


2 Answers 2


The simplest approach would be

  x <- spdf1 - spdf2

  # or, more formally
  y <- erase(spdf1,  spdf2)

See ?'raster-package' (section XIV) for more functions that deal with polygon overlay. These functions use the base-functions of rgeos under the hood, in 'user-level' (as opposed to 'developer-level') functions.

  • What do you mean by "in 'user-level' (as opposed to 'developer-level') functions"?
    – luciano
    Commented Nov 9, 2015 at 21:16
  • rgeos provides geometric operations but does not deal with the attributes of the data. Thus using these functions requires a lot of work to keep everything together. The raster functions simplify this and the behave like similar functions in GIS software, Commented Nov 9, 2015 at 21:27
  • Yeah, but this can lead to Error in SpatialPolygonsDataFrame(part2, x@data[match(row.names(part2), : row.names of data and Polygons IDs do not match
    – jebyrnes
    Commented Jan 12, 2018 at 16:51
  • That would be a bug. Can you provide an example of that? Commented Jan 13, 2018 at 18:31

A workaround would be to re-add the attributes after doing the clip, while converting from SpatialPolygons to SpatialPolygonsDataFrame.

sp3 <- gDifference(spdf1, spdf2, byid = TRUE)
row.names(sp3) <- row.names(spdf1)

spdf3 <- SpatialPolygonsDataFrame(sp3, data = spdf1@data)


   variable1 variable2
p1       232       235
p2       242       464

plot(spdf3, col="red")

enter image description here

  • This looks like a problem I have, only the output clip in my particular instance produces rownames from spdf1 that do not exist (as a simple gsub to get rid of the 2nd digit in row.names should take care of the rownames match, no?)
    – jebyrnes
    Commented Jan 12, 2018 at 17:11
  • Error in SpatialPolygonsDataFrame(clip, data = as.data.frame(all_spdfs_together@data)) : Object length mismatch: clip has 1718 Polygons objects, but as.data.frame(all_spdfs_together@data) has 86 rows
    – jebyrnes
    Commented Jan 12, 2018 at 17:18
  • Sure - I have a bunch of polygons of things in the ocean. Some are incorrectly placed on land, or overlap with land. I want to clip that out so I only have areas that are in the ocean. I have a coastline shapefile to compare against. Here's the code - gist.github.com/jebyrnes/c2e8d2b6c82849dad3a813d952ab8bb0
    – jebyrnes
    Commented Jan 12, 2018 at 20:12
  • 1
    nevermind - the raster::erase solution works (it hadn't before for some odd reason)
    – jebyrnes
    Commented Jan 12, 2018 at 20:32

Your Answer

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

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