# Using simple features (sf) in R, how do I erase polygons overlapping with another layer

With a regular spatialDataFrame I could use `erase(x, y)` or `x-y`. How can I perform the equivalent using sf objects?

As an example:

``````##Create layers based on this answer: http://stackoverflow.com/questions/26620373/spatialpolygons-creating-a-set-of-polygons-in-r-from-coordinates
set.seed(145)
square1 <- t(replicate(50, {
o <- runif(2)
c(o, o + c(0, 0.1), o + 0.1, o + c(0.1, 0), o)
}))
ID <- paste0('sq', seq_len(nrow(square1)))

polys1 <- SpatialPolygons(mapply(function(poly, id) {
xy <- matrix(poly, ncol=2, byrow=TRUE)
Polygons(list(Polygon(xy)), ID=id)
}, split(square1, row(square1)), ID))

polys.df <- SpatialPolygonsDataFrame(polys1, data.frame(id=ID, row.names=ID))

square2 <- t(replicate(5, {
o <- runif(2)
c(o, o + c(0, 0.1), o + 0.1, o + c(0.1, 0), o)
}))
ID <- paste0('sq', seq_len(nrow(square2)))

polys2 <- SpatialPolygons(mapply(function(poly, id) {
xy <- matrix(poly, ncol=2, byrow=TRUE)
Polygons(list(Polygon(xy)), ID=id)
}, split(square2, row(square2)), ID))
polys.df2 <- SpatialPolygonsDataFrame(polys2, data.frame(id=ID, row.names=ID))

#plot the layers
plot(polys.df, col=rainbow(50, alpha=0.5))

#Compare methods - erase and difference
erasedPoly <- erase(polys.df, polys.df2)
plot(erasedPoly, col=rainbow(50, alpha=0.5))

polys.df_sf <- st_as_sf(polys.df)
polys.df_sf2 <- st_as_sf(polys.df2)
diffPoly <- st_difference(polys.df_sf, polys.df_sf2)
plot(diffPoly, col=rainbow(50, alpha=0.5), max.plot = 1)
``````

Erase here does what I want, but st_difference doesn't. My workaround, if there is no in-built function in sf, will be to convert back to a Spatial PolygonDataFrame.

The solution was to union the geometries first, like so:

``````diffPoly <- st_difference(polys.df_sf, st_union(polys.df_sf2))
``````
• Nice! Is this technique useful in many cases we have overlapping geometries, perhaps? May 17, 2017 at 11:44
• A note of caution: when using this solution on a dataset of projected polygons, this produced artifacts and null geometries.
– Matt
Sep 15, 2021 at 17:17

`rmapshaper` library is super helpful. It's also faster.

``````diffPoly <- rmapshaper::ms_erase(polys.df_sf, polys.df_sf2)
``````
• I'll give it a go. The limitation on data size within R might be a bit of a pain. Sep 2, 2019 at 6:23
• I tried both solutions, and this one produced clean polygons without geometric errors. and @plnnr is right, it's fast. As much as it's annoying to add yet another library to the workflow, this solution seems to be superior.
– Matt
Sep 15, 2021 at 17:16
• @Matt, thanks for testing. Accepted this answer Sep 15, 2021 at 18:06