基于Shapefile的R语言道路网络两点驾驶距离计算报错解决

需求说明 #

需要计算以下两点的路网实际驾驶距离：

• 加拿大国家电视塔（CN Tower）：坐标(-79.61203, 43.68312)
• 多伦多皮尔逊机场：坐标(-79.38709, 43.64256)

数据准备（已完成步骤） #

已下载加拿大道路网Shapefile并筛选出安大略省数据，代码如下：

library(sf)
library(sfnetworks)
# 下载并读取安大略省道路网
url <- "https://www12.statcan.gc.ca/census-recensement/2011/geo/RNF-FRR/files-fichiers/lrnf000r22a_e.zip"
temp_dir <- tempdir()
temp_file <- file.path(temp_dir, "lrnf000r22a_e.zip")
download.file(url, temp_file)
unzip(temp_file, exdir = temp_dir)
# 筛选安大略省（PRUID_R='35'）数据
ontario_roads <- st_read(file.path(temp_dir, "lrnf000r22a_e.shp"), 
                        query="select * from lrnf000r22a_e where PRUID_R ='35'")

注：路网数据的投影为NAD83 / Statistics Canada Lambert（单位：米）

原代码报错问题 #

原计算代码执行时触发错误：

Error in UseMethod("st_geometry") : no applicable method for 'st_geometry' applied to an object of class "c('tbl_df', 'tbl', 'data.frame')"

原代码片段：

net <- as_sfnetwork(ontario_roads)
# 定义起点终点
q1 <- st_point(c(-79.61203, 43.68312))
q2 <- st_point(c(-79.38709, 43.64256))
q1 <- st_sfc(q1, crs = st_crs(ontario_roads))
q2 <- st_sfc(q2, crs = st_crs(ontario_roads))
path <- st_network_paths(net, q1, q2)
distance <- sum(st_length(path))

错误原因 #

1. CRS不匹配：输入坐标为WGS84经纬度（EPSG:4326），直接设置为路网投影会导致坐标定位错误，需先创建WGS84坐标系的点，再转换到路网投影。
2. 路径结果处理错误：st_network_paths返回的是包含路径节点/边ID的DataFrame，不是几何对象，无法直接用st_length计算长度。

修正后的完整代码 #

library(sf)
library(sfnetworks)
library(tidygraph)

# 1. 构建有效道路网络（修复拓扑问题）
net <- as_sfnetwork(ontario_roads, directed = FALSE) %>% 
  convert(to_spatial_smooth)

# 2. 正确转换起点终点坐标
# 创建WGS84坐标系的点
start_point <- st_sfc(st_point(c(-79.61203, 43.68312)), crs = 4326)
end_point <- st_sfc(st_point(c(-79.38709, 43.64256)), crs = 4326)
# 转换到路网的投影坐标系
start_proj <- st_transform(start_point, st_crs(ontario_roads))
end_proj <- st_transform(end_point, st_crs(ontario_roads))

# 3. 计算基于道路长度的最短路径
path_result <- st_network_paths(net, 
                                from = start_proj, 
                                to = end_proj,
                                weights = "length")

# 4. 提取路径几何并计算总距离
# 从网络中获取路径对应的边
path_edges <- slice(activate(net, "edges"), path_result$edge_path[[1]])
# 计算总长度（单位：米）
total_distance <- sum(st_length(path_edges))
# 转换为公里（可选）
total_distance_km <- units::set_units(total_distance, "km")

# 输出结果
cat("驾驶距离（米）：", total_distance, "\n")
cat("驾驶距离（公里）：", total_distance_km, "\n")

关键说明 #

• directed = FALSE：将道路网设置为无向，符合实际驾驶双向通行场景
• convert(to_spatial_smooth)：修复路网可能存在的拓扑错误（如节点不连通、线段断裂）
• weights = "length"：基于道路实际长度计算最短路径，与几何长度一致
• 通过slice和edge_path提取路径对应的道路段，确保获取正确的几何对象用于长度计算

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