更新日期 20190506

1 问题

已知巴黎的地理坐标为东经2°20'14''，北纬48°50'11''，华盛顿的地理坐标为西经 77°03'56''，北纬 38°55'17''，求两城间的距离（单位千米）。

2 分析：

已知经纬度，求两点间的地理距离，不能直接用经纬度求平面距离，而是要计算球面距离（Great Circle Distance）。

如果对精度要求不高，可以用低精度公式，计算过程相对简单。如果对经度要求较高， 则要用高精度公式。高精度公式一般用于天文计算、大地测量等。在距离较远的城市之间， 两种方法所得结果可相差几十千米。

当然，现在全部通过计算机程序计算，都是直接使用高精度方法了。本文参考比利时天文学家Jean Meeus (1991)的《天文算法》81-86页的公式，给出相应的R代码。

2.1 将经纬度的度、分、秒转换为十进制

定义函数deg2dec完成这一转换:

deg2dec <- function(h,m,s){   
   if(h < 0){
       m = - m
       s = -s
   }
   res = h + m/60 + s/3600
   return(res)
}

## 巴黎
L1 = deg2dec(-2,20,14)
phi1 = deg2dec(48, 50, 11)

## 华盛顿
L2 = deg2dec(77,03,56)
phi2 = deg2dec(38,55,17)

2.2 低精度公式

由于地球的扁率较小，可近似看作球体，利用球面三角公式即可求出距离，R代码如下:

geodist <- function(L1, phi1, L2, phi2){
    Ri = 6371
    d = (Ri * acos(sin(phi1 * (pi / 180)) *
                      sin(phi2 * (pi / 180)) +
                      cos(phi1 * (pi / 180)) *
                      cos(phi2 * (pi / 180)) *
                      cos((L1 - L2) * (pi / 180))))
    res <- round(d, 3)    
    return(res)
}

geodist(L1, phi1, L2, phi2)

## [1] 6165.597

因此，如果假设地球是球体，华盛顿和巴黎之间的距离是6165.597km

2.3 高精度公式

假设地球是椭球，a是长轴，f是扁率，公式如下：

hageodist <- function(L1, phi1, L2, phi2){
   a = 6378.14
   f = 1/298.257
   F = (phi1 + phi2)/2
   G = (phi1 - phi2)/2
   ramda <- (L1 - L2)/2
   S = (sin(G * pi / 180)^2)*(cos(ramda * pi / 180)^2) + 
       (cos(F * pi / 180)^2)*(sin(ramda * pi / 180)^2) 
   C = (cos(G * pi / 180)^2)*(cos(ramda * pi / 180)^2) + 
       (sin(F * pi / 180)^2)*(sin(ramda * pi / 180)^2) 
   omega = atan(sqrt(S / C))
   R = sqrt(S * C)/omega
   D = 2 * omega * a
   H1 = (3 * R - 1)/(2 * C)
   H2 = (3 * R + 1)/(2 * S)
   res = D * (1 + f * H1 * (sin(F * pi / 180)^2)*
               (cos(G * pi / 180)^2) -
               f * H2 * (cos(F * pi/180)^2)*
               (sin(G * pi / 180)^2))   
   return(round(res, 3))
}

hageodist(L1, phi1, L2, phi2)

## [1] 6181.628

因此，华盛顿和巴黎之间的距离是6181.628km。高精度和低精度公式的结果相差16km。

2.4 注意事项

Meeus这本书中，东经为负，西经为正，公式采用国际天文学联合会IAU1976规定的相关参数，其坐标系未指定。

注意，Meeus的算法中，经度与常用的WGS84坐标系正好相反。在WGS84坐标系中，东经为正，西经为负。现在通常所说的经纬度是用GPS获得的WGS84经纬度，所以最保险的方法是以下R脚本计算：

3 推荐的R程序包和脚本

3.1 sp包

library(sp)
longitude <- as.numeric(c(L1, L2))
latitude  <- as.numeric(c(phi1, phi2))
cities <- c("Paris", "Washington")
location <- cbind(longitude = -longitude, latitude)
row.names(location) <- cities
location <- data.frame(location)
coordinates(location) <- ~longitude+latitude
proj4string(location) <- CRS("+proj=longlat +datum=WGS84") # 明确使用WGS84坐标系
spDists(location)

##          [,1]     [,2]
## [1,]    0.000 6181.626
## [2,] 6181.626    0.000

3.2 geosphere包

xy <- rbind(paris =c(-L1, phi1), washington =c(-L2, phi2))
library(geosphere)
distm(xy)

##         [,1]    [,2]
## [1,]       0 6181632
## [2,] 6181632       0

4 在地图上绘制Great Circle Distance

library(maps)
library(geosphere)
map('world',col="grey", fill=FALSE, bg="white", lwd=0.05,interior = FALSE)
points(x = -L1, y = phi1, pch = 19)
points(x = -L2, y = phi2, pch = 19)
text(x = -L1, y = phi1 + 5, labels = "Paris", col = "blue")
text(x = -L2, y = phi2 + 5, labels = "Washington", col = "blue")

connection <- gcIntermediate(c(-L1, phi1), c(-L2, phi2), 
    n=100, addStartEnd = TRUE, breakAtDateLine=FALSE)             
lines(connection, col="red", lwd=2)

5 参考

• https://www.r-graph-gallery.com/how-to-draw-connecting-routes-on-map-with-r-and-great-circles/

• https://www.jessesadler.com/post/great-circles-sp-sf/