在R中处理两个不同时间点sf对象的重叠与非重叠多边形

核心问题分析 #

出现GEOS exception是因为st_sym_difference生成的要素存在无效几何（比如自相交、零面积多边形等拓扑问题），直接调用无参数的st_intersection()处理内部重叠时，GEOS无法解析这些错误拓扑。

方案一：修复几何后分区域合并 #

这是最稳妥的方式，分三步生成包含所有区域的结果：

1. 预处理：修复几何并消除内部重叠
   先对原始图层做几何有效性修复，再处理内部重叠：

options(stringsasfactors = FALSE)
library(ggplot2)
library(sf)
library(dplyr)
library(tidyr)

# 读取并修复几何
t1 <- st_read('T1.gpkg') %>% st_make_valid()
t2 <- st_read('T2.gpkg') %>% st_make_valid()

# 消除图层内部重叠
t1_no_overlap <- st_intersection(t1)
t2_no_overlap <- st_intersection(t2)

2. 拆分出所有独立区域
   通过合并两个图层的边界，拆分出所有唯一的多边形单元：

# 合并两个图层的几何边界
combined_boundary <- st_union(t1_no_overlap, t2_no_overlap)
# 拆分为单个多边形
all_regions <- st_cast(combined_boundary, "POLYGON") %>% st_sf()

3. 关联两个时间点的属性
   对每个拆分后的区域，分别匹配t1和t2的属性，无匹配的设为NA：

# 匹配t1属性
all_regions <- all_regions %>%
  st_join(t1_no_overlap %>% rename(t1_cat = cat), left = TRUE) %>%
  group_by(geometry) %>%
  slice(1) %>% # 去重，避免同一区域匹配多个t1要素
  ungroup()

# 匹配t2属性
all_regions <- all_regions %>%
  st_join(t2_no_overlap %>% rename(t2_cat = cat), left = TRUE) %>%
  group_by(geometry) %>%
  slice(1) %>%
  ungroup()

方案二：分区域计算后合并 #

直接计算交集、t1独有区、t2独有区，再合并结果：

# 预处理（同方案一）
t1 <- st_read('T1.gpkg') %>% st_make_valid() %>% st_intersection()
t2 <- st_read('T2.gpkg') %>% st_make_valid() %>% st_intersection()

# 计算重叠区域，保留双属性
overlap <- st_intersection(t1, t2) %>%
  mutate(t1_cat = cat.x, t2_cat = cat.y) %>%
  select(t1_cat, t2_cat, geometry)

# 计算t1独有区域，t2属性设NA
t1_only <- st_difference(t1, st_union(t2)) %>%
  mutate(t1_cat = cat, t2_cat = NA) %>%
  select(t1_cat, t2_cat, geometry)

# 计算t2独有区域，t1属性设NA
t2_only <- st_difference(t2, st_union(t1)) %>%
  mutate(t1_cat = NA, t2_cat = cat) %>%
  select(t1_cat, t2_cat, geometry)

# 合并所有区域
final_result <- bind_rows(overlap, t1_only, t2_only) %>% st_make_valid()

验证结果 #

可以用ggplot查看最终结果：

ggplot(final_result) +
  geom_sf(aes(fill = interaction(t1_cat, t2_cat)), alpha = 0.7) +
  labs(fill = "T1 & T2 属性组合")

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