Общая площадь под несколькими перекрывающимися полигонами

используя растровый подход и данные вашего примера:

• собрать различные простые функции в один пространственный фрейм данных (конечно, это можно сделать вручную):

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....

• создайте базовый растр с общей протяженностью полигонов (полученной с помощью st_bbox и желаемым разрешением. Этот пустой базовый растр позже заполняется значениями полигонов.

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

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

• создайте выходной растр, rasterizeвключив полигоны в базовый растр, sumподсчитав количество непустых значений (по одному на полигон, охватывающий конкретную ячейку):

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

• извлеките и сведите в таблицу значения ячеек:

    table(values(raster_output))

##   1   2   3 
## 325  78   7

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

визуальный контроль:

    plot(raster_output)

Для векторного метода я бы подошел к этому, объединив данные за все годы в один слой с переменной, идентифицирующей год, а затем выполнил самопересечение. В результате будет создан индекс, который можно будет использовать для обратной ссылки на объединенный годичный слой, чтобы получить подробную информацию о каждом перекрытии, которую вы можете объединить в новое поле, в котором перечислены годы. Затем это можно объединить, суммируя полученный слой научной фантастики по интересующим вас переменным. В случае действительно больших слоев этап пересечения может занять некоторое время.

#assign a value to each year

year1 <- st_multipolygon(list(polygon1, polygon2)) %>%
  st_sfc() %>%
    st_union() %>%
      st_as_sf() %>%
        mutate(Year=2024)%>%
  st_set_geometry("geometry")
  

year2 <- st_multipolygon(list(polygon3, polygon4, polygon5))%>%
  st_sfc() %>%
  st_union() %>%
  st_as_sf() %>%
  mutate(Year=2023)%>%
  st_set_geometry("geometry")

year3 <- st_multipolygon(list(polygon6, polygon7))%>%
  st_sfc() %>%
  st_union() %>%
  st_as_sf() %>%
  mutate(Year=2022) %>%
    st_set_geometry("geometry")




###combine into a single layer
combo_years <- bind_rows(list(year1, year2, year3)) 


#self intersect, give an index of overlaps "origins"
combo_years_intersect <- st_intersection(combo_years)

#use the index nubers to the get the year values from the original data set and concantenate into a list using purr functions, could alos use an apply or loop to do the same
combo_years_intersect_list <- combo_years_intersect %>%
  st_collection_extract("POLYGON") %>%  #this not necessary in simplified example but in messy shapes sometimes linestrings an points created
  mutate(Year_List=map_chr(origins, \(x) paste0(combo_years$Year[x], collapse = "; ") ))



#group and summarise by the nuber of overlaps and the the overalpping years
combo_years_intersect_summarised <- combo_years_intersect_list %>%group_by(n.overlaps, Year_List) %>%
  summarise() %>%
    ungroup() %>%
    mutate(area=st_area(.))




mapview::mapview(combo_years_intersect_summarised, zcol = "n.overlaps")



mapview::mapview(combo_years_intersect_summarised, zcol = "Year_List")



##can be piped together   



combo_years_summary_piped <- combo_years %>%
  st_intersection() %>%
  mutate(Year_List=map_chr(origins, \(x) paste0(combo_years$Year[x], collapse = "; ") ))%>%
  group_by(n.overlaps, Year_List) %>%
  summarise()%>%
  ungroup() %>%
  mutate(area=st_area(.))



##can calc other values at same time


combo_years_summary_many <- combo_years %>%
  st_intersection() %>%
  mutate(Year_List=map_chr(origins, \(x) paste0(combo_years$Year[x], collapse = "; ") ), 
         Year_Min=map_dbl(origins, \(x) min(combo_years$Year[x], na.rm=T) ),
         Year_Max=map_dbl(origins, \(x) min(combo_years$Year[x], na.rm=T) ))%>%
  group_by(n.overlaps, Year_List,Year_Min, Year_Max) %>%
  summarise()%>%
  ungroup() %>%
  mutate(area=st_area(.))




mapview::mapview(combo_years_summary_many, zcol = "n.overlaps")




mapview::mapview(combo_years_summary_many, zcol = "Year_Min")
  


##can now summarise total area by overlapping treatments 



OVERLAP_SUMMARY_AREAS <- combo_years_summary_many %>%
  st_drop_geometry() %>%
    group_by(n.overlaps) %>%
      summarise(total_area=sum(area))

#or

OVERLAP_SUMMARY_AREAS <- combo_years_summary_many %>%
  st_drop_geometry() %>%
  group_by(Year_List) %>%
  summarise(total_area=sum(area))