Find intersection of multiple polygons by group

user_5137706 0
Viewed 464

So finding the union of a more than two polygons turns out to be straight forward using st_union and group_by() from dplyr. Finding the intersection with st_intersection turns out to not be as simple since st_intersection does only work for finding intersection of two objects. What's the best way to find the intersection of multiple polygons by group? Here's some example data

s1 <- rbind(c(1, 1), c(10, 1), c(10, 10), c(1, 10), c(1, 1))
s2 <- s1 + 4
s3 <- s1 - 4

group <- c(rep(1, 3), rep(2, 3))

df <- data.frame("group" = c(rep(1, 3), rep(2, 3)))

df <- data.frame(cbind(df, rep((c(st_sfc(st_polygon(list(s1))), st_sfc(st_polygon(list(s2))), st_sfc(st_polygon(list(s3))))), 2)))

What I want is to replace the geometry column with the intersection in each group, that is the object created by


p <- list(s1 = s1, s2 = s2, s3 = s3)

p  <-  lapply(p, function(x) st_sfc(st_polygon(list(x))))

intersection <- accumulate(p, st_intersection)$s3

intersection <- st_sfc(st_polygon(intersection))

So the final data would look like


df <- data.frame("group" = c(rep(1, 3), rep(2, 3)))

df <- data.frame(st_sf(cbind(df, rep(st_sfc(st_polygon(intersection)), 6))))
sf dplyr intersection r
2 Answers

Managed to do it with a loop but would prefer a dplyr solution

agg <- as.data.frame(matrix(NA, 0, 2))
for (i in 1:nrow(df)) {
  group <- df$group[i]
  vec <- df$geometry[df$group == group]
  my.list <- as.list(vec)
  output <- Reduce(st_intersection, my.list) %>% st_geometry()
  my.df <- data.frame("group" = group)
  res <- as.data.frame(st_sf(cbind(my.df, output)))
  
  agg <- rbind(agg, res)
}
user_51377060

First up, st_intersection() can work with only a single R object, for example if you supply it with an sf collection it will return it's self-intersection; check out the sf page here for more details, and another answer as an example here. The resulting collection includes features with attributes including n.overlaps (number of layers) and origins (which polygons contribute to each overlapping area).

So for this example, we employ st_intersection() and keep only the geometries where n.overlaps == 3 (as there are 3 geometries in each group). Basically the steps are: make df into an sf collection, then we use split() to split the sf collection by group, find the self-intersection of each group, join both groups back together, then keep only the areas with 3 overlaps.

# example data from question ~~~~~

library(sf)
library(tidyverse)

s1 <- rbind(c(1, 1), c(10, 1), c(10, 10), c(1, 10), c(1, 1))
s2 <- s1 + 4
s3 <- s1 - 4

group <- c(rep(1, 3), rep(2, 3))
df <- data.frame("group" = c(rep(1, 3), rep(2, 3)))
df <- data.frame(cbind(df, rep((c(st_sfc(st_polygon(list(s1))), st_sfc(st_polygon(list(s2))), st_sfc(st_polygon(list(s3))))), 2)))


# solution ~~~~~

df %>%
  st_as_sf %>% 
  split(., .$group) %>% 
  map(st_intersection) %>% 
  bind_rows %>% 
  filter(
    n.overlaps == 3
  )

# Simple feature collection with 2 features and 3 fields
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 5 ymin: 5 xmax: 6 ymax: 6
# CRS:           NA
#     group n.overlaps origins                       geometry
# 1.3     1          3 1, 2, 3 POLYGON ((5 5, 5 6, 6 6, 6 ...
# 4.3     2          3 1, 2, 3 POLYGON ((5 5, 5 6, 6 6, 6 ...

Unfortunately, this method isn't entirely dplyr as it uses split() and map() to iterate through each group, however it is pipe-able. I thought group_by() might have worked in place of split() %>% map() %>% bind_rows(), but unfortunately I couldn't get st_intersection() to respect group_by().

Hopefully this is still useful, and perhaps someone can find a way to swap out split() and map() for group_by().

user_134787490
Related