Is there an efficient way to calculate all pairwise intersection of polygons in sf (without using for loops) in R?

Question

I want to compute the pairwise intersection of polygons from a variable in a data frame.

Here is what I'm looking for :

set.seed(131)
library(sf)
m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))
p = st_polygon(list(m))
n = 5
l = vector("list", n)
for (i in 1:n)
  l[[i]] = p + 2 * runif(2)
s = st_sfc(l)
s.f = st_sf(s)
s.f$id = c(1,1,2,2,3)
plot(s.f.2, col = s.f.2$id) 
    
s.f.2 = s.f %>% group_by(id) %>%   summarise(geometry = sf::st_union(s)) # %>% summarise(area = st_area(s))

s.f.2$area = st_area(s.f.2)


all.ids.pol = unique(s.f.2$id)
df = data.frame(NULL)
for (i in 1:length(all.ids.pol)) {
  for (j in 1:length(all.ids.pol)) {
    ol1 = st_intersection(s.f.2[c(all.ids.pol[i],all.ids.pol[j]),])
    if (dim(ol1)[1]>2 | i == j) {
      #add in an area count column (area of each arable poly in intersect layer)
      ol1$areaover <- st_area(ol1$geometry)
    } else {ol1$areaover = 0}
    #run the intersect function, converting the output to a tibble in the process
    if (i == j) {
      ol1.1 <- as_tibble(ol1)[1,]  
    } else {ol1.1 <- as_tibble(ol1)[2,]  }
    
    id.names = c(all.ids.pol[i],all.ids.pol[j])
    my.df =data.frame(area = ol1.1$areaover,id1 = id.names[1],id2 =id.names[2])
    df = rbind(df,my.df)
  }    
}
intersected.areas = df %>% 
  arrange(id1) %>%
  mutate(area.ha = units::set_units(area, ha)) %>% 
  select(-area) %>% 
  pivot_wider(names_from = id1, values_from = area.ha) %>%  
  arrange(id2) %>%
  rename(ID = id2)

In the end I want a matrix that shows in rows and columns the "ID" and each cell in the matrix showing the area of intersection between the polygons.

# A tibble: 3 × 4
     ID   `1`   `2`   `3`
  <dbl>  [ha]  [ha]  [ha]
1     1 1.59  0.501     0
2     2 0.501 1.86      0
3     3 0     0         1

So the code works, but I'm wondering if there exists a function in sf that can do that in a more efficient way than using for loop (this becomes very unpractical for large datasets)

Answer 1

library(sf)
library(tidyverse)
s.f.2 %>%
  st_intersection(s.f.2) %>%
  mutate(intersect_area = st_area(.)) %>%
  st_drop_geometry() %>%
  pivot_wider(id_cols = id, names_from = id.1, 
              values_from = intersect_area,
              values_fill = 0)

# # A tibble: 3 x 4
#      id   `1`   `2`   `3`
#   <dbl> <dbl> <dbl> <dbl>
# 1     1 1.59  0.501     0
# 2     2 0.501 1.86      0
# 3     3 0     0         1