3

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)

basic polygons

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)

number of overlaps

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

basic difference

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)

progress so far

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

3
  • 1
    You should probably spell out the problem here, in the event that any linked question is deleted then this question becomes useless. Give example data if possible, and some example of the code you've tried.
    – Spacedman
    Commented Apr 19, 2020 at 11:41
  • 1
    It is an honor to receive attention from @Spacedman. I have added a reproducible example.
    – ess
    Commented Apr 19, 2020 at 14:01
  • I think the voronoi polygon part necessarily has to be within the original polygon (in this case the squares), so st_within(partition, b) should tell you that, and then you match that to the other bits and join them. My intuition about the "within" requirement might be wrong though, and no time to experiment at the moment.... But nice question, if people will undo their close votes...
    – Spacedman
    Commented Apr 19, 2020 at 17:57

2 Answers 2

4

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)

answer


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

0

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
}
 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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