已知经纬度计算两点间地理距离

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完成这一转换:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
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是扁率,公式如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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包

1
2
3
4
5
6
7
8
9
10
11
12
13
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包

1
2
3
4
5
6
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

1
2
3
4
5
6
7
8
9
10
11
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 参考