土地覆被SHP整合道路缓冲及旧道路多边形替换技术问询

步骤1：数据读取与预处理（补充完整流程） #

假设你已完成道路缓冲区创建，以下是完整数据加载与缓冲区生成代码：

library(terra)

# 读取土地覆被和道路线矢量
lc <- vect("land_cover.shp")
roads <- vect("roads.shp")

# 生成道路缓冲区（按实际需求调整宽度）
road_buf <- buffer(roads, width = 10)
# 给缓冲区添加统一类别字段（与土地覆被字段名保持一致）
road_buf$class <- "new road"

步骤2：整合新道路缓冲区到土地覆被 #

你之前合并失败的核心原因是未过滤旧道路多边形，导致新旧道路重叠。正确做法是先移除旧道路，再合并新缓冲区：

# 过滤土地覆被中的旧道路（假设旧道路类别值为"road"）
lc_no_road <- lc[lc$class != "road", ]

# 统一字段名（若土地覆被类别字段名不同，需先调整匹配）
# names(road_buf)[names(road_buf) == "class"] <- "你的土地覆被类别字段名"

# 合并过滤后的土地覆被与新道路缓冲区
lc_combined <- rbind(lc_no_road, road_buf)

步骤3：模拟ArcGIS Eliminate工具替换旧道路多边形 #

要实现删除旧道路并替换为邻接最长的土地覆被，需按以下流程处理：

# 提取旧道路多边形
old_roads <- lc[lc$class == "road", ]
# 提取所有非道路多边形
non_road <- lc[lc$class != "road", ]

# 初始化合并结果容器
merged_lc <- non_road

# 遍历每个旧道路多边形
for (i in 1:nrow(old_roads)) {
  current_road <- old_roads[i, ]
  # 找到与当前道路邻接的多边形
  neighbors <- adjacent(non_road, current_road, "touching", pairs = FALSE)
  
  if (length(neighbors) > 0) {
    neighbor_polys <- non_road[neighbors, ]
    # 计算每个邻接多边形与旧道路的共享边界长度
    boundary_lengths <- perimeter(intersect(neighbor_polys, current_road))
    # 选择边界最长的邻接多边形
    max_idx <- which.max(boundary_lengths)
    target_poly <- neighbor_polys[max_idx, ]
    
    # 合并旧道路到目标多边形，保留原类别属性
    merged_poly <- union(target_poly, current_road)
    merged_poly$class <- target_poly$class
    
    # 更新合并结果：移除原目标多边形，添加合并后的新多边形
    merged_lc <- merged_lc[-which(merged_lc$ID == target_poly$ID), ] # 替换为你的唯一标识字段
    merged_lc <- rbind(merged_lc, merged_poly)
  } else {
    warning(paste("旧道路多边形", i, "无邻接多边形，已跳过"))
  }
}

# 最终整合新道路缓冲区与替换后的土地覆被
final_lc <- rbind(merged_lc, road_buf)

关键注意事项 #

• 确保所有矢量数据坐标系一致，可用crs()检查，project()转换
• 替换代码中的class、ID字段为你数据中的实际字段名
• 缓冲区宽度需匹配数据比例尺，避免出现不合理的道路范围

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