Total area under multiple overlapping polygons

I have shape files from multiple years (eight total years) representing insecticide treatments over an area of land. I can easily find the total area treated, and compare overlap among two individual years. However, I'm stuck in developing a concise code to determine the total area treated two times, three times, for times, etc.

Simplified image attached as an example, each year is a different colour. Some of the area is only treated once (single colour), some is treated twice (overlapping colours) and some is treated 2+ times.

How do I go about finding the total area treated 1, 2, 3, 4, (etc.) times?

I cannot create example code like the map image (or share my shape files). However, here is a very simplified example of area treated per year with some overlap among years (as requested):

library(tidyverse)
library(sf)

polygon1 <- st_polygon(list(
  matrix(c(0, 0, 1, 0, 1, 1, 0, 1, 0, 0), ncol = 2, byrow = TRUE)
))

polygon2 <- st_polygon(list(
  matrix(c(2, 2, 3, 2, 3, 3, 2, 3, 2, 2), ncol = 2, byrow = TRUE)
))

polygon3 <- st_polygon(list(
  matrix(c(0.5, 0.5, 1.5, 0.5, 1.5, 1.3, 0.5, 1.3, 0.5, 0.5), ncol = 2, byrow = TRUE)
))

polygon4 <- st_polygon(list(
  matrix(c(2, 0.5, 2.5, 0.5, 2.5, 2.5, 2, 2.5, 2, 0.5), ncol = 2, byrow = TRUE)
))

polygon5 <- st_polygon(list(
  matrix(c(1.2, 1.5, 1.8, 1.5, 1.8, 1.8, 1.2, 1.8, 1.2, 1.5), ncol = 2, byrow = TRUE)
))

polygon6 <- st_polygon(list(
  matrix(c(0.8, 0.8, 2.2, 0.8, 2.2, 1.2, 0.8, 1.2, 0.8, 0.8), ncol = 2, byrow = TRUE)
))

polygon7 <- st_polygon(list(
  matrix(c(1.3, 2.2, 2.1, 2.2, 2.1, 2.8, 1.3, 2.8, 1.3, 2.2), ncol = 2, byrow = TRUE)
))


year1 <- st_multipolygon(list(polygon1, polygon2)) |> 
  st_sfc()

year2 <- st_multipolygon(list(polygon3, polygon4, polygon5))|> 
  st_sfc()

year3 <- st_multipolygon(list(polygon6, polygon7))|> 
  st_sfc()


ggplot() +
  geom_sf(data = year1, fill = "#009E73", color = 'black', alpha = 0.5) +
  geom_sf(data = year2, fill = "#E69F00", color = 'black', alpha = 0.5) +
  geom_sf(data = year3, fill = "#0072B2", color = 'black', alpha = 0.5) +
  theme_classic()

How would could I determine the total area covered three times, the total area covered only twice and the total area covered only once.

Solution

using the raster approach and your example data:

gather the various simple features into one spatial dataframe (could be done manually of course):

library(sf)


    the_polygons <- 
      do.call(rbind,
              lapply(list('year1', 'year2', 'year3'),
                  FUN = \(o_name) data.frame(year = o_name,
                                             geometry = get(o_name)
                                             )
                  )
      ) |>
      st_sf()

## > the_polygons
## Simple feature collection with 3 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 0 xmax: 3 ymax: 3
## CRS:           NA
##    year                       geometry
## 1 year1 MULTIPOLYGON (((0 0, 1 0, 1...
## 2 year2 MULTIPOLYGON (((0.5 0.5, 1....
## 3 year3 MULTIPOLYGON (((0.8 0.8, 2....

create a base raster with the overall extent of the polygons (obtained with st_bbox and desired resolution. This empty base raster is later populated with the polygon values.

    library(terra) ## the workhorse for raster, as {sf} is for features

    raster_base <- the_polygons |> st_bbox() |> ext() |> rast(resolution = .1)

create the output raster by rasterizeing the polygons into the base raster, suming the count of non-empty values (one per polygon covering the particular cell):


    raster_output <- rasterize(the_polygons, raster_base, fun = sum)

extract and tabulate the cell values:


    table(values(raster_output))

##   1   2   3 
## 325  78   7

(multiply cell count by cell size to obtain area at your resolution)

visual control:


    plot(raster_output)