[PYTHON] Ungefährer Abstand zwischen zwei Punkten auf der Oberfläche eines rotierenden Ellipsoids (auf der Erdoberfläche)

Hallo. Die Berechnung des Abstands zwischen zwei Punkten auf der rotierenden elliptischen Oberfläche ist das umgekehrte Problem der Grundlinie [^ 1], aber ich habe die ungefähre Berechnung [^ 2] unter der Bedingung des kurzen Abstands versucht. Im folgenden Beispiel für die Bewertung der Approximationsgenauigkeit beträgt der Fehler der Entfernung von 12,457 km 1,3 mm.

[^ 1]: Dies kann mit GeographicLib (`Inverse ()`) berechnet werden. [^ 2]: Dies ist auch die Methode, die für die Anfangswertberechnung für die Newton-Methode in der inversen Problemberechnung von GeographicLib verwendet wird.

$ ./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)