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.

Simon
  • 835
  • 1
  • 7
  • 13

2 Answers2

11

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

diffPoly <- st_difference(polys.df_sf, st_union(polys.df_sf2))
Simon
  • 835
  • 1
  • 7
  • 13
  • Nice! Is this technique useful in many cases we have overlapping geometries, perhaps? – Kazuhito May 17 '17 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 '21 at 17:17
5

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

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