如何统计多年份杀虫剂处理多边形的不同覆盖次数对应面积？

核心思路 #

将所有年份的处理地块整合为带年份标识的空间数据集，通过空间相交操作分割出互不重叠的最小区域，统计每个区域被多少个年份覆盖，最后按覆盖次数分组计算总面积。

针对示例代码的实现步骤 #

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()
# ---------------------- 原示例代码结束 ----------------------

# 1. 将各年份数据转为带年份标识的sf对象
year1_sf <- st_sf(year = 1, geometry = year1)
year2_sf <- st_sf(year = 2, geometry = year2)
year3_sf <- st_sf(year = 3, geometry = year3)

# 2. 合并所有年份的空间数据
all_years <- bind_rows(year1_sf, year2_sf, year3_sf)

# 3. 确保几何有效性（避免后续相交操作报错）
all_years <- st_make_valid(all_years)

# 4. 分割为互不重叠的区域，自动生成覆盖次数（n.overlaps）
split_regions <- st_intersection(all_years)

# 5. 计算每个区域的面积，并按覆盖次数汇总总面积
coverage_summary <- split_regions |>
  mutate(area = st_area(geometry)) |>
  as_tibble() |>
  group_by(n.overlaps) |>
  summarise(total_area = sum(area)) |>
  rename(coverage_count = n.overlaps)

# 查看结果
print(coverage_summary)

适配8个年份shapefile的通用步骤 #

1. 用st_read()逐个读取每个年份的shapefile，例如：

   year_2015 <- st_read("path/to/year_2015.shp") |> mutate(year = 2015)
   year_2016 <- st_read("path/to/year_2016.shp") |> mutate(year = 2016)
   # 以此类推读取所有8个年份的数据
   

2. 用bind_rows()合并所有年份的sf对象；
3. 后续步骤与示例代码的第3-5步完全一致。

关键说明 #

• st_intersection()会自动将所有重叠区域分割为互不重叠的小块，并通过n.overlaps列标记每个小块被多少个原始多边形（即多少个年份）覆盖；
• st_make_valid()是必要步骤，可修复shapefile中常见的几何无效问题（如自相交、环方向错误），避免后续操作报错；
• 最终结果中的coverage_count就是覆盖次数（1次、2次、3次...），total_area对应该次数的总面积。

内容的提问来源于stack exchange，提问作者Sara