1

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?

  • 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 Apr 19 at 11:41
  • 1
    It is an honor to receive attention from @Spacedman. I have added a reproducible example. – ess Apr 19 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 Apr 19 at 17:57
2

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

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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