"Dabei ist w im Bogenmaß einzusetzen"

[Strawk]

Hallo!

Versuche gerade, die Vincenty-Formel in Python zu übersetzen. Bis zum groben Abstand "D" sieht es ganz vielversprechend aus, aber dazu fehlt noch eine Angabe von "w" im Bogenmaß. Wie kann ich einen Wert "w" "im Bogenmaß einsetzen"?

Hier der vorläufige Code:

Code: Alles auswählen

# -*- coding: utf-8 -*-
"""
Created on Fri Apr 26 10:45:34 2019

@author: Karl Kraft
"""
from __future__ import print_function, division
import numpy as np
import math
from math import radians

def vincenty(lat1, lat2, lon1, lon2):
    """ calculates the distance in kilometers between two geographical points
    """
    '''
    lat1 = radians(lat1)
    lat2 = radians(lat2)
    lon1 = radians(lon1)
    lon2 = radians(lon2)
    '''
        
    f = 1/298.257223563
    # print(f)
    a = 6378.137
    
    F = ((lat1+lat2)/2)

    G = ((lat1-lat2)/2)

    l = ((lon1-lon2)/2)

    S = ((math.sin(G))**2) * ((math.cos(l))**2) + ((math.cos(F))**2) * ((math.sin(l))**2)

    C = ((math.cos(G))**2) * ((math.cos(l))**2) + ((math.sin(F))**2) * ((math.sin(l))**2)

    w = np.arctan(math.sqrt(S/C))

    D = 2 * 0.0151440 * a

    print(w)    
    
    return D

if __name__ == "__main__":
    '''
    westlicher Ort: Position 1
    östlicher Ort: Position 2
    '''
    '''
    # Madrid (1) - Budapest (2)
    lat1 = 40.4167754
    lat2 = 47.497912
    lon1 = -3.7037902
    lon2 = 19.040235
    '''
    '''
    # New York (1) - Aachen (2)
    lat1 = 40.712784
    lat2 = 50.7753455
    lon1 = -74.005941
    lon2 = 6.0838868
    '''
    '''
    # Berlin (1) - Pretoria (2)
    lat1 = 52.520007
    lat2 = -25.746111
    lon1 = 13.404954
    lon2 = 28.188056
    '''
    
    # Pretoria (1) - Moskau (2)
    lat1 = -25.746111
    lat2 = 55.755826
    lon1 = 28.188056
    lon2 = 37.617300
    
    '''
    # Aachen (1) - Mönchengladbach (2)
    lat1 = 50.7753455
    lat2 = 51.1804572
    lon1 = 6.0838868
    lon2 = 6.4428041
    '''
    result = vincenty(lat1, lat2, lon1, lon2)
    print("km: %.1f" % (result), "\t")

Grüße
Strawk

[Sirius3]

Du mußt nichts einsetzen, da Du ja `w` als Ergebnis bekommst. arctan2 wäre da auch besser.
Und wie immer, Zahlen nicht irgendwo im Code verstecken, sondern als Konstanten am Anfang definieren.

[Strawk]

Artan2 braucht aber 2 Werte. Und meinst du D=2*a ?

[Sirius3]

Genau, 2 Werte sind besser als einer, weil man durch 0 nicht teilen kann.

Warum benutzt Du numpy und math vermischt?

[Strawk]

Und woher soll der zweite Wert kommen? Weil ich es eben nicht besser kann.

[Sirius3]

Du hast doch zwei Variablen.

[Strawk]

Ich habe arctan von Wurzel aus S/C. Also eine.

[Strawk]

Modifizierter Code:

Code: Alles auswählen

# -*- coding: utf-8 -*-
"""
Created on Fri Apr 26 10:45:34 2019

@author: Karl Kraft
"""
from __future__ import print_function, division
import numpy as np

