Non overlapping buffers/morphological tessellation for polygons in sf

Question

I am trying to replicate the tessellation functionality from the momepy package (see source code here, lines 61 to 550) in R, using sf. This is basically a kind of mixture between a buffer and a voronoi polygon in which the size of the buffer around each unit (buildings) is determined by a limit but also by the buffers around other units. It's very clear here.

I've found a similar question asked here and it gets very close to where I want it to go, but not really. In concrete, I'm using the following code. The data I'm using can be found here.

``` r
library(sf)
library(tmap)
library(tidyverse)


# mostly from here https://github.com/ropensci/stplanr/issues/362

p <- st_read("export.geojson")


p <- p %>%
  select(id, building, name, geometry) %>%
  filter(grepl("MULTIPOLYGON", st_geometry_type(geometry)))

p <- st_transform(p, "EPSG:25831")


centroids <- st_centroid(p)



voronoi_polygons_c <- st_voronoi(x = do.call("c", st_geometry(centroids))) %>%
  st_collection_extract() %>%
  st_set_crs(st_crs(p))


voronoi_polygons_for_buildings_c <- voronoi_polygons_c %>% 
  st_set_crs(st_crs(p)) %>% 
  # now I need to merge the polygons associated with the same line
  # https://github.com/r-spatial/sf/issues/1030
  st_cast("MULTIPOLYGON", ids = unlist(st_intersects(voronoi_polygons_c, centroids))) %>% 
  st_union(by_feature = TRUE)

example_buffer <- st_buffer(p, dist = 20)


my_solution_c <- map2(
  .x = st_geometry(example_buffer), 
  .y = voronoi_polygons_for_buildings_c, 
  .f = st_intersection
) %>% 
  st_sfc() %>%
  st_set_crs("EPSG:25831")


# this is the result

tm_shape(p) + 
  tm_polygons(col = "blue",
              border.col = "red",
              lwd = 2) + 
  tm_shape(my_solution_c) + 
  tm_borders(col = "white")

This is how it looks like. The thing is, as you can see, that in the case of large buildings, the voronoi polygon around the centroid is smaller than the building itself.

I have tried another possibility, which is casting the polygons to multipoint objects and then creating the voronoi polygons with them, as follows (this is what the github page I've referenced suggests doing).

mps <- p %>%
  st_cast(to = "POINT") %>%
  group_by(id) %>%
  summarise()

voronoi_polygons <- st_voronoi(x = do.call("c", st_geometry(mps))) %>%
  st_collection_extract() %>%
  st_set_crs(st_crs(p))

voronoi_polygons_for_buildings <- voronoi_polygons %>% 
  st_set_crs(st_crs(p)) %>% 
  # now I need to merge the polygons associated with the same line
  # https://github.com/r-spatial/sf/issues/1030
  st_cast("MULTIPOLYGON", ids = unlist(st_intersects(voronoi_polygons, mps))) %>% 
  st_union(by_feature = TRUE)

example_buffer <- st_buffer(p, dist = 20)


my_solution <- map2(
  .x = st_geometry(example_buffer), 
  .y = voronoi_polygons_for_buildings, 
  .f = st_intersection
) %>% 
  st_sfc() %>%
  st_set_crs("EPSG:25831")

# this is the result

tm_shape(p) + 
  tm_polygons(col = "blue",
              border.col = "red",
              lwd = 2) + 
  tm_shape(my_solution) + 
  tm_borders(col = "white")

But the result is even worse as it makes a polygon around each vertex of each building, thus.

Answer 1

Found a solution using the {raster} and {stars} packages. It is a bit slow but it returns the result I'm looking for. Mostly taken from here

https://stackoverflow.com/questions/65571860/how-to-calculate-the-distance-of-raster-to-the-nearest-polygon

The idea is to create a raster stack with the distance from every pixel to every polygon and get the minimum. Then, transform to a multipolygon object, project, and join the building ids. Any suggestions on how to improve it are very welcome.

library(sf)
library(tmap)
library(raster)
library(stars)
library(tidyverse)

p <- st_read("export.geojson")

p <- p %>%
  select(id, building, name, geometry) %>%
  filter(grepl("MULTIPOLYGON", st_geometry_type(geometry)))

qtm(p)

p <- st_transform(p, "EPSG:25831")

# Raster method - from here https://stackoverflow.com/questions/65571860/how-to-calculate-the-distance-of-raster-to-the-nearest-polygon


get_tessellations <- function(x){
  bbox <- st_bbox(x)
  r <- raster(xmn =  bbox[[1]], xmx = bbox[[3]], ymn = bbox[[2]], ymx = bbox[[4]], res=0.5) # taken from raster
  r[] <- sample(1:10,ncell(r), replace=T)
  st <- stack()
  p$id <- as.factor(x$id) # transform id to factor - you need to create an unique id column beforehand
  for(i in x$id){
    r_pol <- rasterize(x[x$id==i,], r, field="id")
    rd <- distance(r_pol)
    st <- stack(st,rd)
    
  }
  r_min_which_pol <- which.min(st) # nearest polygon of each cell
  r_min_which_pol <- stars::st_as_stars(r_min_which_pol)
  
  tessellation <- st_as_sf(r_min_which_pol, as.points = FALSE, merge = TRUE) %>%
    group_by(layer) %>%
    summarise() # create one multipolygon object for each layer
  st_crs(tessellation) <- st_crs(x)
  x$layer <- seq.int(1,nrow(x)) # this is to make the merge - the layer object will be equivalent to the row number of original polygon
  tessellation <- tessellation %>%
    left_join(select(as_tibble(p), id, layer), by = "layer") # assign unique id of polygons to each tessellation
  return(tessellation)
}

tes <- get_tessellations(p)

tm_shape(p) + 
  tm_polygons(col = "blue",
              border.col = "red",
              lwd = 2) + 
  tm_shape(tes) + 
  tm_borders(col = "blue")

The result being as follows.