Hello. The calculation of the distance between two points on the spheroidal surface is the inverse problem of the geodesic [^ 1], but I tried the approximate calculation [^ 2] under the short distance condition. In the approximation accuracy evaluation example below, the error of the distance of 12.457km is 1.3mm.
[^ 1]: This can be calculated using GeographicLib (```Inverse () `` `). [^ 2]: This is the method that is also used for the initial value calculation for Newton's method in the inverse problem calculation of GeographicLib.
$ ./geodesic_inverse_approx.py 60.0 60.1 0.1
input: lat1= 60.000000 lat2= 60.100000 lam12= 0.100000
output: az1= 26.526 az2= 26.612 s12= 12456.777
error: az1= 1.2e-05 az2= 1.2e-05 s12= 1.3e-03
geodesic_inverse_approx.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import division
from __future__ import print_function
import sys
import numpy as np
from easydict import EasyDict as edict
from geographiclib.geodesic import Geodesic
RE = 6378137.0 # IS-GPS
FE = (1/298.257223563) # IS-GPS
E2 = FE * (2 - FE)
DEGREE = np.pi/180
geod = Geodesic.WGS84
# "Algorithms for geodesics" Charles F. F. Karney (2013)
# 5. Starting point for Newton’s method
# solving for the great circle on the auxiliary sphere
def geodesic_inverse_approx(beta1, beta2, lam12):
cosBeta = (beta1.cos + beta2.cos)/2
wm = np.sqrt(1-E2*cosBeta**2)
omg12 = angle_approx(lam12/wm)
az1, r = angle(beta1.cos*beta2.sin - beta1.sin*beta2.cos*omg12.cos, beta2.cos*omg12.sin)
az2, _ = angle(-beta1.sin*beta2.cos + beta1.cos*beta2.sin*omg12.cos, beta1.cos*omg12.sin)
sig12 = arg_approx(beta1.sin*beta2.sin + beta1.cos*beta2.cos*omg12.cos, r)
s12 = sig12 * wm
return az1, az2, s12
def reducedLatitude(lat):
chi = np.tan(lat/2)
beta, _ = angle(1-chi**2, 2*chi*(1-FE))
return beta
def arg(*args):
if len(args) == 1:
c, s = args[0].cos, args[0].sin
else:
c, s = args[0], args[1]
return np.arctan2(s, c)
def arg_approx(c, s):
x = s/c
return x*(1-x**2/3)
def angle(*args):
if len(args) == 1:
x = np.tan(args[0]/2)
return angle(1-x**2, 2*x)
r = np.hypot(args[0], args[1])
c, s = args[0]/r, args[1]/r
return edict({'cos': c, 'sin': s}), r
def angle_approx(x):
return edict({'cos': 1-x**2/2, 'sin': x})
def print_geodesic_inverse(lat1, lat2, lam12):
print("input: lat1=%10.6f lat2=%10.6f lam12=%13.6f" % (lat1, lat2, lam12))
beta1 = reducedLatitude(lat1*DEGREE)
beta2 = reducedLatitude(lat2*DEGREE)
az1, az2, s12 = geodesic_inverse_approx(beta1, beta2, lam12*DEGREE)
az1, az2, s12 = arg(az1)/DEGREE, arg(az2)/DEGREE, s12*RE
print("output: az1=%7.3f az2=%7.3f s12=%10.3f" % (az1, az2, s12))
g = edict(geod.Inverse(lat1, 0, lat2, lam12))
print("error: az1=%8.1e az2=%8.1e s12=%10.1e" % (az1-g.azi1, az2-g.azi2, s12-g.s12))
return
def main():
args = sys.argv[1:]
if len(args) == 3:
print_geodesic_inverse(*map(lambda x:np.float(x), args))
return
if __name__ == '__main__':
main()
exit(0)
Recommended Posts