在计算地球上两点之间的距离时,由于地球是近似为一个三维椭球体,一般是不能直接地理经纬度相减来计算.而将经纬度转到直角坐标下,在对xyz分别做差,并计算平方和最后开根,也不能反应真实地表距离.在讨论球面距离时,一般是获取两点之间的球面最短距离,也就是同时过两点以及球心的大圆上连接两点的较短的圆弧的距离.

在 Julia 中 Geodesy.jl 提供了表示地理空间的一些数据结构以及不同坐标之间的转换函数. 比如 LLA 就可以表示地理坐标下的经纬度高程, ECEF 则是地心地固坐标下的 xyz. 整体而言这个包是比较简单的,不过也算够用.但是其中的distance ref,是直接转到 ECEF 然后计算 norm, 这个其实在大部分时候都不是程序使用者想要的. 比如下面计算距离的代码, distance 地球椭球体(默认wgs84)的三个轴长.

using Geodesy
a = LLA(0,  0,0)
b = LLA(0,180,0)

distance(a,b)   # 1.2756274e7
distance(a,b)/2 # 6.378137e6

c = LLA(0, 90,0)
d = LLA(0,-90,0)

distance(c,d)   # 1.2756274e7
distance(c,d)/2 # 6.378137e6

e = LLA(-90,0,0)
f = LLA( 90,0,0)

distance(e,f)   # 1.2713504628490359e7
distance(e,f)/2 # 6.356752314245179e6

在地球上,一般是不太可能直接沿着直线(穿过部分岩石圈,甚至过地心),到达另一点的. 常规操作应该是由地表一点沿着地表或者近地表到达另一点.而在球面上这所有可能路径中,距离最短就是大圆上较短的圆弧.

计算球面的大圆距离的方式,主要有 Haversine formula, Vincenty’s formulae, 还有GeographicLib. Haversine 最直接,不过适用的是球体, Vincenty 的考虑椭球扁率,但是在某些区域无法收敛,最后的 GeographicLib 提供的算法应该是既考虑地球形状又具有很好的收敛性的. 但是 GeographicLib 本身是 C++ 库,虽然也提供诸如 Python 的包,但是不包括 Julia. 有一些个人的转写,不过没有到Julia管理系统中. 而Vincenty的收敛问题,主要是对趾点附近,在考虑大范围的时候通常都会遇到. 最后考虑地球的扁率较小, Haversine 是一个可以接受的算法.

对于像 Haversine 一样较为通用的算法, Julia 还是有现成的库可以使用的. Distances.jl 中就提供了 haversine.

最初我是在 discourse.julialang.org 中查到相关的描述的. 根据标记为 Solution 的回答,也就是如下的代码:

using Distances
julia> l1 = (-27.468937, 153.023628)
julia> l2 = (-27.465933, 153.025900)
julia> haversine(l1, l2, 6372.8)
0.39054922275889736

我一直以为 haversine 输入是纬度在前,经度在后的 TupleArray.而且有一段程序就是基于这个认识,并且运行了一段时间了. 最近新写的代码中又需要计算这个距离,但是结果总是和预想的不一样. 于是找到如下的一个确定算例,来验证 haversine 的结果. 输入当然还是有区别,最后直接到 Github 上找代码,发现上面直接写着:

       # longitudes
        Δλ = deg2rad(y[1] - x[1])

问题就是输出应当是经度在前,维度在后. 之前的程序能够勉强运行的原因是输入的两个坐标经纬度都互换了,输出还是能够一定的远近程度.而 discourse.julialang.org 中的 l1, l2 是两个十分相近的位置,所以结果基本上也差不多. 我最初想是不是加上检测语句确定经度在-180到360之间,然后纬度在-90到90之间.不过对于这个函数,实际不是针对地球而是球体的.应该是如果用在地球上的话,由用户来判断.

借用网上的代码的时候,最好还是看下函数体,然后找几个算例试试,看结果是不是想要的.

发表评论

电子邮件地址不会被公开。 必填项已用*标注