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?
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...