在R中计算不规则多边形内两点不越界的最短路径距离

要计算不规则多边形内部两点间不穿越边界的最短路径，你可以借助sf生态下的空间网络或成本表面工具实现，以下是两种可行方案：

方案1：基于空间网络（沿边界的最短路径） #

这种方法将多边形边界转为无向网络，把内部点融合到网络的最近边界节点，再用图算法计算最短路径，路径会严格沿多边形边界延伸，不会穿越边界。

library(sf)
library(sfnetworks)
library(igraph)

# 修正多边形创建（转为POLYGON类型）
mypolygon = data.frame(x = c(-5,-5, 5, 5,3,-3,4, -5),  
                       y = c(-5,5,5,-5,0,0,-4, -5))
p_poly = st_as_sf(mypolygon, coords = c("x", "y")) %>%
  summarise(do_union = FALSE) %>%
  st_cast("POLYGON")

# 创建内部两点
point_df <- data.frame(x = c(2,4.8), y = c(-4,-4))
points_sf <- st_as_sf(point_df, coords = c("x", "y"), crs = st_crs(p_poly))

# 构建多边形边界网络
poly_boundary <- st_cast(p_poly, "LINESTRING")
network <- as_sfnetwork(poly_boundary, directed = FALSE)

# 将内部点融合到网络（匹配最近边界点）
network_with_points <- st_network_blend(network, points_sf)

# 提取两点对应的网络节点ID
node_indices <- V(network_with_points)[st_is_within_distance(points_sf, network_with_points, dist = 1e-6)]

# 计算最短路径
shortest_path <- shortest_paths(network_with_points, node_indices[1], node_indices[2], output = "path")

# 转换路径为sf对象并绘图
path_sf <- st_as_sf(network_with_points, "edges")[unlist(shortest_path$vpath[[1]][-1]), ] %>%
  st_combine() %>%
  st_cast("LINESTRING")

# 可视化结果
plot(st_geometry(p_poly), border = "black", lwd = 2)
plot(st_geometry(points_sf), pch = 19, cex = 2, col = "red", add = TRUE)
plot(st_geometry(path_sf), col = "yellow", lwd = 3, add = TRUE)

方案2：基于成本表面（平滑内部路径） #

如果需要更平滑的内部路径（不严格沿边界），可以用gdistance包创建成本表面，将多边形外部设为不可通行区域，再计算两点间的最短路径。

library(sf)
library(gdistance)
library(raster)

# 复用之前的多边形和两点对象
# 创建栅格（分辨率可根据精度需求调整）
r <- raster(extent(p_poly), res = 0.1)
r[] <- 1

# 将多边形外部区域设为不可通行（成本无穷大）
poly_mask <- st_as_sfc(st_bbox(p_poly)) %>%
  st_difference(p_poly)
r[poly_mask] <- Inf

# 创建过渡层并校正地理距离
tr <- transition(r, function(x) 1/mean(x), 8)
tr <- geoCorrection(tr)

# 计算最短路径并转为sf对象
shortest_path_gd <- shortestPath(tr, points_sf[1,], points_sf[2,], output = "SpatialLines")
shortest_path_gd_sf <- st_as_sf(shortest_path_gd)

# 可视化结果
plot(st_geometry(p_poly), border = "black", lwd = 2)
plot(st_geometry(points_sf), pch = 19, cex = 2, col = "red", add = TRUE)
plot(st_geometry(shortest_path_gd_sf), col = "yellow", lwd = 3, add = TRUE)

两种方案对比 #

• 空间网络方案：路径严格贴合多边形边界，适合需要精确沿边界移动的场景，计算速度快。
• 成本表面方案：路径更平滑，可在多边形内部自由通行，适合允许内部路径的场景，精度由栅格分辨率决定。

内容的提问来源于stack exchange，提问作者Alejandro Molina Moctezuma