TransWikia.com

How to create multiring buffer with sf in R

Geographic Information Systems Asked by Bernd V. on May 28, 2021

I have a polygon of an area which I like to buffer multiple times with different distances for intersection with another layer buffer(ring)-wise later (this part already works!)

What I do not achieve, is to create a single sf object from the list of buffers (to be able to export it as a layer).
The intersections from another layer with those buffers (done in the same loop) DOES get combined to a sf object by do.call(rbind ,…)!

I’m lost in the complexity of those object types.

Example:

#Polygon
x <- structure(list(AREA = 1489841.35074, geometry = structure(list(
    structure(list(structure(c(4371395.39166986, 4371395.39166986, 
    4372634.07252284, 4372634.07252284, 4371395.39166986, 5478810.68590502, 
    5480013.45040013, 5480013.45040013, 5478810.68590502, 5478810.68590502
    ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
    epsg = NA_integer_, proj4string = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam +units=m +no_defs"), .Names = c("epsg", 
"proj4string"), class = "crs"), class = c("sfc_POLYGON", "sfc"
), precision = 0, bbox = structure(c(4371395.39166986, 5478810.68590502, 
4372634.07252284, 5480013.45040013), .Names = c("xmin", "ymin", 
"xmax", "ymax"), class = "bbox"))), .Names = c("AREA", "geometry"
), row.names = c(NA, -1L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(NA_integer_, class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = "AREA"))

library(sf)

#distances for buffers                                                                                                                                                                                                                                                                                                                                                                               "aggregate", "identity"), .Names = "AREA"))                                                                                                                                                                                                                                                                                                                        "aggregate", "identity"), .Names = "AREA")
dist <- c(50,500,1000,3000)

#list holding the buffer rings
lstbuffers <-list()
for (i in seq(dist)){
  #first buffer
  if (i == 1){
    lstbuffers[[i]]<- st_buffer(x,dist[i])

  }else{ 
  #next buffers as rings around the previous
    lstbuffers[[i]] <- st_difference(st_buffer(x,dist[i]),st_buffer(x, dist[i-1]))
  }
  lstbuffers[[i]]$radius <-dist[i] #add the radius
}  

buffers <- do.call(rbind,lstbuffers) #DOES NOT WORK, I get a matrix of lists :(

So, I can do operations with my lstbuffers[i] list items, but fail to combine them into one layer.

3 Answers

When you do an st_difference(a,b) the result has all the attributes from both features. So if a has a $NAME and b has a $LABEL then the difference will have both a $NAME and a $LABEL attribute. If the features have attributes of the same name, eg AREA, then R will create a new name for the second one with a number - so you get $AREA and $AREA.1.

Your first feature in the list comes from a buffer operation on a single feature, so it only has the $AREA attribute. Further features come from a difference operation and so have $AREA and $AREA.1 features, and so can't be rbind-ed.

Solution: remove $AREA from x before doing anything:

`x$AREA=NULL`

then do all the buffering and differencing and the rbind should work.

Correct answer by Spacedman on May 28, 2021

I created a function to return a multi-buffered polygon layer from a sf object in 3 steps. There's plenty room for improvement (many for loops, etc..) but works for the purpose of producing buffers of equal distance.

# list of buffers polygons at different lenght

st_multibuffer <- function(data=x, from, to, by, nQuadSegs = 30) {
# get a multi buffered polygon. Require an sf object
  library(sf)
  library(dplyr)

  seq.buffer <-seq(from,to, by) # allocate a  vector of length equal to number of buffer required

  if(inherits(x,"sf")) {

    # create a list of buffers 
    k <- vector('list', length(seq.buffer))
    for (i in 1:length(seq.buffer)) {
      k[[i]] <- st_buffer(x =x, dist = seq.buffer[i], nQuadSegs = nQuadSegs)
      k[[i]]$idx <- i
      k[[i]]$distance <- seq.buffer[i]
      st_agr(k[[i]]) <- "constant"
    }

    # clip from the bigger to the lower 

    l <- vector('list', length(seq.buffer))

    for (i in length(seq.buffer):1) {

      if(i==1){
        l[[i]] <- st_difference(k[[i]], st_geometry(x))
      } 
      else {
        l[[i]] <- st_difference(k[[i]], st_geometry(k[[i - 1]]))
      }
    }
    # Join all multipolygon

    if (length(seq.buffer) == 2) {
      temp <- rbind(l[[1]], l[[2]])
    } else {
      temp <- rbind(l[[1]], l[[2]])
      for (m in 3:length(seq.buffer)) {
       temp <- rbind(temp, l[[m]])
      }
    }

    return(temp)

  } else {

    stop("x is not a sf object")

  }

}

Answered by user1591727 on May 28, 2021

Here is a function that ask for three parameters:

  • x --> sf_layer
  • n --> number of rings
  • d --> distance between rings (m)

Moreover this add and ID to each ring based on the distance transformed to km.

#-------------------------------------------------------------------------------
library(sf); library(mapview)

# Create a ringed buffer with a number of rings and a distance between them
Multiring <- function(x,n,d){
  buffers <- list(); names <- list(); nd <- d
  for (i in 1:n){
    buffers[[i]] <- st_as_sf(st_union(st_buffer(x,nd)))
    buffers[[i]]$ID <- paste0("Buffer ", round(nd/1000,1), " km")
    nd <- nd+d
  }
  
  jlayers <- function(x){ 
    if (length(x)==1){ 
      # if length is == 1 , return 1 layer 
      xm <- x[[1]] 
    } else { 
      for (i in 1:(length(x)-1)){ 
        if(i==1){xm <- x[[1]]} 
        xm <- rbind(xm, x[[i+1]]) 
      } 
    } 
    return(xm) 
  }
  
  return(jlayers(buffers))
}

mapview(layer, 3, 5000)

enter image description here

Answered by César Arquero on May 28, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP