迁移至sf包后st_intersection/st_difference无法正确裁剪多边形

问题描述 #

需要裁剪/差分超出边界的多边形，因rgdal包即将停用，正迁移至sf包。原rgdal代码运行正常，可生成裁剪后的多边形效果；但改用sf包后，无论是使用st_intersection还是st_difference，均无法正确计算差分，多边形未被裁剪。期望达到与rgdal版本相同的效果。

原rgdal代码 #

WorldCountry <- sf::st_read("data/modified.countries.geo.json")
countries <- c("Canada", "Mexico")
data_Map <- WorldCountry[WorldCountry$name %in% countries, ]
data_Map <- as_Spatial(data_Map)
for (i in seq(length(spa_polys), 3)) {
   spa_polys@proj4string <- data_Map@proj4string
   geo_diff <- gDifference(spa_polys[i], data_Map)
   iso@polygons[[i]] <- geo_diff@polygons[[1]]
   iso@polygons[[i]]@ID = as.character(i)
   rm(geo_diff)
 }

rgdal版本效果：可正确裁剪超出边界的多边形

现有sf包代码 #

WorldCountry <- sf::st_read("data/modified.countries.geo.json")
countries <- c("Canada", "Mexico")
data_Map <- WorldCountry[WorldCountry$name %in% countries, ]
data_Map <- as_Spatial(data_Map)
data_Map_sf <- st_as_sf(data_Map, coords = c("longitude", "latitude"), crs = 4326)

for (i in seq(length(spa_polys), 1)) {
  spa_poly_sf <- st_as_sf(spa_polys[i])
  st_crs(spa_poly_sf) <- st_crs(data_Map_sf)
  if (!st_is_longlat(spa_poly_sf)) {
    spa_poly_sf <- st_transform(spa_poly_sf, st_crs(data_Map_sf))
  }
  geo_cut <- st_intersection(spa_poly_sf, data_Map_sf)
  iso$geometry[[i]] <- geo_cut
  iso$ID[i] <- as.character(i)
  rm(geo_cut)
}

sf版本效果：无法正确裁剪，多边形未被处理

解决建议 #

• 避免格式反复转换：原代码中把data_Map转成sp格式再转回sf属于多余操作，直接保留sf格式即可，减少转换误差：

  data_Map_sf <- WorldCountry[WorldCountry$name %in% countries, ]
  # 无需转成as_Spatial再转回sf
  

• 确保CRS匹配正确：不要强行赋值CRS，先确认原始数据的坐标系，再统一转换：

  spa_poly_sf <- st_as_sf(spa_polys[i])
  # 先设置spa_polys的原始CRS（示例为4326，需替换为实际值）
  if (is.na(st_crs(spa_poly_sf))) {
    st_crs(spa_poly_sf) <- 4326
  }
  # 统一转换到与data_Map_sf一致的CRS
  spa_poly_sf <- st_transform(spa_poly_sf, st_crs(data_Map_sf))
  

• 对齐原代码逻辑：原rgdal用的是gDifference（差分，即裁剪掉重叠部分），而sf代码误用了st_intersection（交集，保留重叠部分），需替换为st_difference：

  geo_cut <- st_difference(spa_poly_sf, data_Map_sf)
  

• 处理特殊几何情况：差分后可能生成空几何或多部分几何，需判断后再赋值：

  if (!st_is_empty(geo_cut)) {
    # 若为多部分多边形，可选择保留第一个部分或全部
    if (st_geometry_type(geo_cut) == "MULTIPOLYGON") {
      iso$geometry[[i]] <- st_cast(geo_cut, "POLYGON")[1]
    } else {
      iso$geometry[[i]] <- geo_cut$geometry
    }
  }
  

• 修正循环序列：原rgdal代码循环范围是seq(length(spa_polys), 3)，sf代码写成了seq(length(spa_polys), 1)，索引范围不一致，需保持和原代码一致的循环逻辑。

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