基于多边形分割GPX轨迹并计算区域耗时与距离的技术求助

问题描述 #

我有多份.gpx轨迹文件，每份轨迹均跨越多个多边形处理区域，同时具备对应多边形的shapefile文件。已编写R代码实现轨迹分割与耗时、距离计算，但仍存在两个问题：

• 当多边形内仅存在单个轨迹点时，无法统计该点对应的停留时间（即过渡时间），导致数据丢失，需优化代码避免该问题；
• 需批量处理约100份.gpx文件，最终生成包含原文件名、各多边形总耗时与总距离的数据集，用于评估捕获率的投入力度。

现有代码 #

library(sf)         # 空间数据处理
library(lubridate)  # 时间处理
library(dplyr)      # 数据清洗
library(geosphere)  # 距离计算
library(sp)         # 空间操作
library(tidyverse)  # 工具集（可选）
library(trackeR)    # GPX解析（可选）

# 加载轨迹
gpx_points <- sf::read_sf("try1/Track_2023-10-20 FRINGE_cjh.gpx", layer = "track_points")

# 加载多边形
polygons <- st_read("try1/polygons_combo.shp")  

# 统一坐标系
gpx_points <- st_transform(gpx_points, st_crs(polygons))

# 创建轨迹线
gpx_line <- gpx_points %>%
  arrange(time) %>%
  summarise(do_union = FALSE) %>%
  st_cast("LINESTRING")

# 分割轨迹线
segments <- st_intersection(gpx_line, polygons)

# 匹配轨迹点到多边形
gpx_in_poly <- st_join(gpx_points, polygons, join = st_within)

# 计算各多边形的耗时与距离
gpx_in_poly1 <- gpx_in_poly %>%
  arrange(time) %>%
  group_by(SiteName) %>%
  mutate(
    next_time = lead(time),
    next_geom = lead(geometry),
    time_diff = as.numeric(difftime(next_time, time, units = "secs")),
    dist_m = geosphere::distHaversine(
      st_coordinates(geometry),
      st_coordinates(next_geom)
    )
  ) %>%
  # 过滤时间间隔大于60秒的片段
  filter(time_diff <= 60) %>%
  summarise(
    total_time_sec = sum(time_diff, na.rm = TRUE),
    total_distance_m = sum(dist_m, na.rm = TRUE)
  )

gpx_in_poly1
 
# 快速可视化
library(ggplot2)
ggplot() +
  geom_sf(data = polygons, fill = NA, color = "grey40") +
  geom_sf(data = gpx_in_poly, aes(color = factor(SiteName)), size = 0.8) +
  scale_color_viridis_d(na.value = "black") +
  theme_minimal()

优化方案 #

1. 解决单个轨迹点的停留时间统计问题 #

原代码中单个轨迹点因无后续点导致time_diff为NA，无法统计耗时。优化思路是：先标记每个多边形内的轨迹点数量，对单个点的情况，根据其在整条轨迹中的位置（起点、终点、中间点）分配合理的停留时间。

修改后的核心计算代码：

gpx_in_poly1 <- gpx_in_poly %>%
  arrange(time) %>%
  # 全局标记前后点信息，方便跨分组调用
  mutate(
    prev_time = lag(time),
    next_time_global = lead(time)
  ) %>%
  group_by(SiteName) %>%
  mutate(
    point_count = n(),
    # 多点组常规计算时间差
    time_diff = ifelse(point_count > 1, 
                       as.numeric(difftime(lead(time), time, units = "secs")),
                       NA),
    dist_m = ifelse(point_count > 1,
                    geosphere::distHaversine(st_coordinates(geometry), st_coordinates(lead(geometry))),
                    NA),
    # 单点组按轨迹位置分配时间
    single_point_time = ifelse(point_count == 1,
                               case_when(
                                 # 轨迹起点，取与下一个点的时间差
                                 is.na(prev_time) ~ as.numeric(difftime(next_time_global, time, units = "secs")),
                                 # 轨迹终点，取与上一个点的时间差
                                 is.na(next_time_global) ~ as.numeric(difftime(time, prev_time, units = "secs")),
                                 # 中间单点，取前后点时间差的平均值
                                 TRUE ~ as.numeric(difftime(next_time_global, prev_time, units = "secs")) / 2
                               ),
                               NA)
  ) %>%
  filter(is.na(time_diff) | time_diff <= 60) %>%
  summarise(
    total_time_sec = sum(time_diff, na.rm = TRUE) + sum(single_point_time, na.rm = TRUE),
    total_distance_m = sum(dist_m, na.rm = TRUE)
  )

2. 批量处理.gpx文件 #

通过定义单个文件处理函数，结合purrr工具批量遍历所有.gpx文件，合并结果并保留原文件名：

# 定义单个gpx文件的处理函数
process_gpx <- function(gpx_path, polygons) {
  # 读取并预处理轨迹点
  gpx_points <- sf::read_sf(gpx_path, layer = "track_points") %>%
    st_transform(st_crs(polygons)) %>%
    arrange(time)
  
  # 匹配轨迹点到多边形
  gpx_in_poly <- st_join(gpx_points, polygons, join = st_within)
  
  # 计算耗时与距离（包含单点处理）
  result <- gpx_in_poly %>%
    mutate(
      prev_time = lag(time),
      next_time_global = lead(time)
    ) %>%
    group_by(SiteName) %>%
    mutate(
      point_count = n(),
      time_diff = ifelse(point_count > 1, 
                         as.numeric(difftime(lead(time), time, units = "secs")),
                         NA),
      dist_m = ifelse(point_count > 1,
                      geosphere::distHaversine(st_coordinates(geometry), st_coordinates(lead(geometry))),
                      NA),
      single_point_time = ifelse(point_count == 1,
                                 case_when(
                                   is.na(prev_time) ~ as.numeric(difftime(next_time_global, time, units = "secs")),
                                   is.na(next_time_global) ~ as.numeric(difftime(time, prev_time, units = "secs")),
                                   TRUE ~ as.numeric(difftime(next_time_global, prev_time, units = "secs")) / 2
                                 ),
                                 NA)
    ) %>%
    filter(is.na(time_diff) | time_diff <= 60) %>%
    summarise(
      total_time_sec = sum(time_diff, na.rm = TRUE) + sum(single_point_time, na.rm = TRUE),
      total_distance_m = sum(dist_m, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    # 添加原文件名
    mutate(gpx_filename = basename(gpx_path))
  
  return(result)
}

# 批量处理流程
# 获取目标文件夹下所有gpx文件路径
gpx_files <- list.files(path = "try1", pattern = "\\.gpx$", full.names = TRUE)

# 加载多边形数据并统一坐标系
polygons <- st_read("try1/polygons_combo.shp") %>%
  st_transform(4326) # 适配GPX默认的WGS84坐标系

# 批量处理并合并所有结果
final_dataset <- purrr::map_dfr(gpx_files, ~process_gpx(.x, polygons))

# 查看结果
head(final_dataset)

# 保存结果为CSV文件
write.csv(final_dataset, "gpx_polygon_stats.csv", row.names = FALSE)

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