計算距離

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功能（包直接讀取這些數據類型爲）