如何用sf或R其他包合并GPS轨迹估算真实道路位置？

我有多条同一道路路段的GPS通行轨迹数据，希望合并这些轨迹以得到真实道路位置的最优估算，能否用sf包（或其他R包）实现？

以下是可复现代码示例：

所需R包 #

library(sf)
library(dplyr)
library(magrittr)
library(rnaturalearthdata)
library(ggplot2) # 后续可视化用

生成「真实道路」样本 #

从coastline数据集中提取一段线段作为真实道路：

road_real_coo <-  st_coordinates(coastline50[16,1])[-c(1:31, 67:71),1:2] 
road_real <-  road_real_coo %>% 
    as.data.frame() %>% 
    st_as_sf(coords=c("X", "Y"), crs=4326) %>% 
    group_by() %>% 
    summarise(do_union=FALSE) %>% 
    st_cast("LINESTRING") 

真实道路可视化：

模拟GPS通行轨迹 #

从真实道路顶点采样并添加随机扰动模拟GPS误差，生成4条轨迹：

set.seed(123)    
road_gps_1 <- road_real_coo[sort(sample(1:nrow(road_real_coo), 10)),] %>% 
  magrittr::add(runif(20,-.05,.05)) %>% 
  as.data.frame() %>% 
  st_as_sf(coords=c("X", "Y"), crs=4326) %>% 
  group_by() %>% 
  summarise(do_union=FALSE) %>% 
  st_cast("LINESTRING")
road_gps_2 <- road_real_coo[sort(sample(1:nrow(road_real_coo), 12)),] %>% 
  magrittr::add(runif(24,-.05,.05)) %>% 
  as.data.frame() %>% 
  st_as_sf(coords=c("X", "Y"), crs=4326) %>% 
  group_by() %>% 
  summarise(do_union=FALSE) %>% 
  st_cast("LINESTRING")
road_gps_3 <- road_real_coo[sort(sample(1:nrow(road_real_coo), 8)),] %>% 
  magrittr::add(runif(16,-.05,.05)) %>% 
  as.data.frame() %>% 
  st_as_sf(coords=c("X", "Y"), crs=4326) %>% 
  group_by() %>% 
  summarise(do_union=FALSE) %>% 
  st_cast("LINESTRING")
road_gps_4 <- road_real_coo[sort(sample(1:nrow(road_real_coo), 16)),] %>% 
  magrittr::add(runif(32,-.05,.05)) %>% 
  as.data.frame() %>% 
  st_as_sf(coords=c("X", "Y"), crs=4326) %>% 
  group_by() %>% 
  summarise(do_union=FALSE) %>% 
  st_cast("LINESTRING")
road_gps <- rbind(road_gps_1,road_gps_2,road_gps_3,road_gps_4)

GPS轨迹与真实道路对比可视化：

plot(road_real, reset=F, lwd=2)
plot(road_gps, add=T, col = 2:5)

目标是合并这些GPS轨迹，得到类似手动拟合的估算道路：

补充场景细节 #

• 需处理约20万-30万条路段；
• 每个路段的通行轨迹数量为2至2000条不等；
• 每条轨迹的行驶方向随机；
• 每条轨迹包含2至100个GPS点；
• 少数路段存在自重叠情况；
• 不匹配公开道路数据，仅用现有数据估算。

解决方案 #

可以通过轨迹方向统一、点插值对齐、均值拟合的流程实现，结合sf包和其他空间工具包完成，以下是具体步骤：

步骤1：统一所有轨迹的行驶方向 #

由于轨迹方向随机，先将所有轨迹对齐到同一方向（以真实道路的起点-终点为基准）：

# 获取真实道路的起点和终点坐标
real_start <- st_coordinates(road_real)[1,]
real_end <- st_coordinates(road_real)[nrow(st_coordinates(road_real)),]

# 定义反转轨迹方向的函数
reverse_line <- function(line) {
  coords <- st_coordinates(line)
  reversed_coords <- coords[nrow(coords):1,]
  reversed_line <- reversed_coords %>%
    as.data.frame() %>%
    st_as_sf(coords=c("X", "Y"), crs=st_crs(line)) %>%
    summarise(do_union=FALSE) %>%
    st_cast("LINESTRING")
  return(reversed_line)
}

# 统一每条GPS轨迹的方向
road_gps_aligned <- road_gps %>%
  rowwise() %>%
  mutate(
    line_start = st_coordinates(geometry)[1,],
    dist_to_real_start = sqrt((line_start[1]-real_start[1])^2 + (line_start[2]-real_start[2])^2),
    dist_to_real_end = sqrt((line_start[1]-real_end[1])^2 + (line_start[2]-real_end[2])^2),
    geometry = ifelse(dist_to_real_start > dist_to_real_end, reverse_line(geometry), geometry)
  ) %>%
  st_sf() %>%
  select(-line_start, -dist_to_real_start, -dist_to_real_end)

步骤2：对轨迹进行均匀采样对齐 #

不同轨迹点数不同，用st_line_sample对每条轨迹进行均匀采样，统一点数：

# 对每条轨迹采样20个均匀分布的点
sample_n <- 20
gps_points <- road_gps_aligned %>%
  rowwise() %>%
  mutate(sampled = list(st_line_sample(geometry, n=sample_n) %>% st_cast("POINT"))) %>%
  st_drop_geometry() %>%
  tidyr::unnest(sampled) %>%
  st_as_sf()

# 给每个采样点分配位置索引（第1个点、第2个点...）
gps_points <- gps_points %>%
  group_by(rowid = rep(1:nrow(road_gps_aligned), each=sample_n)) %>%
  mutate(point_idx = 1:n()) %>%
  ungroup()

步骤3：计算均值坐标生成估算道路 #

按位置索引对所有点的经纬度取均值，再连接成线：

# 计算每个位置索引的平均坐标
estimated_road_points <- gps_points %>%
  group_by(point_idx) %>%
  summarise(
    X = mean(st_coordinates(sampled)[,1]),
    Y = mean(st_coordinates(sampled)[,2])
  ) %>%
  st_as_sf(coords=c("X", "Y"), crs=4326)

# 将均值点连接成线
estimated_road <- estimated_road_points %>%
  summarise(do_union=FALSE) %>%
  st_cast("LINESTRING")

步骤4：可视化验证 #

ggplot() +
  geom_sf(data=road_real, color="black", size=1.5) +
  geom_sf(data=road_gps, color=rgb(0.8,0.2,0.2,0.5), size=0.8) +
  geom_sf(data=estimated_road, color="blue", size=1.2, linetype="dashed") +
  theme_minimal()

大规模数据优化建议 #

针对20-30万条路段的场景，需提升处理效率：

• 用data.table替代dplyr进行分组计算，处理大数据量速度更快；
• 用purrr的向量化函数批量处理单路段轨迹，避免循环；
• 对存在自重叠的路段，先用dbscan包做空间聚类，分离重叠区域后再分别拟合；
• 轨迹数量少的路段（2-5条）直接取几何均值线；轨迹数量多的路段（>100条）用ks包做核密度估计，提取高密度中心线替代简单均值。

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