计算距离

Question

我确实有GPS数据的数据帧，我需要之间的两行（当前和以前）来计算距离

  id     time  lat  long heartrate altitude 
1 20130424.tcx 2013-04-24T04:53:22Z 50.024818 14.522724  151  <NA> 
2 20130424.tcx 2013-04-24T04:53:26Z 50.024818 14.522724  96  <NA> 
3 20130424.tcx 2013-04-24T04:53:30Z 50.024818 14.522724  104  <NA> 
4 20130424.tcx 2013-04-24T04:53:34Z 50.024818 14.522724  107  <NA> 
5 20130424.tcx 2013-04-24T04:53:38Z 50.024818 14.522724  108  <NA> 
6 20130424.tcx 2013-04-24T04:53:42Z 50.024818 14.522724  112 372.0 
7 20130424.tcx 2013-04-24T04:53:46Z 50.024818 14.522724  151 372.0 
8 20130424.tcx 2013-04-24T04:53:47Z 50.024677 14.522874  151 356.0 
9 20130424.tcx 2013-04-24T04:53:50Z 50.024677 14.522874  118 356.0 
10 20130424.tcx 2013-04-24T04:53:54Z 50.024677 14.522874  118 356.0 
11 20130424.tcx 2013-04-24T04:53:58Z 50.024464 14.522917  147 358.0 
12 20130424.tcx 2013-04-24T04:54:02Z 50.024464 14.522917  144 358.0 
13 20130424.tcx 2013-04-24T04:54:06Z 50.024269 14.522853  150 367.0 
14 20130424.tcx 2013-04-24T04:54:10Z 50.024269 14.522853  152 367.0 
15 20130424.tcx 2013-04-24T04:54:13Z 50.024002 14.522874  152 380.0 

我能够将数据加入到本身并取得前一行每一行（有可能是简单的解决方案）：

library(sqldf) 
mydft = mydf[-nrow(mydf),] 
mydft$id = seq_along(mydft$id) +1 
mydf$id = seq_along(mydf$id) 
mydft2 <- sqldf("select a.*, b.lat as lat2, b.long as long2 from mydf a left join mydft b using (id)") 

如何我现在可以计算出从列lat, long, lat2, long2距离是多少？我尝试过描述的方法here：

R <- 6371 # Earth mean radius [km] 
mydft2$delta.long <- (mydft2$long2 - mydft2$long) 
mydft2$delta.lat <- (mydft2$lat2 - mydft2$lat) 
mydft2$a <- sin(mydft2$delta.lat/2)^2 + cos(mydft2$lat) * cos(mydft2$lat2) * sin(mydft2$delta.long/2)^2 
mydft2$c <- 2 * asin(min(1,sqrt(mydft2$a))) 
mydft2$d = R * c 

但是这个回退只能列出NAs。

Answer 1

简单实现的一种方法是使用R设施进行空间分析。对于这个例子，我们可以在优秀的sp包中使用spDistN1函数。

的第一步是将数据转换成SpatialPoints（或SpatialPointsDataFrame），我假定在该示例中，您的点在地理上的突起（longlat与WSG84基准面）

txt <- "   id     time  lat  long heartrate altitude 
20130424.tcx 2013-04-24T04:53:22Z 50.024818 14.522724  151  <NA> 
20130424.tcx 2013-04-24T04:53:26Z 50.024818 14.522724  96  <NA> 
20130424.tcx 2013-04-24T04:53:30Z 50.024818 14.522724  104  <NA> 
20130424.tcx 2013-04-24T04:53:34Z 50.024818 14.522724  107  <NA> 
20130424.tcx 2013-04-24T04:53:38Z 50.024818 14.522724  108  <NA> 
20130424.tcx 2013-04-24T04:53:42Z 50.024818 14.522724  112 372.0 
20130424.tcx 2013-04-24T04:53:46Z 50.024818 14.522724  151 372.0 
20130424.tcx 2013-04-24T04:53:47Z 50.024677 14.522874  151 356.0 
20130424.tcx 2013-04-24T04:53:50Z 50.024677 14.522874  118 356.0 
20130424.tcx 2013-04-24T04:53:54Z 50.024677 14.522874  118 356.0 
20130424.tcx 2013-04-24T04:53:58Z 50.024464 14.522917  147 358.0 
20130424.tcx 2013-04-24T04:54:02Z 50.024464 14.522917  144 358.0 
20130424.tcx 2013-04-24T04:54:06Z 50.024269 14.522853  150 367.0 
20130424.tcx 2013-04-24T04:54:10Z 50.024269 14.522853  152 367.0 
20130424.tcx 2013-04-24T04:54:13Z 50.024002 14.522874  152 380.0" 

gpsdat <- read.table(text = txt, header = TRUE, na.strings = "<NA>") 
str(gpsdat) 

require(sp) 
coordinates(gpsdat) <- ~ long + lat 
proj4string(gpsdat) <- CRS("+proj=longlat +datum=WGS84") 

在第二步骤中，我们现在可以使用spDistN1功能

sapply(seq_along(gpsdat[-1, ]), function(i) 
     spDistsN1(pts = gpsdat[i, ], pt = gpsdat[i+1, ], longlat = TRUE)) 

[1] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
[7] 0.019004 0.000000 0.000000 0.023875 0.000000 0.022155 
[13] 0.000000 0.029716 

取决于我们使用GPS数据的类型，可以使用readOGR功能（包直接读取这些数据类型为）