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) 

Recommended Posts

Distance approximative entre deux points à la surface d'un ellipsoïde en rotation (à la surface de la terre)
Approximation par la méthode des moindres carrés d'un cercle à deux points fixes
[Rust] Lire les données csv de latitude et de longitude pour trouver la distance entre deux points
Calculer la probabilité de valeurs aberrantes sur les moustaches de la boîte
Modélisation-estimation de Bayes de la différence entre les deux groupes-
Points Python du point de vue d'un programmeur en langage C
Une réflexion sur la visualisation du champ d'application du modèle de prédiction
Remarque sur le comportement par défaut de collate_fn dans PyTorch
Un résumé approximatif des différences entre Windows et Linux
Script Python qui compare le contenu de deux répertoires
Partagez le volume physique d'un PC à double démarrage entre les systèmes d'exploitation.
Sous Linux, l'horodatage d'un fichier est un peu dépassé.
Trouvez le rang de la matrice dans le monde XOR (rang de la matrice sur F2)
Une commande pour vérifier facilement la vitesse du réseau sur la console
Recherche de quartier pour un ensemble de points (coordonnées géographiques) sur une sphère
Trouvez la distance (en tenant compte de la rondeur de la terre) de la latitude et de la longitude.
Obtenez le nombre de lecteurs d'articles sur Mendeley en Python