在R语言中计算折线沿线两点间的距离

嘿，这问题我熟！既然你已经有了蜿蜒公路的折线shapefile和落在这条线上的经纬度点，用R的sf包就能完美解决沿折线算距离的需求，不用再纠结大圆直线距离啦。下面是一步步的实操指南：

步骤1：加载必备工具包 #

首先得把处理空间数据的核心包装上（如果没装的话先跑install.packages(c("sf", "dplyr"))）：

library(sf)
library(dplyr)

步骤2：读取并预处理空间数据 #

读取折线shapefile #

先把你的公路shapefile读进来，确保坐标系是WGS84（也就是EPSG:4326，对应你的十进制度数经纬度），如果不是的话转一下：

# 读取shapefile
road <- st_read("你的公路文件路径.shp")

# 转成WGS84坐标系（如果原数据不是的话）
road <- st_transform(road, crs = 4326)

# 如果是多段线（MULTILINESTRING），合并成单条折线（LINESTRING）
road_line <- road %>% 
  st_combine() %>% 
  st_cast("LINESTRING")

读取并转换点数据 #

把csv里的经纬度点转换成sf格式的空间点对象，同样对齐坐标系：

# 读取csv
points_df <- read.csv("你的点数据文件.csv")

# 转成sf点对象，替换成你实际的经度、纬度列名
points_sf <- st_as_sf(
  points_df,
  coords = c("longitude", "latitude"), # 这里改成你csv里的经纬度列名
  crs = 4326
)

步骤3：转成投影坐标系（关键！） #

因为WGS84是地理坐标系，直接算距离会有误差，我们要转成平面投影坐标系（比如UTM，单位是米），根据你的公路所在区域选对应的EPSG代码（比如美国东部选EPSG:32618，西部选EPSG:32610，不确定的话可以查一下所在区域的UTM带）：

# 替换成你所在区域的EPSG代码
target_crs <- 32618

road_proj <- st_transform(road_line, crs = target_crs)
points_proj <- st_transform(points_sf, crs = target_crs)

步骤4：计算每个点到折线起点的沿线距离 #

用st_project函数就能直接算出每个点沿折线到起点的累计距离，结果单位是米：

# 给点数据添加沿线距离列
points_proj <- points_proj %>% 
  mutate(
    distance_from_start = st_project(road_proj, .)$distance
  )

步骤5：计算任意两点的沿折线通行距离 #

现在只要取两个点的distance_from_start的绝对值差，就是它们沿公路的通行距离啦！比如你要算id为"point_A"和"point_B"的距离：

# 替换成你实际的点ID列名和ID值
point_a_dist <- points_proj$distance_from_start[points_proj$id == "point_A"]
point_b_dist <- points_proj$distance_from_start[points_proj$id == "point_B"]

# 沿折线的通行距离（单位：米）
route_distance <- abs(point_a_dist - point_b_dist)
print(route_distance)

小提示 #

• 如果你需要公里为单位，只要把结果除以1000就行
• 确保所有点确实都在折线上，否则st_project会计算点到折线的最近点距离，不过你说所有点都在折线上，所以没问题

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