r - 在 R 中按顺序将点捕捉到线

Question

我有一组 GPS 点和一个线串(代表一条公交线路)，GPS 点应该属于(两者都是有序的)。所以我使用了一个函数来捕捉线串的点:

library(dplyr)
library(sf)
library(readr)

# Function to snap points to the closest line

snap_points_to_line <- function(points, line) {

  # alinhar as pradas gps com a linha
  points_align <- st_nearest_points(points, line) %>%
    st_cast("POINT")

  # pegar so os pontos pares
  points_new_geometry <- points_align[c(seq(2, length(points_align), by = 2))]

  points_align_end <- points %>%
    st_set_geometry(points_new_geometry)

}

# GPS Points
gps <- structure(list(id = 1:11, 
                      lon = c(-38.477035, -38.477143, -38.478701, 
                              -38.479795, -38.480923, -38.481078, 
                              -38.481885, -38.484545, -38.486156, 
                              -38.486813, -38.486506), 
                      lat = c(-3.743078, -3.743019, -3.742566, 
                              -3.742246, -3.741844, -3.741853, 
                              -3.741596, -3.740711, -3.740076, 
                              -3.739399, -3.73886)), 
                 class = "data.frame", 
                 row.names = c(NA,-11L))

gps

   id       lon       lat
1   1 -38.47704 -3.743078
2   2 -38.47714 -3.743019
3   3 -38.47870 -3.742566
4   4 -38.47980 -3.742246
5   5 -38.48092 -3.741844
6   6 -38.48108 -3.741853
7   7 -38.48188 -3.741596
8   8 -38.48454 -3.740711
9   9 -38.48616 -3.740076
10 10 -38.48681 -3.739399
11 11 -38.48651 -3.738860

# Download line
line <- read_rds(gzcon(url("https://github.com/kauebraga/dissertacao/raw/master/junk/line_so.rds")))

# Make snap
gps_snap <- snap_points_to_line(gps %>% st_as_sf(coords = c("lon", "lat"), crs = 4326), line)

大多数情况下，快照工作正常。但也有一些情况下公交线路会掉头和 一些点被捕捉到了道路的错误一侧 因为 GPS 位置可能有误差。在下图中，道路南侧的三个点应该在北侧:



