R sf: st_intersection on a list of polygons

Question

I have a list of polygons (and multi-polygons) named p_1, p_2, ..., p_n. And I would like to obtain the area in which they all intersect. As st_intersection() does not accept lists as arguments, I tried the following three approaches. None of them provides a satisfactory solution, which is why I am looking for alternative, more efficient techniques.

(i) I could loop through the list

for(i in P) p_1 <- st_intersection(p_1, i)

where P is a list containing polygons p_2 to p_n. But that is rather slow.

(ii) A do.call() approach, i.e.

p <- do.call(st_intersection, P)

where P is a list of polygons p_1 to p_n, only computes the intersection between the first two polygons in the list.

(iii) I could combine the polygons into one sf object and then run st_intersection():

p <- do.call(c, P) %>% 
   st_sf() %>% 
   st_intersection()

It works but is slow. Presumeably because it also derives a lot of other polygons apart from the common intersection of all polygons in P.

None of the three approaches provides a satisfactory solution. Looping through a hierarchy of pairwise comparisons in a parallelized framework might be faster. However, I assume there to be a simpler and more efficient solution than that.

Any comments and suggestions are welcome.

A note to the person who closed this question yesterday: do not close this question. Comment or send me a private message, if you personally have a problem with it. But do not close it.

user_12500315 · Answer

I don't think the overhead of iterating through a list is a problem here: finding the intersection of multiple polygons is just computationally expensive. However, the method of sequentially applying a function to members of a list (effectively what you were trying to do with do.call) is easily managed using purrr::accumulate:

You don't have a reproducible example for folks here to test possible solutions, and creating sf polygons from scratch involves some work, so that may have been why your previous question was closed - I don't know.

Anyway, lets create three overlapping squares in a list and draw them:

library(sf)
library(purrr)

# create square
s1 <- rbind(c(1, 1), c(10, 1), c(10, 10), c(1, 10), c(1, 1))
p  <- list(s1 = s1, s2 = s1 + 4, s3 = s1 - 4)
p  <-  lapply(p, function(x) st_sfc(st_polygon(list(x))) )

plot(p[[1]], xlim = c(-5, 15), ylim = c(-5, 15))
plot(p[[2]], add = TRUE)
plot(p[[3]], add = TRUE)

Our goal is to find the intersection of all three squares, which of course is the tiny square in the center. Using purrr, this is as easy as:

intersection <- accumulate(p, st_intersection)$s3

So when we add our result, coloured red, we get:

plot(intersection, col = "red", add = TRUE)

In terms of performance, accumulate is only about 10% faster than the raw loop, so you may need to parallelize this if performance is a big problem. Also, if there is a possibility that there is no intersection between all the polygons, you can find your smallest polygon and use st_intersects to ensure that all the polygons actually intersect it. This is a much quicker calculation provided that there is a fair chance of there being no unqiue intersection.

R sf: st_intersection on a list of polygons

1 Answers