在 R 中按顺序对齐点
Snap points to line in order in R
我有一组 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 中其他空间工具的帮助将不胜感激。谢谢！
设置数据
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点。它们现在被捕捉到道路的正确一侧。
我有一组 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 中其他空间工具的帮助将不胜感激。谢谢！
设置数据
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点。它们现在被捕捉到道路的正确一侧。