从2000年Blocks到2010年Tracts插值：计算重叠占比遇几何错误

需求说明 #

搭建从2000年Blocks到2010年Tracts的插值流程，生成每个Block在各Tract中的占比数据。例如Block A一半位于Tract X、一半位于Tract Y时，生成两条对应记录，占比列显示50%。

问题场景 #

参照相关示例操作时触发错误，已用st_is_valid验证所有多边形均为有效，但仍报错。

尝试代码 #

library("sf")
library("tigris")
library("dplyr")

maricopa_blocks00 <- blocks(
      "AZ", 
      "Maricopa", 
      year = 2000
    )
    
    
maricopa_tracts10 <- tracts(
  "AZ", 
  "Maricopa", 
  year = 2010
)

intersect_pct <- st_intersection(maricopa_tracts10, maricopa_blocks00) %>%
  mutate(intersect_area = st_area(.))

错误信息 #

Error in stopifnot(): ! Problem while computing intersect_area = st_area(.). Caused by error in wk_handle.wk_wkb(): ! Loop 0 is not
valid: Edge 106 is degenerate (duplicate vertex)

解决方法 #

错误源于细微拓扑问题（重复顶点/退化边），st_is_valid无法检测到这类问题，需用st_make_valid()修复。以下是修正后的完整流程：

修正代码 #

library("sf")
library("tigris")
library("dplyr")

# 加载并预处理Block数据：修复拓扑+预计算面积
maricopa_blocks00 <- blocks("AZ", "Maricopa", year = 2000) %>%
  st_make_valid() %>%
  mutate(block_area = st_area(.))

# 加载并预处理Tract数据：修复拓扑
maricopa_tracts10 <- tracts("AZ", "Maricopa", year = 2010) %>%
  st_make_valid()

# 确保两个图层坐标系一致（tigris默认一致，显式处理更稳妥）
maricopa_blocks00 <- st_transform(maricopa_blocks00, st_crs(maricopa_tracts10))

# 计算交集并生成占比列
intersect_pct <- st_intersection(maricopa_blocks00, maricopa_tracts10) %>%
  mutate(
    intersect_area = st_area(.),
    overlap_pct = as.numeric(intersect_area / block_area) * 100  # 转换为百分比格式
  ) %>%
  # 保留关键字段，可根据需求调整
  select(GEOID, GEOID10, intersect_area, block_area, overlap_pct)

关键说明 #

1. 拓扑修复：st_make_valid()自动处理重复顶点、退化边等细微拓扑问题，这是解决报错的核心步骤。
2. 预计算面积：提前计算Block的总面积，避免后续重复计算，提升运行效率。
3. 坐标系对齐：显式转换CRS确保两个图层空间参考一致，避免潜在的交集计算错误。
4. 占比计算：通过intersect_area / block_area得到占比，乘以100转换为百分比格式。

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