如何用R的sf包创建几何列为两数据框几何差的空间数据框？

问题描述 #

我用R的sf包加载了两个空间数据框：第一个包含61个美国县的空间数据，第二个包含每个县对应的子区域数据。两个数据框有相同的识别列及对应值，每个县都有唯一对应的子区域。我想要得到第三个空间数据框，其几何列为县与对应子区域的几何差。

尝试代码如下：

library(tidyverse)
library(magrittr)
library(sf)

# 加载shapefile
dfc <- st_read('data/counties.shp') %>%      # 包含61个美国县的shapefile
  arrange(., st, cty)                        # 用于唯一识别观测的列
dfp <- st_read('data/county_parts.shp') %>%  # 包含61个美国县子区域的shapefile
  arrange(., st, cty)                        # 用于唯一识别观测的列

# 计算每对对应观测的几何差
geom <- map(seq(nrow(dfc)), function(r) 
            st_difference(dfc[r, 'geo_cty'], dfp[r, 'geo_sub'])

（注：geo_cty为县几何列，geo_sub为县子区域几何列。按st和cty排序后，dfc与dfp的每行r对应相同的(st, cty)值，且geo_cty包含geo_sub。）

这段代码能正常运行，但无法将生成的geom列表转为空间数据框的几何列。尝试以下合并代码时报错：

df <- bind_cols(dfc, select(dfp, geo_sub), geom)

错误信息：

Error in `stop_vctrs()`:
! Can't recycle `..1` (size 61) to match `..58` (size 0).
Run `rlang::last_error()` to see where the error occurred.

解决方法 #

问题根源 #

map生成的geom是sf对象的列表，而bind_cols无法直接将列表与sf数据框合并——sf数据框的几何列是特殊的sfc（简单要素几何列）类型，并非普通列表，直接合并会触发维度不匹配的错误。

最佳实现方式：向量化运算 #

sf包支持向量化空间运算，无需循环即可直接对整列计算几何差，全程保持sf数据框结构，代码更简洁高效：

library(tidyverse)
library(sf)

# 加载并排序数据，确保行序对应
dfc <- st_read('data/counties.shp') %>% 
  arrange(st, cty)
dfp <- st_read('data/county_parts.shp') %>% 
  arrange(st, cty)

# 生成包含几何差的目标空间数据框
df_result <- dfc %>%
  # 保留需要的属性列（可根据需求调整）
  select(st, cty, geo_cty) %>%
  # 匹配添加子区域几何列
  mutate(geo_sub = dfp$geo_sub) %>%
  # 向量化计算几何差
  mutate(geo_diff = st_difference(geo_cty, geo_sub)) %>%
  # 将几何差列设为数据框的活跃几何列
  st_set_geometry("geo_diff") %>%
  # 可选：删除不需要的中间列
  select(-geo_cty, -geo_sub)

循环方式的兼容处理 #

如果坚持使用原map逻辑，需先将列表转为sfc类型，再整合到sf数据框：

# 补全原map代码的闭合括号（原代码缺少闭合）
geom <- map(seq(nrow(dfc)), function(r) 
            st_difference(dfc[r, 'geo_cty'], dfp[r, 'geo_sub']))

# 将列表转为sfc几何列
geom_sfc <- do.call(c, geom)

# 构建结果空间数据框
df_result <- dfc %>%
  mutate(geo_diff = geom_sfc) %>%
  st_set_geometry("geo_diff") %>%
  select(st, cty, geo_diff) # 保留需要的列

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