创建具有测地距离的新列答案

【问题标题】：Create new column with geodesic distances创建具有测地距离的新列
【发布时间】：2015-06-05 20:31:32
【问题描述】：

如何向我的数据框中添加一个新列，其中包含每个观测值的测地距离？我的数据框中的行是对租用汽车的观察，包括旅行的起点（经度 + 纬度）和终点（经度 + 纬度）。

测地线距离使用以下 R 函数计算：

haversine<- function(long1, lat1, long2, lat2) {

  stopifnot(is.numeric(long1),
            is.numeric(lat1),
            is.numeric(long2),
            is.numeric(lat2),
            long1 > -180,
            long1 < 180,
            lat1 > -180,
            lat1 < 180,
            long2 > -180,
            long2 < 180,
            lat2 > -180,
            lat2 < 180  
    )

  long1 <- long1*pi/180
  lat1 <- lat1*pi/180
  long2 <- long2*pi/180
  lat2 <- lat2*pi/180

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

用这个公式可以计算距离：

distance <- haversine(longitude1,latitude1,longitude2,latitude2)

由于我的数据框包含超过 100,000 个观测值，我想在我的数据框中添加一列，其中包含汽车旅行的所有测地线距离。

我试过了： df["GeoDist"] <- haversine(LongitudeStart, LatitudeStart, LongitudeEnd, LatitudeEnd)

这导致了一个额外的列，但是该列中的所有距离都为零。当我使用前面的公式手动计算时，这些距离不为零。

【问题讨论】：

标签： r haversine

【解决方案1】：

您可以只使用 geosphere 包中的 distHaversine()。或者，要获得更准确的结果，请使用 distVincentyEllipsoid()。（这两个函数都以度为单位获取坐标，以米为单位返回距离。）

## Example data representing voyages of 1 and 2 degrees directly north from equator.
df <- data.frame(x_begin = c(0,0), y_begin = c(0,0), 
                 x_end   = c(0,0), y_end   = c(1,2))

library(geosphere)
distHaversine(df[, 1:2],  df[, 3:4])        ## Assumes a spherical Earth
# [1] 111319.5 222639.0
distVincentyEllipsoid(df[,1:2],  df[,3:4])  ## Much more accurate*
## [1] 110574.4 221149.5

## * According to Wikipedia, 1 degree of latitude at equator == 110.574 km 
## (http://en.wikipedia.org/wiki/Longitude#Length_of_a_degree_of_longitude)

【讨论】：

【解决方案2】：

您需要mapply，因为您正在传递向量参数。

# reproducible example
# 38.8895° N,  77.0352° W  # 
# 47.6204° N, 122.3491° W  # washington monument to space needle
# 40.6892° N,  74.0444° W  # 
# 37.8197° N, 122.4786° W  # statue of liberty to golden gate bridge
df <- data.frame(y1=c(  38.8895,  40.6892),
                 x1=c( -77.0352, -74.0444),
                 y2=c(  47.6204, 37.8197 ),
                 x2=c( -122.3491,-122.4786 ))

with(df, mapply(haversine, lat1=y1, long1=x1, lat2=y2, long2=x2) )
# [1] 3738.587 4129.629 # looks right to me

【讨论】：