经纬度计算距离
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){
s = -s
res = h + m/60 + s/3600
return(res)
## 巴黎
2,20,14)
48, 50, 11)
## 华盛顿
77,03,56)
38,55,17)2.2 低精度公式由于地球的扁率较小，可近似看作球体，利用球面三角公式即可求出距离，R代码如下:
geodist <- function(L1, phi1, L2, phi2)
Ri = 6371
180)) *
180)) +
180)) *
180)) *
180))))
3)
return(res)
## [1] 6165.597因此，如果假设地球是球体，华盛顿和巴黎之间的距离是6165.597km
2.3 高精度公式假设地球是椭球，a是长轴，f是扁率，公式如下：
hageodist <- function(L1, phi1, L2, phi2){
6378.14
1/298.257
F = (phi1+phi2)/2
2
2
180)^2)*(cos(ramda*pi/180)^2)+
F*pi/180)^2)*(sin(ramda*pi/180)^2)
180)^2)*(cos(ramda*pi/180)^2)+
F*pi/180)^2)*(sin(ramda*pi/180)^2)
R = sqrt(S*C)/omega
2*omega*a
3*R-1)/(2*C)
3*R+1)/(2*S)
res = D*(1 +
f * H1*(sin(F*pi/180)^2)*
(cos(G*pi/180)^2)-
F*pi/180)^2)*
180)^2))
return(round(res,3))
## [1] 6181.628因此，华盛顿和巴黎之间的距离是6181.628km。高精度和低精度公式的结果相差16km。
2.4 注意事项Meeus这本书中，东经为负，西经为正，公式采用国际天文学联合会（IAU）1976规定的相关参数，其坐标系未指定。
注意，Meeus的算法中，经度与常用的WGS84坐标系正好相反。在WGS84坐标系中，东经为正，西经为负。现在通常所说的经纬度是用GPS获得的WGS84经纬度，所以最保险的方法是以下R脚本计算：
3 推荐的R程序包和脚本3.1 sp包library(sp)
latitude <- as.numeric(c(phi1, phi2))
“Paris”, “Washington”)
row.names(location) <- cities
coordinates(location) <- ~longitude+latitude
# 明确使用WGS84坐标系
proj4string(location) <- CRS(“+proj=longlat +datum=WGS84”)
spDists(location)## [,1] [,2]
## [2,] 6181.626 0.0003.2 geosphere包xy <- rbind(paris =c(-L1, phi1),
washington =c(-L2, phi2))
library(geosphere)
## [,1] [,2]
## [2,] 6181632 04 在地图上绘制Great Circle Distancelibrary(maps)
library(geosphere)
map(‘world’,col=”grey”, fill=FALSE,
bg=”white”, lwd=0.05,
interior = FALSE)
points(x = -L1, y = phi1, pch = 19)
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/
