R语言sf包st_difference处理投影空间对象异常问题求助

问题根源 #

你遇到的问题本质是地理坐标系与平面坐标系的拓扑运算逻辑差异：

• 当你设置CRS为4326（WGS84经纬度）时，sf会把数据当作球面/椭球面上的地理要素处理，拓扑运算（st_difference、st_intersection）会基于球面几何计算，这类计算容易出现精度偏差，尤其是线与多边形的裁剪场景，结果会偏离直观预期。
• 无CRS时，sf默认按平面笛卡尔坐标系处理，拓扑运算基于平面几何，结果符合你预期的裁剪效果。

解决方法 #

把地理坐标系的要素转换为**投影坐标系（平面坐标系）**完成拓扑运算，之后可按需转回原坐标系。

操作步骤 #

1. 选择适配数据的投影坐标系：针对美国本土数据，推荐用EPSG:5070（NAD83 Conus Albers，等面积投影），或者根据区域选对应UTM分区（比如EPSG:26910对应美国西部），EPSG:3857（Web墨卡托）也可通用但距离有轻微失真。
2. 用st_transform()转换坐标系，完成运算后再转回原CRS（如果需要经纬度格式）。

修正后的测试代码 #

if(!"pacman" %in% installed.packages()[ , "Package"])
  install.packages("pacman")
require(pacman)
p_load(dplyr, sf, ggplot2, tidyverse, units)

# 创建测试要素
m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))
p = st_polygon(list(m)) %>% st_sfc()
m = rbind(c(0.5,0.5), c(0.5,1.5))
l1 = st_linestring(m) %>% st_sfc()

# 设置地理坐标系
st_crs(p) <- 4326
st_crs(l1) <- 4326

# 转换为投影坐标系（示例用EPSG:5070）
p_proj <- st_transform(p, 5070)
l1_proj <- st_transform(l1, 5070)

# 执行拓扑运算
inter1_proj <- st_intersection(l1_proj, p_proj)
diff_proj <- st_difference(l1_proj, inter1_proj)

# 转回原地理坐标系（按需选择）
inter1 <- st_transform(inter1_proj, 4326)
diff <- st_transform(diff_proj, 4326)

# 绘图验证
ggplot(p)+
  geom_sf()+
  geom_sf(data = l1, col = "blue")+
  geom_sf(data = inter1, col = "red", lwd = 3, alpha = 0.25, lty = 1)+
  geom_sf(data = diff, col = "green", lwd = 2, alpha = 0.5, lty = 5)

针对NHGIS实际数据的建议 #

• 计算县内铁路长度时，务必在投影坐标系下用st_length()计算，这样得到的是真实的平面距离，避免经纬度直接计算长度的失真问题。
• 若需要按县汇总，可结合group_by()和sum()对投影后的铁路长度进行汇总，结果更准确。

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