使用st_intersection聚合空间多边形数据失败求助

我在将印度1970年区县shapefile与带土壤数据的2020年区县shapefile进行空间聚合时遇到困难。已通过readOGR()读取文件并转换为sf对象（soilshape.sf为2020年带土壤变量的区县，dist1970.sf为1970年区县），能正常运行以下代码查看空间重叠情况：

# 成功绘制1970区县叠加在2020区县上的地图
ggplot(data = soilshape.sf) +
  geom_sf(aes(fill = Zn_s)) +
  geom_sf(data = dist1970.sf, lwd = .5, color = "#000000", fill = NA)

但在将2020年土壤数据按面积权重聚合到1970年区县时，出现两个异常：

• 386个1970年区县中有105个缺失土壤信息（sum(is.na(dist1970.zn@data$wZn_s))返回结果为105）
• 绘制聚合后地图时报错，运行代码：

ggplot(data = dist1970.zn.sf) +
  geom_sf(aes(fill = wZn_s)) 

得到错误信息：

Error in `geom_sf()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 1st layer.
Caused by error in `scale_type.units()`:
! Variable of class 'units' found, but 'units' package is not attached.
  Please, attach it using 'library(units)' to properly show scales with units.
Run `rlang::last_trace()` to see where the error occurred.

我使用的聚合代码如下（已尝试st_make_valid()修复多边形、切换sf_use_s2(TRUE/FALSE)，问题依旧）：

# 空间叠加找重叠区域
sf_use_s2(FALSE)
overlay <- st_intersection(soilshape.sf, dist1970.sf)

# 计算每个交集的面积
overlay$intersection_area <- st_area(overlay)

# 按1970年区县ID分组计算加权均值
result <- overlay %>%
  dplyr::group_by(ID) %>%
  dplyr::summarize(wZn_s = sum(Zn_s * intersection_area) / sum(intersection_area))

# 合并回1970年区县的sp对象再转sf
dist1970.zn <- merge(dist1970, result, by.x ='ID', by.y ='ID', all.x=TRUE)
dist1970.zn.sf <- st_as_sf(dist1970.zn)

1. 地图绘制错误的快速修复 #

报错核心原因是wZn_s带有面积单位（st_area()返回带单位的数值），ggplot无法识别该类型。有两种修复方式：

• 直接加载units包：

  library(units)
  # 重新运行绘图代码
  ggplot(data = dist1970.zn.sf) +
    geom_sf(aes(fill = wZn_s)) 
  

• 计算加权均值时去除单位（转换为纯数值）：

  result <- overlay %>%
    dplyr::group_by(ID) %>%
    dplyr::summarize(wZn_s = sum(Zn_s * as.numeric(intersection_area)) / sum(as.numeric(intersection_area)))
  

2. 大量区县缺失数据的原因排查与解决 #

原因1：1970年区县与2020年区县无空间交集 #

部分1970年区县可能完全不在2020年区县的覆盖范围内，导致无交集数据。可以通过绘图验证：

# 找出无匹配的ID
missing_ids <- dist1970.sf$ID[!dist1970.sf$ID %in% result$ID]
# 绘制这些区县，查看是否在2020区县范围外
ggplot(data = dist1970.sf %>% filter(ID %in% missing_ids)) +
  geom_sf(fill = "red") +
  geom_sf(data = soilshape.sf, fill = NA, color = "gray")

如果确实存在此类区县，需要补充对应区域的2020年土壤数据。

原因2：2020年区县的土壤数据本身存在NA #

检查soilshape.sf中Zn_s的缺失情况：

sum(is.na(soilshape.sf$Zn_s))

如果Zn_s本身有大量缺失，交集区域的加权计算自然会得到NA。需要先处理2020年土壤数据的缺失值（如插值、补全）。

原因3：空间叠加的拓扑或精度问题 #

即使使用st_make_valid()，仍可能存在细微拓扑错误。可以尝试用st_join结合面积加权的替代方案：

# 关联1970区县与所有相交的2020区县
joined <- st_join(dist1970.sf, soilshape.sf, join = st_intersects)
# 计算每个交集区域占1970区县的面积比例
joined$area_prop <- as.numeric(st_area(st_intersection(joined, soilshape.sf)) / st_area(joined))
# 按ID分组计算加权均值，自动忽略NA
result_new <- joined %>%
  dplyr::group_by(ID) %>%
  dplyr::summarize(wZn_s = weighted.mean(Zn_s, area_prop, na.rm = TRUE))
# 合并回原sf对象
dist1970.zn.sf <- dplyr::left_join(dist1970.sf, result_new, by = "ID")

原因4：1970年区县合并时的拓扑错误 #

原始代码使用sp包的unionSpatialPolygons合并多边形，可能产生无效拓扑。建议改用sf包的原生方法重新合并：

# 重新读取1970年原始数据并合并
setwd("/Users/bevis.16/Library/CloudStorage/OneDrive-TheOhioStateUniversity/BoxLeah/Projects/GR Inequality/Data/Geospatial")
distVDS <- readOGR(dsn="Districts_For_VDS", layer="india70again") %>% st_as_sf()

# 用sf的st_union合并同一ID的多边形
dist1970.sf <- distVDS %>%
  dplyr::group_by(DISTCODE) %>%
  dplyr::summarize(geometry = st_union(geometry)) %>%
  dplyr::rename(ID = DISTCODE) %>%
  st_set_crs(st_crs(soilshape.sf))

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