R – Clipping Buffers to Features in R Using sf Package

rsf

I'm trying to create an sf object of polygons which are buffer regions clipped to the polygons used to create the buffers. Using st_difference returns an sf object with all combinations of both sf objects. Example below:

library(sf)
library(tidyverse)

p1<-st_polygon(list(rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))))
p2<-st_polygon(list(rbind(c(3,1),c(3,2),c(4,2),c(4,1),c(3,1))))

p1p2<-st_sfc(p1,p2)

d<-data.frame(name=c("Poly1","Poly2"))

st_geometry(d)<-p1p2

dBuff<-st_buffer(d,dist=1.5)

dBuffClip<-st_difference(dBuff,d)

dBuffClip

I want only these ggplot()+geom_sf(data=dBuffClip[1,],fill="lightblue"):

enter image description here

And none of these ggplot()+geom_sf(data=dBuffClip[3,],fill="lightblue"):
enter image description here

I've realized the I can filter dBuffClip using filter(name==name.1) but this is very inefficient and crashes my machine because of the size of the data. Does somebody know of an efficient way to do this?

Best Answer

If you explode the geometry so you have POLYGON and not MULTIPOLYGON (using sf::st_cast) then, you should have a one to one match between feature and buffer. Now, you can apply an iterator (eg., for, lapply) to subset the feature and corresponding buffer and apply sf::st_difference. Results can be combined using rbind. This may take awhile with a big dataset but, will keep things memory safe.

Here, lets add the sf package and read in some internal data, which is projected into albers (so we have a real distance). We will then grab 10 random polygons.

library(sf)
  d = st_transform(st_read(system.file("shape/nc.shp", 
                   package="sf"), quiet=TRUE),  
                   st_crs(5070))
    d <- d[sample(1:nrow(d), 10),]

We then buffer our polygons to 5000m and cast them to POLYGON to ensure that we do not end up with MULTIPOLYGON geometry. This may not always be necessary but, makes things safe.

dBuff <- st_cast(st_buffer(d, dist=5000), "POLYGON")  
  plot(st_geometry(dBuff), col="red")
    plot(st_geometry(d), col="white", add=TRUE)

Now we can use an iterator to go through each feature and coresponding buffer to create polygon with the source feature erased. For the iterator we use lapply and to combine the results we wrap it in do.call. The suppressWarnings call is to simply suppress an irrelevant warning message.

dBuffClip <- do.call("rbind", lapply(1:nrow(d), function(i) {
  suppressWarnings(st_difference(dBuff[i,], d[i,]))}))
    plot(st_geometry(dBuffClip))

With your example

d <- st_as_sf(st_sfc(st_polygon(list(rbind(c(1,1),c(1,2), c(2,2),
              c(2,1),c(1,1)))), st_polygon(list(rbind(c(3,1),c(3,2),
              c(4,2),c(4,1),c(3,1))))))

dBuff <- st_cast(st_buffer(d, 1.5), "POLYGON")
  dBuffClip <- do.call("rbind", lapply(1:nrow(d), function(i) {
        suppressWarnings(st_difference(dBuff[i,], d[i,]))}))

plot(st_geometry(dBuffClip), col=c("green","red"))

example result

Related Question