张金龙的博客分享 http://blog.sciencenet.cn/u/zjlcas 物种适应性、分布与进化

博文

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

已有 32706 次阅读 2010-9-21 14:52 |个人分类:科研笔记|系统分类:科研笔记|关键词:学者| R语言, 距离, 地理, 纬度, 经度

更新日期 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://m.sciencenet.cn/blog-255662-365318.html

上一篇:如何计算两日期相差的天数?
下一篇:CRAN Task View: 在R进行系统发育比较分析翻译

3 高建国 杨华磊 yjw413

该博文允许注册用户评论 请点击登录 评论 (19 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-28 18:57

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部