将 GPS 点插值到 R 中的其他点集
Interpolate GPS points to other set of points in R
我有两个数据集：一个 GPS 点（来自公共汽车位置）的时间戳大多每 30 秒收集一次，另一个是落在 GPS 轨迹内的中间点（来自公共汽车站）。
# GPS Points
gps_points <- structure(list(id_gps = c(4138176L, 4136334L, 4134534L, 4132685L,
4130891L, 4129035L, 4127232L, 4125387L, 4123620L, 4121861L, 4120114L,
4118381L, 4116721L, 3380373L, 3374532L, 3369036L, 3363258L, 3357540L,
3351543L, 3345549L, 3339777L, 3333210L, 3326793L, 3319251L, 3312822L,
3306501L), hora = structure(c(1535953786, 1535953816, 1535953846,
1535953876, 1535953906, 1535953936, 1535953966, 1535953996, 1535954026,
1535954056, 1535954086, 1535954116, 1535954146, 1535954176, 1535954206,
1535954236, 1535954266, 1535954296, 1535954326, 1535954356, 1535954386,
1535954416, 1535954446, 1535954476, 1535954506, 1535954536), class = c("POSIXct",
"POSIXt"), tzone = "UTC"), lon = c(-38.500763, -38.501413, -38.50252,
-38.503505, -38.504694, -38.505441, -38.506651, -38.507328, -38.507965,
-38.509063, -38.509735, -38.51022, -38.511546, -38.511778, -38.512788,
-38.513633, -38.514568, -38.51495, -38.515331, -38.515878, -38.516438,
-38.516628, -38.517129, -38.517651, -38.518056, -38.518358),
lat = c(-3.80892, -3.807633, -3.805113, -3.802854, -3.800343,
-3.79881, -3.796178, -3.79474, -3.793426, -3.791048, -3.789561,
-3.78856, -3.78569, -3.785216, -3.783108, -3.781245, -3.778751,
-3.777118, -3.775673, -3.773774, -3.771845, -3.771159, -3.769336,
-3.767198, -3.765478, -3.764019)), row.names = c(NA, -26L
), class = "data.frame", .Names = c("id_gps", "hora", "lon",
"lat"))
# Stops
stops <- structure(list(stop_id = c(4873, 3215, 5083, 3346, 3363, 3362,
3542, 3543, 3540, 4629, 3528), lon = c(-38.516766, -38.515311,
-38.513903, -38.512154, -38.511001, -38.509844, -38.508943, -38.50816,
-38.507062, -38.505798, -38.504044), lat = c(-3.771828, -3.77695,
-3.781432, -3.785157, -3.787631, -3.790069, -3.791997, -3.793663,
-3.796027, -3.798711, -3.802504)), class = c("tbl_df", "tbl",
"data.frame"), row.names = c(NA, -11L), .Names = c("stop_id",
"lon", "lat"))
head(gps_points)
id_gps hora lon lat
1 4138176 2018-09-03 05:49:46 -38.50076 -3.808920
2 4136334 2018-09-03 05:50:16 -38.50141 -3.807633
3 4134534 2018-09-03 05:50:46 -38.50252 -3.805113
4 4132685 2018-09-03 05:51:16 -38.50350 -3.802854
5 4130891 2018-09-03 05:51:46 -38.50469 -3.800343
6 4129035 2018-09-03 05:52:16 -38.50544 -3.798810
head(stops)
stop_id lon lat
1 4873 -38.51677 -3.771828
2 3215 -38.51531 -3.776950
3 5083 -38.51390 -3.781432
4 3346 -38.51215 -3.785157
5 3363 -38.51100 -3.787631
6 3362 -38.50984 -3.790069
GPS 点为红色，停止点为蓝色
我想估计（通过 GPS 数据集上的线性插值）使用 R 估计与每个停靠点相关联的时间戳。 所需的输出将是带有新列的停靠点数据集识别内插时间戳。
我现在正在做的解决方法包括每 5 秒将 GPS 点插入事件（使用 this method），然后计算每个站点最近的 GPS 事件（使用 RANN::nn2）。 GPS 数据集非常大，所以这在计算上是不可行的，而且我仍然没有得到与每个站点关联的 "exact" 时间戳。
# workaround
# Crete combination of timestamps for each 5 seconds
full.time <- with(gps_points,seq(gps_points$hora[1],tail(gps_points$hora,1),by=5))
library(zoo)
# convert to zoo object
df.zoo <- zoo(gps_points[,c("lon", "lat")],gps_points$hora)
# interpolate; result is also a zoo object
result <- na.approx(df.zoo,xout=full.time)
# transform zoo to df
zoo.to.data.frame <- function(x, index.name="hora") {
stopifnot(is.zoo(x))
xn <- if(is.null(dim(x))) deparse(substitute(x)) else colnames(x)
setNames(data.frame(index(x), x, row.names=NULL), c(index.name,xn))
}
gps_points_interpolated <- zoo.to.data.frame(result) %>% as_tibble()
# Create temp_id for stops
stops <- stops %>%
mutate(temp_id = 1:n())
# To each GPS point, what's the closest stop?
opa <- RANN::nn2(select(stops, lon, lat), select(gps_points_interpolated, lon, lat), 1)
vamos <- gps_points_interpolated %>%
mutate(temp_id = opa$nn.idx, dist = opa$nn.dists*111320)
# Bring back stop_id, lon e lat of each stop
vamos <- left_join(vamos, stops, by = "temp_id", suffix = c(".gps", ".stop")) %>%
# Select columns
select(stop_id, lon = lon.stop, lat = lat.stop, hora, dist)
# Select the observations that have minimun distance to each stop
vamos_fim <- vamos %>%
group_by(stop_id) %>%
slice(which.min(dist))
head(vamos_fim)
# A tibble: 6 x 5
# Groups: stop_id [6]
stop_id lon lat hora dist[,1]
<dbl> <dbl> <dbl> <dttm> <dbl>
1 3215 -38.5 -3.78 2018-09-03 05:58:21 34.1
2 3346 -38.5 -3.79 2018-09-03 05:56:21 39.9
3 3362 -38.5 -3.79 2018-09-03 05:54:36 37.1
4 3363 -38.5 -3.79 2018-09-03 05:55:26 37.9
5 3528 -38.5 -3.80 2018-09-03 05:51:21 38.7
6 3540 -38.5 -3.80 2018-09-03 05:52:51 34.6
也欢迎使用线性插值以外的方法。谢谢。
我设法使用 zoo::na.approx 函数自行解决了这个问题。首先，您需要计算两个数据集中从同一点开始的连续点之间的累积距离。在我的示例中，我的 GPS 点按顺序排在第一位，因此我将选择第一个 GPS 点作为两个数据集的 "starting point"。
# I just realized that the datasets are growing in opposite direction, so I'll flip the stops
stops <- map_df(stops, rev)
# Function to calculate distance from previous point
get.dist <- function(lon, lat) geosphere::distHaversine(tail(cbind(lon,lat),-1),head(cbind(lon,lat),-1))
# Calculate cumulative distance of gps points (points must be ordered by time)
gps_points <- gps_points %>%
mutate(dist = c(0, cumsum(get.dist(lon, lat))))
# Input first GPS point and calculate cumulative distance of stops (also must be ordered)
stops <- gps_points %>%
# Select only the first point
slice(1) %>%
# Select columns to match the stops dataset
mutate(stop_id = NA) %>%
select(stop_id, lon, lat) %>%
# Input the stop points
rbind(stops) %>%
# Calculate cumulative dist
mutate(dist = c(0, cumsum(get.dist(lon, lat))))
# Interpolate ------------------------------
x <- gps_points$hora
y <- gps_points$dist
# to which position we want to interpolate? to the stops!
xout <- stops$dist
interp <- as.POSIXct(zoo::na.approx(x ,
y,
xout = xout,
ties = "ordered",
rule = 2),
origin = "1970-01-01")
# Put it together
stops_interp <- stops %>%
# Input the interpolated times
mutate(hora = interp) %>%
# Delete the first row that was inputed from the GPS
slice(-1)
stops_interp
stop_id lon lat dist hora
1 3528 -38.50404 -3.802504 801.8295 2018-09-03 05:51:21
2 4629 -38.50580 -3.798711 1266.8446 2018-09-03 05:52:18
3 3540 -38.50706 -3.796027 1596.9692 2018-09-03 05:52:51
4 3543 -38.50816 -3.793663 1887.0162 2018-09-03 05:53:43
5 3542 -38.50894 -3.791997 2091.8550 2018-09-03 05:54:05
6 3362 -38.50984 -3.790069 2328.6658 2018-09-03 05:54:38
7 3363 -38.51100 -3.787631 2628.9530 2018-09-03 05:55:27
8 3346 -38.51215 -3.785157 2932.6796 2018-09-03 05:56:18
9 5083 -38.51390 -3.781432 3390.5980 2018-09-03 05:57:15
10 3215 -38.51531 -3.776950 3913.4700 2018-09-03 05:58:22
11 4873 -38.51677 -3.771828 4506.1116 2018-09-03 05:59:52
第hora列表示插值时间。
更进一步
但是，这种插值是基于点之间的欧氏距离。如果我有 GPS 点和停靠点的预期路径（如道路）怎么办？在这种情况下，我建议使用 rgeos::gProject 函数来计算路径中 GPS 和停靠点的累积距离。
我有两个数据集：一个 GPS 点（来自公共汽车位置）的时间戳大多每 30 秒收集一次，另一个是落在 GPS 轨迹内的中间点（来自公共汽车站）。
# GPS Points
gps_points <- structure(list(id_gps = c(4138176L, 4136334L, 4134534L, 4132685L,
4130891L, 4129035L, 4127232L, 4125387L, 4123620L, 4121861L, 4120114L,
4118381L, 4116721L, 3380373L, 3374532L, 3369036L, 3363258L, 3357540L,
3351543L, 3345549L, 3339777L, 3333210L, 3326793L, 3319251L, 3312822L,
3306501L), hora = structure(c(1535953786, 1535953816, 1535953846,
1535953876, 1535953906, 1535953936, 1535953966, 1535953996, 1535954026,
1535954056, 1535954086, 1535954116, 1535954146, 1535954176, 1535954206,
1535954236, 1535954266, 1535954296, 1535954326, 1535954356, 1535954386,
1535954416, 1535954446, 1535954476, 1535954506, 1535954536), class = c("POSIXct",
"POSIXt"), tzone = "UTC"), lon = c(-38.500763, -38.501413, -38.50252,
-38.503505, -38.504694, -38.505441, -38.506651, -38.507328, -38.507965,
-38.509063, -38.509735, -38.51022, -38.511546, -38.511778, -38.512788,
-38.513633, -38.514568, -38.51495, -38.515331, -38.515878, -38.516438,
-38.516628, -38.517129, -38.517651, -38.518056, -38.518358),
lat = c(-3.80892, -3.807633, -3.805113, -3.802854, -3.800343,
-3.79881, -3.796178, -3.79474, -3.793426, -3.791048, -3.789561,
-3.78856, -3.78569, -3.785216, -3.783108, -3.781245, -3.778751,
-3.777118, -3.775673, -3.773774, -3.771845, -3.771159, -3.769336,
-3.767198, -3.765478, -3.764019)), row.names = c(NA, -26L
), class = "data.frame", .Names = c("id_gps", "hora", "lon",
"lat"))
# Stops
stops <- structure(list(stop_id = c(4873, 3215, 5083, 3346, 3363, 3362,
3542, 3543, 3540, 4629, 3528), lon = c(-38.516766, -38.515311,
-38.513903, -38.512154, -38.511001, -38.509844, -38.508943, -38.50816,
-38.507062, -38.505798, -38.504044), lat = c(-3.771828, -3.77695,
-3.781432, -3.785157, -3.787631, -3.790069, -3.791997, -3.793663,
-3.796027, -3.798711, -3.802504)), class = c("tbl_df", "tbl",
"data.frame"), row.names = c(NA, -11L), .Names = c("stop_id",
"lon", "lat"))
head(gps_points)
id_gps hora lon lat
1 4138176 2018-09-03 05:49:46 -38.50076 -3.808920
2 4136334 2018-09-03 05:50:16 -38.50141 -3.807633
3 4134534 2018-09-03 05:50:46 -38.50252 -3.805113
4 4132685 2018-09-03 05:51:16 -38.50350 -3.802854
5 4130891 2018-09-03 05:51:46 -38.50469 -3.800343
6 4129035 2018-09-03 05:52:16 -38.50544 -3.798810
head(stops)
stop_id lon lat
1 4873 -38.51677 -3.771828
2 3215 -38.51531 -3.776950
3 5083 -38.51390 -3.781432
4 3346 -38.51215 -3.785157
5 3363 -38.51100 -3.787631
6 3362 -38.50984 -3.790069
GPS 点为红色，停止点为蓝色
我想估计（通过 GPS 数据集上的线性插值）使用 R 估计与每个停靠点相关联的时间戳。 所需的输出将是带有新列的停靠点数据集识别内插时间戳。
我现在正在做的解决方法包括每 5 秒将 GPS 点插入事件（使用 this method），然后计算每个站点最近的 GPS 事件（使用 RANN::nn2）。 GPS 数据集非常大，所以这在计算上是不可行的，而且我仍然没有得到与每个站点关联的 "exact" 时间戳。
# workaround
# Crete combination of timestamps for each 5 seconds
full.time <- with(gps_points,seq(gps_points$hora[1],tail(gps_points$hora,1),by=5))
library(zoo)
# convert to zoo object
df.zoo <- zoo(gps_points[,c("lon", "lat")],gps_points$hora)
# interpolate; result is also a zoo object
result <- na.approx(df.zoo,xout=full.time)
# transform zoo to df
zoo.to.data.frame <- function(x, index.name="hora") {
stopifnot(is.zoo(x))
xn <- if(is.null(dim(x))) deparse(substitute(x)) else colnames(x)
setNames(data.frame(index(x), x, row.names=NULL), c(index.name,xn))
}
gps_points_interpolated <- zoo.to.data.frame(result) %>% as_tibble()
# Create temp_id for stops
stops <- stops %>%
mutate(temp_id = 1:n())
# To each GPS point, what's the closest stop?
opa <- RANN::nn2(select(stops, lon, lat), select(gps_points_interpolated, lon, lat), 1)
vamos <- gps_points_interpolated %>%
mutate(temp_id = opa$nn.idx, dist = opa$nn.dists*111320)
# Bring back stop_id, lon e lat of each stop
vamos <- left_join(vamos, stops, by = "temp_id", suffix = c(".gps", ".stop")) %>%
# Select columns
select(stop_id, lon = lon.stop, lat = lat.stop, hora, dist)
# Select the observations that have minimun distance to each stop
vamos_fim <- vamos %>%
group_by(stop_id) %>%
slice(which.min(dist))
head(vamos_fim)
# A tibble: 6 x 5
# Groups: stop_id [6]
stop_id lon lat hora dist[,1]
<dbl> <dbl> <dbl> <dttm> <dbl>
1 3215 -38.5 -3.78 2018-09-03 05:58:21 34.1
2 3346 -38.5 -3.79 2018-09-03 05:56:21 39.9
3 3362 -38.5 -3.79 2018-09-03 05:54:36 37.1
4 3363 -38.5 -3.79 2018-09-03 05:55:26 37.9
5 3528 -38.5 -3.80 2018-09-03 05:51:21 38.7
6 3540 -38.5 -3.80 2018-09-03 05:52:51 34.6
也欢迎使用线性插值以外的方法。谢谢。
我设法使用 zoo::na.approx 函数自行解决了这个问题。首先，您需要计算两个数据集中从同一点开始的连续点之间的累积距离。在我的示例中，我的 GPS 点按顺序排在第一位，因此我将选择第一个 GPS 点作为两个数据集的 "starting point"。
# I just realized that the datasets are growing in opposite direction, so I'll flip the stops
stops <- map_df(stops, rev)
# Function to calculate distance from previous point
get.dist <- function(lon, lat) geosphere::distHaversine(tail(cbind(lon,lat),-1),head(cbind(lon,lat),-1))
# Calculate cumulative distance of gps points (points must be ordered by time)
gps_points <- gps_points %>%
mutate(dist = c(0, cumsum(get.dist(lon, lat))))
# Input first GPS point and calculate cumulative distance of stops (also must be ordered)
stops <- gps_points %>%
# Select only the first point
slice(1) %>%
# Select columns to match the stops dataset
mutate(stop_id = NA) %>%
select(stop_id, lon, lat) %>%
# Input the stop points
rbind(stops) %>%
# Calculate cumulative dist
mutate(dist = c(0, cumsum(get.dist(lon, lat))))
# Interpolate ------------------------------
x <- gps_points$hora
y <- gps_points$dist
# to which position we want to interpolate? to the stops!
xout <- stops$dist
interp <- as.POSIXct(zoo::na.approx(x ,
y,
xout = xout,
ties = "ordered",
rule = 2),
origin = "1970-01-01")
# Put it together
stops_interp <- stops %>%
# Input the interpolated times
mutate(hora = interp) %>%
# Delete the first row that was inputed from the GPS
slice(-1)
stops_interp
stop_id lon lat dist hora
1 3528 -38.50404 -3.802504 801.8295 2018-09-03 05:51:21
2 4629 -38.50580 -3.798711 1266.8446 2018-09-03 05:52:18
3 3540 -38.50706 -3.796027 1596.9692 2018-09-03 05:52:51
4 3543 -38.50816 -3.793663 1887.0162 2018-09-03 05:53:43
5 3542 -38.50894 -3.791997 2091.8550 2018-09-03 05:54:05
6 3362 -38.50984 -3.790069 2328.6658 2018-09-03 05:54:38
7 3363 -38.51100 -3.787631 2628.9530 2018-09-03 05:55:27
8 3346 -38.51215 -3.785157 2932.6796 2018-09-03 05:56:18
9 5083 -38.51390 -3.781432 3390.5980 2018-09-03 05:57:15
10 3215 -38.51531 -3.776950 3913.4700 2018-09-03 05:58:22
11 4873 -38.51677 -3.771828 4506.1116 2018-09-03 05:59:52
第hora列表示插值时间。
更进一步
但是，这种插值是基于点之间的欧氏距离。如果我有 GPS 点和停靠点的预期路径（如道路）怎么办？在这种情况下，我建议使用 rgeos::gProject 函数来计算路径中 GPS 和停靠点的累积距离。