如何在多边形范围内计算点到城市的最近可达距离？

问题描述 #

我需要找到某个体在只能沿多边形区域内路径移动的约束下，距离最近的城市。现有方法sf::st_nearest_feature()仅计算直线距离，不符合需求：

• 直线距离下，个体（红色三角形）到城市A更近；
• 但限制路径在多边形内时，个体到城市B更近。

示例代码 #

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(ggplot2)
library(ggrepel)
library(rnaturalearth)

background <- ne_countries(scale = 'small', type = 'map_units', returnclass = 'sf') |>
  subset(name %in% c("England", "Wales")) |> 
  st_union()

cities <- data.frame(
  name = c("A", "B"),
  lon = c(-4.3, -3.3),
  lat = c(51.2, 51.45)
) |> 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

individual <- data.frame(id = 1, lon = -4.3, lat = 51.6) |> 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

ggplot() +
  geom_sf(data = background) +
  geom_sf(data = cities, size = 3) +
  geom_sf(data = individual, color = "red", shape = 17, size = 3) +
  coord_sf(xlim = c(-6, -1), ylim = c(50, 52)) +
  geom_text_repel(
    data = cities, 
    aes(geometry = geometry, label = name),
    stat = "sf_coordinates",
  )
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data

当前直线距离结果 #

nearest <- st_nearest_feature(individual, cities)
cities[nearest, "name"]
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -4.3 ymin: 51.2 xmax: -4.3 ymax: 51.2
#> Geodetic CRS:  WGS 84
#>   name          geometry
#> 1    A POINT (-4.3 51.2)

需求：修改距离计算逻辑，得到正确的最近城市B，同时方案需高效适配百万级个体（城市数量较少）的计算场景。

解决方案 #

这个问题本质是计算多边形区域内的最短路径距离，而非直线距离。推荐使用gdistance包（专门处理空间最短路径）结合投影优化来实现高效计算，步骤如下：

1. 环境准备与投影转换 #

地理坐标系（WGS84）不适合距离计算，先转换为英国本地平面坐标系（EPSG:27700，OSGB36），提升计算精度和效率：

library(gdistance)

# 转换投影到平面坐标系
background_proj <- st_transform(background, crs = 27700)
cities_proj <- st_transform(cities, crs = 27700)
individual_proj <- st_transform(individual, crs = 27700)

2. 创建可通行区域的过渡矩阵 #

将多边形区域转换为栅格，标记可通行区域（多边形内为1，外为NA），生成过渡矩阵用于最短路径计算：

# 创建栅格（分辨率可根据需求调整，分辨率越低计算越快）
r <- raster(background_proj, res = 1000)  # 1000米分辨率
r <- mask(r, background_proj)  # 仅保留多边形内区域

# 生成过渡矩阵（8方向移动，权重为距离）
trans <- transition(r, function(x) 1, 8)
trans <- geoCorrection(trans, type = "c")  # 校正地理距离

3. 计算个体到所有城市的区域内最短路径距离 #

由于城市数量少，对每个个体计算到所有城市的最短路径长度，取最小值对应的城市：

# 提取城市和个体的坐标（栅格坐标）
city_coords <- st_coordinates(cities_proj)
ind_coord <- st_coordinates(individual_proj)

# 计算个体到每个城市的最短路径距离
distances <- apply(city_coords, 1, function(city) {
  shortestPath(trans, ind_coord, city, output = "distance")
})

# 找到距离最近的城市
nearest_city <- cities[which.min(distances), ]
nearest_city$name
#> [1] "B"

4. 百万级个体的优化方案 #

针对百万级个体，可通过以下方式提升效率：

• 预缓存城市坐标：提前转换所有城市的栅格坐标，避免重复计算；
• 并行计算：使用future.apply或foreach包对个体批量并行计算，充分利用多核CPU；
• 降低栅格分辨率：如果对精度要求不高，可增大栅格分辨率（如2000米），大幅减少计算量；
• 空间索引：用st_join先将个体按空间范围分组，仅计算组内邻近城市的距离，减少不必要的计算。

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