R sf: Find polygons outside of multiple overlapping polygons

Question

I have a large shapefile that contains 1000s of multiple, overlapping polygons. I am trying to the combined area outside of these multiple, overlapping polygons that do not overlap with any polygon in the layer. These are effectively regions where fires have burned multiple times, I am looking for where fires have only burned once.

I am stuck on how to find the area outside that does not have any overlaps. Can you help?

Here is a reproducible example.

#make some overlapping polygons
m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))
p = st_polygon(list(m))
n = 100
l = vector("list", n)
for (i in 1:n)
  l[[i]] = p + 2 * runif(2)
s = st_sfc(l)

#select just a few of these
s5 <- s[1:5]

#now try to step through and get the non-overlapping areas
counter <- 0
sall.out <- list()
for (i in 1:length(s5)) {
    print(i)
    s5.sel <- s5[i]
    s5.out <- s5[!(s5 == s5.sel)] #select all the polygons outside
    s5.int <- st_intersection(s5.out, s5.sel) #intersect the outside polygons with the one selected

    #step through and find all differences between the selected region and the intersected
    for (j in 1:length(s5.int)) {
        print(j)
        s5.out <- st_difference(s5.sel, s5.int[j])
    
        counter <- counter+1
        sall.out[[counter]] <- s5.out
    }
}

plot(s5)
plot(s5.sel, add=T, col="red")
plot(s5.int, add=T, col="blue")
plot(s5.out, add=T, col="pink")

So, now I have all of the sall.out in a list, but how can I remove the ones that overlap each other and flatten the list?

Thank you.

Answer 1

I suggest that you use the convenient property of st_intersection. From the documentation:

When called with missing y, the sfc method for st_intersection returns all non-empty intersections of the geometries of x; an attribute idx contains a list-column with the indexes of contributing geometries.

This basically "segments" the plane, and returns one geometry per segment. When you convert your polygons to sf instead of sfc, this also means you get n.overlaps and origins columns that describe where each geometry came from in the original input. You can then simply filter and see that the overlapping areas have been removed.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 7.1.0
set.seed(1)
m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))
p = st_polygon(list(m))
n = 100
l = vector("list", n)
for (i in 1:n)
  l[[i]] = p + 2 * runif(2)


s = st_sf(polygon = 1:n, geometry = st_sfc(l))
s5 <- s[1:5, ]
plot(s5["polygon"])

non_overlaps <- st_intersection(s5) %>%
  filter(n.overlaps == 1)

plot(non_overlaps["polygon"])

^{Created on 2020-07-21 by the reprex package (v0.3.0)}