[PYTHON] Distance approximative entre deux points à la surface d'un ellipsoïde en rotation (à la surface de la terre)

Bonjour. Le calcul de la distance entre deux points sur la surface elliptique tournante est le problème inverse de la ligne de sol [^ 1], mais j'ai essayé le calcul approximatif [^ 2] sous la condition de courte distance. Dans l'exemple d'évaluation de la précision d'approximation ci-dessous, l'erreur de la distance de 12,457 km est de 1,3 mm.

[^ 1]: Ceci peut être calculé en utilisant GeographicLib (```Inverse () `` `). [^ 2]: C'est également la méthode utilisée pour le calcul de la valeur initiale de la méthode Newton dans le calcul du problème inverse de 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)