如何保证 GPS 点捕捉到正确的道路一侧？ GPS 点和线串的顺序相同(如果您 st_cast(line, "POINT) 它将给出与 GPS 一起增长的点)，所以我希望应该有办法解决这个问题(我不知道如何!)。

使用 sf 的任何帮助或 R 中的其他空间工具将不胜感激。谢谢!

Answer 1

设置数据

library(sf)
library(dplyr)
library(readr)
library(rgeos)

# GPS Points
gps <- structure(list(id = 1:11, 
                      lon = c(-38.477035, -38.477143, -38.478701, 
                              -38.479795, -38.480923, -38.481078, 
                              -38.481885, -38.484545, -38.486156, 
                              -38.486813, -38.486506), 
                      lat = c(-3.743078, -3.743019, -3.742566, 
                              -3.742246, -3.741844, -3.741853, 
                              -3.741596, -3.740711, -3.740076, 
                              -3.739399, -3.73886)), 
                 class = "data.frame", 
                 row.names = c(NA,-11L))

# convert to sf
gps <- gps %>% st_as_sf(coords = c("lon", "lat"), crs = 4326, remove =F) 

line <- read_rds(gzcon(url("https://github.com/kauebraga/dissertacao/raw/master/junk/line_so.rds"))) 

定义自定义捕捉函数
此函数适用于以下逻辑:要捕捉到的正确路段是需要从前一点沿线串(网络距离)行驶的最短距离的路段。
它执行以下操作:

每个点由给定的 tolerance 缓冲(以米为单位，因此我们已转换为您所在地区的米 CRS)

然后该线与我们的缓冲区相交。这将留下两段交通往来的道路，另一段则不然。这如下图所示:

在某些情况下，我们现在有两个选项可以捕捉，所以我们最初捕捉到两者:

我们选择一个明确的点(只有一个捕捉选项)作为引用点，并计算沿网络到下一个 id 的捕捉选项的距离。

对于每个点 id，与前一个 id 网络距离最小的点就是我们想要的点。

找到正确的点 id 后，我们将其设置为新的引用点并从步骤 4 开始重复。

custom_snap <- function(line, points, tolerance, crs = 29194) {
  points <- st_transform(points, crs)
  line <- st_transform(line, crs)
  # buffer the points by the tolerance
  points_buf <- st_buffer(points, 15)
  # intersect the line with the buffer
  line_intersect <- st_intersection(line, points_buf)
  # convert mutlinestrings (more than one road segment) into linestrings
  line_intersect <- do.call(rbind,lapply(1:nrow(line_intersect),function(x){st_cast(line_intersect[x,],"LINESTRING")}))

  # for each line intersection, calculate the nearest point on that line to our gps point
  nearest_pt <- do.call(rbind,lapply(seq_along(points$id), function(i){
    points[points$id==i,] %>%  st_nearest_points(line_intersect[line_intersect$id==i,]) %>% st_sf %>%
      st_cast('POINT') %>% mutate(id = i)
    }))

  nearest_pt<- nearest_pt[seq(2, nrow(nearest_pt), by = 2),] %>%
    mutate(option = 1:nrow(.))

  # find an unambiguous reference point with only one snap option
  unambiguous_pt <- nearest_pt %>%
    group_by(id) %>%
    mutate(count = n()) %>%
    ungroup() %>%
    filter(count == 1) %>%
    slice(1)

  # calculate network distance along our line to each snapped point
  dists <- rgeos::gProject(as(line,'Spatial'), as(nearest_pt,'Spatial'))
  # join back to nearest points data
  dists <- nearest_pt %>% cbind(dists)

  # we want to recursively do the following:
  # 1. calculate the network distance from our unambiguous reference point to the next id point in the data
  # 2. keep the snapped point for that id that was closest *along the network*  to the previous id
  # 3. set the newly snapped point as our reference point
  # 4. repeat

  # get distances from our reference point to the next point id
  for(i in unambiguous_pt$id:(max(dists$id)-1)){
    next_dist <- which.min(abs(dists[dists$id== i +1,]$dists - dists[dists$id== unambiguous_pt$id,]$dists ))
    next_option <- dists[dists$id== i +1,][next_dist,]$option
    nearest_pt <- nearest_pt %>% filter(id != i+1 | option == next_option)
    unambiguous_pt <- nearest_pt %>% filter(id ==i+1 & option == next_option)
    dists <- nearest_pt %>% cbind(dists = rgeos::gProject(as(line,'Spatial'), as(nearest_pt,'Spatial')))
  }

  # and in the reverse direction
  for(i in unambiguous_pt$id:(min(dists$id)+1)){
    next_dist <- which.min(abs(dists[dists$id== i -1,]$dists - dists[dists$id== unambiguous_pt$id,]$dists ))
    next_option <- dists[dists$id== i -1,][next_dist,]$option
    nearest_pt <- nearest_pt %>% filter(id != i-1 | option == next_option)
    unambiguous_pt <- nearest_pt %>% filter(id ==i-1 & option == next_option)
    dists <- nearest_pt %>% cbind(dists = rgeos::gProject(as(line,'Spatial'), as(nearest_pt,'Spatial')))
  }

  # transform back into lat/lng
  snapped_points <- nearest_pt %>%
    st_transform(4326)

  return(snapped_points)
}

计算要捕捉到哪条线

gps_snap <- custom_snap(line, gps, 15) %>%
  cbind(st_coordinates(.))

在传单中绘制结果

library(leaflet)
# get line coords
line_coords <- line %>%
  st_coordinates(.) 

# plot in leaflet
leaflet() %>%
  leaflet::setView(lng = -38.4798, lat = -3.741829, zoom = 18) %>%
  addProviderTiles('CartoDB.Positron') %>%
  addPolylines(lng = line_coords[,'X'], lat = line_coords[,'Y']) %>%
  addCircles(lng = gps$lon, lat = gps$lat, radius = 3, color ='red') %>%
  addCircles(lng = gps_snap$X, lat = gps_snap$Y, col ='green', radius = 4) 

绿色代表捕捉点，红色代表原始 GPS 点。它们现在被捕捉到道路的正确一侧。