使用st_interpolate_aw时丢失ID：输出多边形比目标少1个的问题

问题背景 #

我拥有两个经st_make_valid()处理、st_is_valid()验证全部有效的sf对象：

• soilshape.sf.drop：2020年带土壤属性数据的矢量图层，维度为693×11
• dist1970.sf：1970年的边界矢量图层，维度为386×3，所有ID均为唯一值

尝试执行以下代码，期望将dist1970.sf的ID关联到插值结果时出现报错：

atest1 = st_interpolate_aw(soilshape.sf["Zn_s"], dist1970.sf, extensive = FALSE) |>
  mutate(ID = dist1970.sf$ID)

报错信息：

Error in `stopifnot()`:
ℹ In argument: `ID = dist1970.sf$ID`.
Caused by error:
! `ID` must be size 385 or 1, not 386.
Run `rlang::last_trace()` to see where the error occurred.

排查发现单独执行插值代码得到的atest2对象维度为385×2，比dist1970.sf少1个多边形。需要解决：

• 为何st_interpolate_aw输出的多边形数量与目标对象不符？
• 如何保留原ID，使插值后的数据可与原图层匹配使用？

原因分析 #

st_interpolate_aw输出数量减少的核心原因是：目标图层dist1970.sf中存在一个与源图层soilshape.sf.drop完全无空间交集的多边形。这类多边形因无法计算权重值（无重叠区域无法分配土壤数据），会被自动从插值结果中剔除。

解决步骤 #

1. 定位丢失的多边形 #

先找出dist1970.sf中与源图层无交集的多边形：

# 计算每个目标多边形与源图层总范围的交集面积
intersect_check = dist1970.sf |>
  mutate(intersect_area = st_area(st_intersection(., st_union(soilshape.sf.drop))))

# 筛选出交集面积为0的多边形
missing_polygon = intersect_check[intersect_check$intersect_area == 0, ]

查看missing_polygon的ID即可确认被剔除的具体多边形。

2. 保留全部多边形并关联原ID #

不要直接用mutate强制绑定ID，而是通过空间左连接保留所有原多边形，确保插值结果与原图层一一对应：

# 执行插值得到385行的结果
interpolated_data = st_interpolate_aw(soilshape.sf["Zn_s"], dist1970.sf, extensive = FALSE)

# 基于几何对象左连接，保留原图层所有386个多边形
atest1 = dist1970.sf |>
  left_join(interpolated_data, by = st_geometry_name(dist1970.sf))

处理后atest1会保持386行，与原dist1970.sf完全匹配，无交集的多边形对应的Zn_s值会显示为NA，同时完整保留原ID列。

3. 可选：填充无交集区域的插值值 #

如果需要给无交集的多边形设置默认值，可按需求填充（示例为填充土壤数据的全局均值）：

atest1 = atest1 |>
  mutate(Zn_s = ifelse(is.na(Zn_s), mean(soilshape.sf$Zn_s, na.rm = TRUE), Zn_s))

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