8

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))
plot(polys.df2, col="black", add = T)

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

0

2 Answers 2

11

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

diffPoly <- st_difference(polys.df_sf, st_union(polys.df_sf2))
2
  • Nice! Is this technique useful in many cases we have overlapping geometries, perhaps?
    – Kazuhito
    Commented 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
    Commented Sep 15, 2021 at 17:17
5

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

diffPoly <- rmapshaper::ms_erase(polys.df_sf, polys.df_sf2)
3
  • I'll give it a go. The limitation on data size within R might be a bit of a pain.
    – Simon
    Commented 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
    Commented Sep 15, 2021 at 17:16
  • @Matt, thanks for testing. Accepted this answer
    – Simon
    Commented Sep 15, 2021 at 18:06

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.