4

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
8

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
    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 at 17:17
2

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
    Sep 2 '19 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 at 17:16
  • @Matt, thanks for testing. Accepted this answer
    – Simon
    Sep 15 at 18:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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