R/great.circle.distance.vector.vector.r

      
  great.circle.distance.vector.vector = function (loc1, loc2, R) {
    if (is.null(R)) R = 6367.436  # radius of earth (geometric mean) km
    # if R=1 then distances in radians
    if (missing(loc2)) loc2 <- loc1
    pi180 = pi/180
    coslat1 = cos(loc1$lat * pi180)
    sinlat1 = sin(loc1$lat * pi180)
    coslon1 = cos(loc1$lon * pi180)
    sinlon1 = sin(loc1$lon * pi180)
    coslat2 = cos(loc2$lat * pi180)
    sinlat2 = sin(loc2$lat * pi180)
    coslon2 = cos(loc2$lon * pi180)
    sinlon2 = sin(loc2$lon * pi180)
    pp = diag(  cbind(coslat1 * coslon1, coslat1 * sinlon1, sinlat1) %*%
         t(cbind(coslat2 * coslon2, coslat2 * sinlon2, sinlat2)))
    
    d = R * acos(ifelse(pp > 1, 1, pp))
    
    return(d)
}
jae0/stmdat documentation built on May 28, 2019, 11 p.m.