如何高效计算空间点到最近道路线的最短距离？（rgdal/rgeos已弃用）

现在rgdal、rgeos已被弃用，推荐使用sf包（当前空间数据处理的标准工具）来完成这个任务，针对9000个点的规模，以下方法高效且易实现：

步骤1：转换数据格式 #

首先将SpatialLinesDataFrame转换为sf对象，sf包的st_as_sf()可以直接完成转换：

library(sf)

# 将道路数据转为sf对象
roads_sf <- st_as_sf(your_spatial_lines_df)
# 你的点数据已经是sf对象，记为points_sf

步骤2：统一坐标系 #

确保点和道路数据的坐标系一致，否则距离计算结果会出错：

# 检查坐标系
st_crs(points_sf)
st_crs(roads_sf)

# 如果不一致，将道路数据转换为点数据的坐标系
if (st_crs(points_sf) != st_crs(roads_sf)) {
  roads_sf <- st_transform(roads_sf, st_crs(points_sf))
}

步骤3：高效计算最近距离 #

直接计算所有点到所有道路的距离会生成大矩阵，占用内存且慢。我们可以先找到每个点对应的最近道路，再计算该点到这条道路的距离：

# 找到每个点最近的道路的索引
nearest_road_idx <- st_nearest_feature(points_sf, roads_sf)

# 计算每个点到对应最近道路的距离（by_element=TRUE确保一对一计算）
points_sf$distance_to_road <- st_distance(points_sf, roads_sf[nearest_road_idx, ], by_element = TRUE)

注意：距离单位转换 #

如果你的数据用的是WGS84（EPSG:4326，经纬度坐标系），st_distance返回的距离单位是度，这没有实际意义。需要先转换为投影坐标系（如UTM分区），再计算米为单位的距离：

# 选择对应区域的UTM EPSG代码（示例为北半球33区，根据你的数据位置调整）
utm_epsg <- 32633

points_utm <- st_transform(points_sf, utm_epsg)
roads_utm <- st_transform(roads_sf, utm_epsg)

nearest_road_idx <- st_nearest_feature(points_utm, roads_utm)
# 转换为数值型，方便后续处理
points_sf$distance_m <- as.numeric(st_distance(points_utm, roads_utm[nearest_road_idx, ], by_element = TRUE))

sf包内置空间索引，st_nearest_feature会利用索引快速定位最近要素，9000个点的计算几乎瞬间完成，完全不用担心效率问题。

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