Splitting overlap between polygons and assign to nearest polygon using R?

Question

My question is how to perform the process described in the answer to this similar question in R, ideally using the sf package. In the linked case this is done using using PostGIS with some sf commands, but I am looking for a native R solution.

Basically: given a set of overlapping (multi)polygons, how do I split the overlapping polygons according to the non-overlapping sections which are closest.

Here is a reproducible example along with my progress so far (with help from the r-spatial vignettes and github).

library(sf)
library(tidyverse)

pol = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
b = st_sfc(pol, pol + c(.8, .2), pol + c(.2, .8))
par(mar = rep(0, 4))
plot(b, col = NA)

i = st_intersection(st_sf(b))
par(mar = rep(0, 4))
cl = sf.colors(3, categorical = TRUE)
plot(st_geometry(b))
plot(st_geometry(i[i$n.overlaps == 3,2]), col = cl[1], add = TRUE)
plot(st_geometry(i[i$n.overlaps == 2,2]), col = cl[2], add = TRUE)

d = st_difference(b)
plot(d, col = cl)

So far I have been able to create the divisions I want at the voronoi partitions among the centroids of each polygon:

independent <- b %>% st_sf %>% st_intersection %>% subset(n.overlaps<=1)
overlap <- b %>% st_sf %>% st_intersection %>% subset(n.overlaps>1) %>% st_union()
partition <- b %>% st_centroid %>% st_union %>% st_voronoi %>% st_cast %>% st_intersection(overlap)

plot(st_geometry(independent), col=cl)
plot(st_geometry(partition), col=cl, add=TRUE)

How do I assign each partition of the overlap to the appropriate polygon from the non-overlapping sections? Is this the best method?

Answer 1

Update:

In addition to the original answer below, I was able to adapt a similar question from the sf github for this purpose:

# credit to https://github.com/r-spatial/sf/issues/824
library(sf)
library(tidyverse)

st_no_overlap <- function(polygons) {

  centroids <- polygons %>% st_centroid

     # Voronoi tesselation
     voronoi <- 
          centroids %>% 
          st_geometry() %>%
          st_union() %>%
          st_voronoi() %>%
          st_collection_extract()

     # Put them back in their original order
     voronoi <-
          voronoi[unlist(st_intersects(centroids,voronoi))]

     # Keep the attributes
     result <- centroids

     # Intersect voronoi zones with buffer zones
     st_geometry(result) <-
          mapply(function(x,y) st_intersection(x,y),
                 #st_buffer(st_geometry(centroids),dist), 
                 polygons$geometry,
                 voronoi,
                 SIMPLIFY=FALSE) %>%
          st_sfc(crs=st_crs(centroids))

     result
}

pol = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
b = st_sfc(pol, pol + c(.8, .2), pol + c(.2, .8))

plot(st_no_overlap(st_sf(geometry=b)), col=cl)

---

With help from @Spacedman and this additional question, here is the answer I came up with:

library(sf)
library(tidyverse)

pol = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
b = st_sfc(pol, pol + c(.8, .2), pol + c(.2, .8))

independent <- b %>% st_sf %>% st_intersection %>% subset(n.overlaps<=1)
overlap <- b %>% st_sf %>% st_intersection %>% subset(n.overlaps>1) %>% st_union()
partition <- b %>% st_centroid %>% st_union %>% st_voronoi %>% st_cast

merge_list <- st_within(partition %>% st_intersection(overlap), b)

merged_list <- lapply(1:length(merge_list), function(i){st_sf(st_intersection(partition[i], b[merge_list[[i]]]))})

new_b <- do.call(rbind, merged_list)
plot(new_b, col=cl)

Answer 2

For reasons I haven't fully dug into, the function in the answer by ess marked correct wasn't working for the inputs I had. Namely, at the intersection step, I received the error:

Error in geos_op2_geom("intersection", x, y, ...) : 
  st_crs(x) == st_crs(y) is not TRUE

Further issues appeared where the number of pieces created by the intersection was far higher than my original input and the st_sfc step at the end also complained. I'm using sf 1.0-7.

I made a few modifications here that still solve the original prompt while avoiding some of the issues described above:

st_no_overlap <- function(polygons) {
  
  centroids <- polygons %>% st_centroid
  
  # Voronoi tesselation
  voronoi <- 
    centroids %>% 
    st_geometry() %>%
    st_union() %>%
    st_voronoi() %>%
    st_collection_extract() # now it's sfc
  
  # Put them back in their original order
  voronoi <-
    voronoi[unlist(st_intersects(centroids,voronoi))]
  
  # Keep the attributes
  result <- centroids
  
  st_geometry(result) <-
    mapply(
      function(x,y) {
        z <- 
          st_intersection(
            x,
            y
          ) %>% 
          # this can create multiple parts, so we union.
          st_union()
      },
      # we need this to produce a list to iterate over
      # in parallel with voronoi elements, so we 
      # convert to sfc
      st_as_sfc(polygons[attributes(polygons)$sf_column]),
      voronoi,
      SIMPLIFY=FALSE
    ) %>% 
    # st_sfc() returned errors, but st_as_sfc() did not
    st_as_sfc(crs = st_crs(centroids))
  
  result
}