def vincenty(lat1, lat2, lon1, lon2):
    """ calculates the distance in kilometers between two geographical points using
        the Vincenty-formula
    """
    '''
    lat1 = radians(lat1)
    lat2 = radians(lat2)
    lon1 = radians(lon1)
    lon2 = radians(lon2)
    '''
        
    f = 1/298.257223563
    
    a = 6378.137
    
    F = ((lat1+lat2)/2)

    G = ((lat1-lat2)/2)

    l = ((lon1-lon2)/2)

    S = ((np.sin(G))**2) * ((np.cos(l))**2) + ((np.cos(F))**2) * ((np.sin(l))**2)

    C = ((np.cos(G))**2) * ((np.cos(l))**2) + ((np.sin(F))**2) * ((np.sin(l))**2)

    w = np.arctan(np.sqrt(S/C))

    D = 2 * w * a

    T = (np.sqrt(S*C))/w
    
    H1 = (3*T-1) / (2*C)
    
    H2 = (3*T+1) / (2*S)
    
    s = D * ((1+f*H1*(np.sin(F))**2) * (np.cos(G)**2) - f * H2 * (np.cos(F)**2) * (np.sin(G)**2))
    
    return s

if __name__ == "__main__":
    '''
    westlicher Ort: Position 1
    östlicher Ort: Position 2
    '''
    
    # Madrid (1) - Budapest (2)
    lat1 = 40.4167754
    lat2 = 47.497912
    lon1 = -3.7037902
    lon2 = 19.040235
    
    '''
    # New York (1) - Aachen (2)
    lat1 = 40.712784
    lat2 = 50.7753455
    lon1 = -74.005941
    lon2 = 6.0838868
    '''
    '''
    # Berlin (1) - Pretoria (2)
    lat1 = 52.520007
    lat2 = -25.746111
    lon1 = 13.404954
    lon2 = 28.188056
    '''
    '''
    # Pretoria (1) - Moskau (2)
    lat1 = -25.746111
    lat2 = 55.755826
    lon1 = 28.188056
    lon2 = 37.617300
    '''
    '''
    # Aachen (1) - Mönchengladbach (2)
    lat1 = 50.7753455
    lat2 = 51.1804572
    lon1 = 6.0838868
    lon2 = 6.4428041
    '''
    result = vincenty(lat1, lat2, lon1, lon2)
    print("km: %.1f" % (result), "\t")

[Sirius3]

so:

Code: Alles auswählen

EARTH_RADIUS = 6378.137
EARTH_FLATTENING = 1/298.257223563

def vincenty(lat1, lat2, lon1, lon2):
    """ calculates the distance in kilometers between two geographical points using
        the Vincenty-formula
    """
    lat1 = np.radians(lat1)
    lat2 = np.radians(lat2)
    lon1 = np.radians(lon1)
    lon2 = np.radians(lon2)
    f = (lat1 + lat2) / 2
    g = (lat1 - lat2) / 2
    l = (lon1 - lon2) / 2
    sin_g2, cos_g2 = np.sin(g) ** 2, np.cos(g) ** 2
    sin_f2, cos_f2 = np.sin(f) ** 2, np.cos(f) ** 2
    sin_l2, cos_l2 = np.sin(l) ** 2, np.cos(l) ** 2
    s = sin_g2 * cos_l2 + cos_f2 * sin_l2
    c = cos_g2 * cos_l2 + sin_f2 * sin_l2
    s_sqrt = s ** 0.5
    c_sqrt = c ** 0.5
    w = np.arctan2(s_sqrt, c_sqrt)
    d = 2 * w * EARTH_RADIUS
    t = 3 * (s_sqrt * c_sqrt) / w    
    h1 = (t - 1) / (2*c)
    h2 = (t + 1) / (2*s)
    return d * ((1 + EARTH_FLATTENING * h1 * sin_f2) * cos_g2
        - EARTH_FLATTENING * h2 * cos_f2 * sin_g2)

[Strawk]

Danke!
Das ist außerordentlich hilfreich!