使用gtfstools过滤ATM的GTFS空间范围时遇sfheaders错误求助

问题背景 #

合并米兰地区多家运营商GTFS数据时，使用filter_by_spatial_extent函数做空间过滤，仅保留米兰区域内服务。其余运营商（包括无shapes.txt的文件）均能正常过滤，但市政运营商ATM的GTFS触发以下错误：

sfheaders - error indexing lines, perhaps caused by un-ordered data?

已提前按shape_id和shape_pt_sequence对ATM的shapes数据排序，且验证过数据结构有效性，但问题仍未解决。

排查与解决步骤 #

1. 检查单顶点shape_id #

sfheaders处理仅含单个坐标点的shape时会报错，先筛选并移除这类数据：

# 提取ATM的shapes数据
atm_shapes <- gtfs_list[["atm"]]$shapes

# 统计每个shape_id的点数
shape_point_count <- atm_shapes[, .N, by = shape_id]

# 找出仅含1个点的shape_id
single_point_shapes <- shape_point_count[N == 1]$shape_id

if(length(single_point_shapes) > 0) {
  message("发现 ", length(single_point_shapes), " 个单顶点shape")
  # 移除这些异常shape
  gtfs_list[["atm"]]$shapes <- atm_shapes[!shape_id %in% single_point_shapes]
}

处理后重新执行filter_by_spatial_extent测试。

2. 验证坐标格式与空间参考一致性 #

GTFS标准要求坐标为WGS84（EPSG:4326），需确认ATM数据与边界数据的空间参考匹配：

# 检查stibm的坐标系
st_crs(stibm)

# 检查ATM shapes的坐标范围（WGS84经度范围：-180~180，纬度：-90~90）
range(atm_shapes$shape_pt_lon)
range(atm_shapes$shape_pt_lat)

若坐标范围异常，需将数据转换为EPSG:4326后再执行过滤。

3. 手动测试sfheaders转换，定位异常shape #

直接用ATM的shapes数据构建sf线对象，定位触发错误的具体shape：

# 批量测试sfheaders转换
tryCatch({
  sf_lines <- sfheaders::sf_linestring(
    atm_shapes,
    x = "shape_pt_lon",
    y = "shape_pt_lat",
    linestring_id = "shape_id"
  )
}, error = function(e) {
  message("sfheaders直接转换报错：", e$message)
})

# 逐个shape_id排查异常
test_shapes <- unique(atm_shapes$shape_id)
for(sid in test_shapes) {
  sub_shapes <- atm_shapes[shape_id == sid]
  try({
    sfheaders::sf_linestring(sub_shapes, x = "shape_pt_lon", y = "shape_pt_lat")
  }, silent = FALSE)
}

找到异常shape_id后，检查其点序列是否存在重复序号、缺失坐标等问题，修复后再执行过滤。

4. 备选：绕过shapes，通过站点过滤 #

若以上方法无效，可先过滤区域内的站点，再关联保留对应的线路与班次：

# 转换ATM站点为sf对象
atm_stops <- gtfs_list[["atm"]]$stops %>%
  sfheaders::sf_point(x = "stop_lon", y = "stop_lat") %>%
  st_set_crs(4326)

# 筛选米兰区域内的站点
valid_stops <- atm_stops[st_intersects(atm_stops, stibm, sparse = FALSE), ]$stop_id

# 保留关联的班次、线路与shape
filtered_atm <- gtfs_list[["atm"]]
filtered_atm$stop_times <- filtered_atm$stop_times[stop_id %in% valid_stops]
filtered_atm$trips <- filtered_atm$trips[trip_id %in% filtered_atm$stop_times$trip_id]
filtered_atm$routes <- filtered_atm$routes[route_id %in% filtered_atm$trips$route_id]
filtered_atm$shapes <- filtered_atm$shapes[shape_id %in% filtered_atm$trips$shape_id]

内容的提问来源于stack exchange，提问作者Antonio de Rosa