Haversine-Formel liefert falsche Werte

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

Hallo!
Die Haversine-Formel, zu finden unter https://en.wikipedia.org/wiki/Haversine_formula in Python zu übertragen, habe ich versucht. Die zurückgelieferten Ergebnisse sind aber falsch. Warum?

Code: Alles auswählen

from __future__ import print_function, division

import numpy as np

def haversine(phi_1, phi_2, lambda_1, lambda_2):
    """ calculates the distance in kilometers between two geographical points
    """
    leftParenthesis = ((phi_2 - phi_1)/2)
    rightParenthesis = ((lambda_2 - lambda_1)/2)

    allUnderSquareRoot = ((np.sin(leftParenthesis))**2) + np.cos(phi_1)*np.cos(phi_2)*(np.sin(rightParenthesis))**2
    squareRootResult = np.sqrt(allUnderSquareRoot)
    arcSin = np.arcsin(squareRootResult)
    r = 6371
    result = 2 * r * arcSin

    return result

if __name__ == "__main__":
    '''
    westlicher Ort: Position 1
    östlicher Ort: Position 2
    '''
    
    '''
    # Madrid (1) - Budapest (2)
    PHI_1 = 40.4167754
    PHI_2 = 47.497912
    LAMBDA_1 = -3.7037902
    LAMBDA_2 = 19.040235
    '''

    # New York (1) - Aachen (2)
    phi_1 = 40.712784
    phi_2 = 50.7753455
    lambda_1 = -74.005941
    lambda_2 = 6.0838868   

    result = haversine(phi_1, phi_2, lambda_1, lambda_2)
    print(result)
Ich programmiere erfolglos, also bin ich nicht.
Benutzeravatar
ThomasL
User
Beiträge: 1366
Registriert: Montag 14. Mai 2018, 14:44
Wohnort: Kreis Unna NRW

Schau mal hier rein, eventuell hilfreich.
viewtopic.php?f=1&t=19980
Ich bin Pazifist und greife niemanden an, auch nicht mit Worten.
Für alle meine Code Beispiele gilt: "There is always a better way."
https://projecteuler.net/profile/Brotherluii.png
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

Hallo!
Das funktioniert, ja. Hilft mir aber nicht, den Fehler in meiner Formel zu finden, da jene Formel ganz anders aufgebaut ist. Schönen Sonntag!

Grüße
Strawk
:)
Ich programmiere erfolglos, also bin ich nicht.
Sirius3
User
Beiträge: 17741
Registriert: Sonntag 21. Oktober 2012, 17:20

@Strawk: der Fehler im verlinkten Beitrag ist der selbe, wie bei Dir. Trigonometrische Funktionen erwarten den Winkel in Radians.
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

Hallo!

Funktioniert jetzt!

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
from math import radians

import numpy as np

def haversine(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)
    
    leftParenthesis = ((lat2 - lat1)/2)
    rightParenthesis = ((lon2 - lon1)/2)

    allUnderSquareRoot = ((np.sin(leftParenthesis))**2) + np.cos(lat1)*np.cos(lat2)*(np.sin(rightParenthesis))**2
    squareRootResult = np.sqrt(allUnderSquareRoot)
    arcSin = np.arcsin(squareRootResult)
    r = 6371
    result = 2 * r * arcSin

    return result

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
    
    result = haversine(lat1, lat2, lon1, lon2)
    print("km: %.1f" % (result), "\t")
Ich programmiere erfolglos, also bin ich nicht.
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

@ThomasL: Das war ein heißer Tipp, der zum Erfolg führte. Danke!
Grüße
Strawk
:wink:
Ich programmiere erfolglos, also bin ich nicht.
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

Hallo!

Ich möchte die Abstände aller weltweiten Flughäfen zu allen berechnen. Das dauert auf meinem Rechner mit der geopy-Modul-Funktion „distance“ (Methode „great_circle“) ca. 10 Minuten bei folgenden Werten:

distances_airports_worldwide.py:
country == "All"
number of calculations: 25,801,336 (25 Mio.)
km: 218,753,071,512.4 (219 Mrd. km)
Computation of distances took 554.27 secs = ca. 9 min 15 secs

Ich habe jetzt bei unten stehendem Code alle Breitengrade (lats) und alle Längengrade (lons) als numpy-arrays vorliegen. Nun möchte ich die Haversine-Funktion über den ganzen Vektor anwenden, um die Rechenzeit erheblich zu verkürzen. Hierzu erbitte ich eure Tipps.

Danke.

Grüße
Strawk

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
from math import radians
# from gpx_project import dist
from distances_airports_worldwide import get_airport_positions

def haversine(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)
    
    leftParenthesis = ((lat2 - lat1)/2)
    rightParenthesis = ((lon2 - lon1)/2)

    allUnderSquareRoot = ((np.sin(leftParenthesis))**2) + np.cos(lat1)*np.cos(lat2)*(np.sin(rightParenthesis))**2
    squareRootResult = np.sqrt(allUnderSquareRoot)
    arcSin = np.arcsin(squareRootResult)
    # before: 6371 km
    middleEarthRadiusKm = 6378.1
    middleEarthRadiusM = middleEarthRadiusKm * 1000
    result = 2 * middleEarthRadiusM * arcSin

    return result

def distvec(lat0, lon0, lats, lons):
    """calculates the distance of all airports worldwide to each other in kilometers
    """
    myResult = haversine(lat0, lon0, lats, lons)
    '''
    lst = [2, 3, 7.9, 3.3, 6.9, 0.11, 10.3, 12.9]
    v = np.array(lst)
    print(v + 2)
    '''
    return myResult
           
if __name__ == "__main__":
    coords = get_airport_positions("data\\airports.dat", "All")
    lats = [t[0] for t in coords]
    lons = [t[1] for t in coords]
    lats = np.array(lats)
    lons = np.array(lons)
    myResult = distvec(-6.081689834590001, 145.391998291, lats, lons)
    print(myResult)
    
    '''
    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 = haversine(lat1, lat2, lon1, lon2)
    print("km: %.1f" % (result), "\t")
    '''
    '''
    result_haversine = haversine(lat1, lat2, lon1, lon2) / 1000
    result_great_circle = dist((lat1, lon1), (lat2, lon2), method="great_circle") / 1000
    rel_error = abs((result_great_circle - result_haversine) / result_great_circle)
    rel_error_percent = rel_error * 100
    
    print("result haversine km: %.1f" % (result_haversine), "\t")
    print("result great_circle km: %.1f" % (result_great_circle), "\t")
    print("result rel_error km: %.7f" % (rel_error), "\t")
    print("result rel_error_percent %.1f" % (rel_error_percent), "\t")
    '''

Ich programmiere erfolglos, also bin ich nicht.
Sirius3
User
Beiträge: 17741
Registriert: Sonntag 21. Oktober 2012, 17:20

Was ist denn `coords` für ein Typ. Die LCs sehen irgendwie falsch aus.
Wo hängt es denn nun noch? Gibt es einen Fehler? Wie weicht das erwartete vom eigentlichen Ergebnis ab?
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

@Sirius3: 'coords' ist vom Typen <class 'list'>. Was bedeutet "LC"s? In der Tat war die Funktionsschnittstelle falsch. Ich habe sie jetzt geändert:

Code: Alles auswählen

myResult = haversine(lat0, lats, lon0, lons)
Es hängt dahingehend:
TypeError: only size-1 arrays can be converted to Python scalars
Für ein Ergebnis ist es noch zu früh. Weder ist es mir bisher gelungen, über 2 numpy-arrays zu rechnen, noch das Dreiecksmatrixproblem zu lösen (Flughafen mit sich selbst, A:B = B:A)
Grüße
Strawk

So sieht coords aus:
[(-6.0816898345900015, 145.391998291), (-5.20707988739, 145.789001465), (-5.826789855957031, 144.29600524902344), (-6.569803, 146.725977), (-9.443380355834961, 147.22000122070312), (-3.5838301181800003, 143.66900634799998), (61.1604995728, -45.4259986877), (64.19090271, -51.678100585900005), (67.0122218992, -50.7116031647), (76.5311965942, -68.70320129390001), (65.66000366210939, -18.07270050048828), (65.28330230712889, -14.401399612426758), (64.295601, -15.2272), (65.952301, -17.426001), (66.05809783935547, -23.13529968261719), (63.985000610352, -22.605600357055998), (65.555801, -23.965), (64.1299972534, -21.9405994415), (66.133301, -18.9167), (63.42430114746094, -20.27890014648437), (46.48500061035156, -84.5093994140625), (50.0564002991, -97.0325012207), (44.639701843299996, -63.4994010925), (51.391899108900006, -56.083099365200006), (49.079833, -125.77558300000001), (68.534401, -89.808098), (49.13249969482422, -68.20439910888672), (48.33060073852539, -70.99639892578125), (64.29889678960001, -96.077796936), (49.950 ... ... ... ... ...
Ich programmiere erfolglos, also bin ich nicht.
Sirius3
User
Beiträge: 17741
Registriert: Sonntag 21. Oktober 2012, 17:20

LCs sind ListComprehensions und die kann man also einfach ersetzen:

Code: Alles auswählen

lats, lons = np.array(coords).T
Jetzt wäre es noch gut, zu wissen, wo genau der Fehler auftritt, denn dort wirst Du merken, dass Du eine Funktion durch die äquivalente Numpy-Variante ersetzen mußt.

`myResult` ist ein schlechter Name, oder gibt es auch ein your_result?
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

@Sirius3: Der Fehler tritt bei Zeile 17 auf, wo mit

Code: Alles auswählen

lat2 = radians(lat2)
ein Array, anstatt eines einzelnen Wertes übergeben wird. Wie lautet denn die äquivalente Numpy-Variante? Ich habe 'myResult' durch zwei erwartete Ergebnisse ersetzt: 'sumMeters' und 'nCalculations '. Das sind die Summe der Meter und die Anzahl der getätigten Distanzberechnungen.
Grüße, Strawk
Ich programmiere erfolglos, also bin ich nicht.
Benutzeravatar
__blackjack__
User
Beiträge: 13080
Registriert: Samstag 2. Juni 2018, 10:21
Wohnort: 127.0.0.1
Kontaktdaten:

Du musst das `radians()` aus `numpy` nehmen. Das aus `math` kennt logischerweise keine Numpy-Arrays.
„All religions are the same: religion is basically guilt, with different holidays.” — Cathy Ladman
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

Hallo ihrs! :)
Der Code sieht jetzt so aus:

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
# from math import radians
# from gpx_project import dist
from distances_airports_worldwide import get_airport_positions

def haversine(lat1, lat2, lon1, lon2):
    """ calculates the distance in kilometers between two geographical points
    """
    # nCalculations = len(lat2)
    
    lat1 = np.radians(lat1)
    lat2 = np.radians(lat2)
    lon1 = np.radians(lon1)
    lon2 = np.radians(lon2)
    
    leftParenthesis = ((lat2 - lat1)/2)
    rightParenthesis = ((lon2 - lon1)/2)

    allUnderSquareRoot = ((np.sin(leftParenthesis))**2) + np.cos(lat1)*np.cos(lat2)*(np.sin(rightParenthesis))**2
    squareRootResult = np.sqrt(allUnderSquareRoot)
    arcSin = np.arcsin(squareRootResult)
    # before: 6371 km
    middleEarthRadiusKm = 6378.1
    middleEarthRadiusM = middleEarthRadiusKm * 1000
    result = 2 * middleEarthRadiusM * arcSin

    return result #, nCalculations

def distvec(lat0, lats, lon0, lons):
    """calculates the distance of all airports worldwide to each other in kilometers
    """
    # lats = lats.astype(float)
    # lons = lons.astype(float)
    print(type(lats))
    print(type(lons))
    
    # lats = np.array(lats).astype(float)
    # lons = np.array(lons).astype(float)
    sumMeters = haversine(lat0, lats, lon0, lons)
    
    return sumMeters
    
    # lst = [2, 3, 7.9, 3.3, 6.9, 0.11, 10.3, 12.9]
    # v = np.array(lst)
    # print(v + 2)
    
    # return myResult
           
if __name__ == "__main__":
    coords = get_airport_positions("data\\airports.dat", "All")
    # print(coords)
    # print(type(coords))
    lats = [t[0] for t in coords]
    lons = [t[1] for t in coords]
    lats = np.array(lats)
    lons = np.array(lons)
    sumMeters = distvec(-6.081689834590001, lats, 145.391998291, lons)
    print(sumMeters)
    
    '''
    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 = haversine(lat1, lat2, lon1, lon2)
    print("km: %.1f" % (result), "\t")
    '''
    '''
    result_haversine = haversine(lat1, lat2, lon1, lon2) / 1000
    result_great_circle = dist((lat1, lon1), (lat2, lon2), method="great_circle") / 1000
    rel_error = abs((result_great_circle - result_haversine) / result_great_circle)
    rel_error_percent = rel_error * 100
    
    print("result haversine km: %.1f" % (result_haversine), "\t")
    print("result great_circle km: %.1f" % (result_great_circle), "\t")
    print("result rel_error km: %.7f" % (rel_error), "\t")
    print("result rel_error_percent %.1f" % (rel_error_percent), "\t")
    '''

und ich erhalte als Ergebnis:
runfile('C:/Users/Karl Kraft/Documents/Programmierung_ausser_PHP/Python/analysis_of_geo_data/geosphere_distances.py', wdir='C:/Users/Karl Kraft/Documents/Programmierung_ausser_PHP/Python/analysis_of_geo_data')
Reloaded modules: distances_airports_worldwide, gpx_project, extremesBoundBoxMap
<class 'numpy.ndarray'>
<class 'numpy.ndarray'>
[8.85139184e-11 1.06832824e+05 1.24619768e+05 ... 7.77830749e+06
5.70006415e+06 5.77467174e+06]
Ich programmiere erfolglos, also bin ich nicht.
Sirius3
User
Beiträge: 17741
Registriert: Sonntag 21. Oktober 2012, 17:20

Das Problem ist nur, dass da gar keine Summe aus Metern gebildet wird ;-). Namenskonvention ist in Python Variablen klein_mit_unterstrich zu schreiben.
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

Ja, genau. Eine einzelne Meter-Berechnung erfolgt. Zum Ziel siehe bitte meinen Beitrag von heute 17.56 Uhr. 2 Probleme bestehen noch: Wie sagt man der haversine-Funktion, dass sie eine Menge berechnen soll, nicht 1 Abstand. Und wie trägt man dabei dem Dreiecksmatrixproblem Rechnung? Alles also wie gehabt. Ideen?

Code aktuell:

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
# from math import radians
# from gpx_project import dist
from distances_airports_worldwide import get_airport_positions

def haversine(lat1, lat2, lon1, lon2):
    """ calculates the distance in kilometers between two geographical points
    """
    # nCalculations = len(lat2)
    
    lat1 = np.radians(lat1)
    lat2 = np.radians(lat2)
    lon1 = np.radians(lon1)
    lon2 = np.radians(lon2)
    
    leftParenthesis = ((lat2 - lat1)/2)
    rightParenthesis = ((lon2 - lon1)/2)

    allUnderSquareRoot = ((np.sin(leftParenthesis))**2) + np.cos(lat1)*np.cos(lat2)*(np.sin(rightParenthesis))**2
    squareRootResult = np.sqrt(allUnderSquareRoot)
    arcSin = np.arcsin(squareRootResult)
    # before: 6371 km
    middleEarthRadiusKm = 6378.1
    middleEarthRadiusM = middleEarthRadiusKm * 1000
    result = 2 * middleEarthRadiusM * arcSin
    result = sum(np.array(result))
    
    return result #, nCalculations

def distvec(lat0, lats, lon0, lons):
    """calculates the distance of all airports worldwide to each other in kilometers
    """
    # lats = lats.astype(float)
    # lons = lons.astype(float)
    # print(type(lats))
    # print(type(lons))
    
    # lats = np.array(lats).astype(float)
    # lons = np.array(lons).astype(float)
    sumMeters = haversine(lat0, lats, lon0, lons)
    
    return sumMeters
    
    # lst = [2, 3, 7.9, 3.3, 6.9, 0.11, 10.3, 12.9]
    # v = np.array(lst)
    # print(v + 2)
    
    # return myResult
           
if __name__ == "__main__":
    coords = get_airport_positions("data\\airports.dat", "All")
    # print(coords)
    # print(type(coords))
    lats = [t[0] for t in coords]
    lons = [t[1] for t in coords]
    lats = np.array(lats)
    lons = np.array(lons)
    sumMeters = distvec(-6.081689834590001, lats, 145.391998291, lons)
    print(sumMeters)
    
    '''
    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 = haversine(lat1, lat2, lon1, lon2)
    print("km: %.1f" % (result), "\t")
    '''
    '''
    result_haversine = haversine(lat1, lat2, lon1, lon2) / 1000
    result_great_circle = dist((lat1, lon1), (lat2, lon2), method="great_circle") / 1000
    rel_error = abs((result_great_circle - result_haversine) / result_great_circle)
    rel_error_percent = rel_error * 100
    
    print("result haversine km: %.1f" % (result_haversine), "\t")
    print("result great_circle km: %.1f" % (result_great_circle), "\t")
    print("result rel_error km: %.7f" % (rel_error), "\t")
    print("result rel_error_percent %.1f" % (rel_error_percent), "\t")
    '''

Ich programmiere erfolglos, also bin ich nicht.
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

Hallo!

Es hat sich was getan. Die Haversine-Funktion ist jetzt in der Lage, einen Vektor aufzunehmen und zu berechnen. Auf diese Weise erhält man alle Abstände von einem gegebenen Flughafen zu allen anderen Flughäfen (hier am Beispiel Deutschlands).

Was seltsam anmutet: Im Ergebnis kommen 3x viel zu hohe Werte vor: 8003.1 km, 928.8 km und 8215.7 km. So groß ist die Ausdehnung Deutschlands bekanntlich nicht. Ich hoffe, dass es sich hier um einen kleinen, einfach zu korrigierenden Fehler handelt.

Fernziel: Alle Abstände von allen Flughäfen zu allen berechnen. Das Dreiecksmatrixproblem kann hierbei erstmal außen vor bleiben. Hat jemand eine Idee, wie man bei gegebenen Koordinaten der Flughäfen alle Abstände aller Flughäfen zu allen in einem oder mehreren Vektoren ausdrücke kann. Die Berechnung wäre dann 100x schneller, als mit einer Doppelschleife.

Grüße
Strawk

Aktueller 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
from distances_airports_worldwide import get_airport_positions

def haversine(lat1, lat2, lon1, lon2):
    """ calculates the distance in meters between geographical points
    """
    # Haversine formula
    lat1 = np.radians(lat1)
    lat2 = np.radians(lat2)
    lon1 = np.radians(lon1)
    lon2 = np.radians(lon2)
    leftParenthesis = ((lat2 - lat1)/2)
    rightParenthesis = ((lon2 - lon1)/2)
    allUnderSquareRoot = ((np.sin(leftParenthesis))**2) + np.cos(lat1)*np.cos(lat2)*(np.sin(rightParenthesis))**2
    squareRootResult = np.sqrt(allUnderSquareRoot)
    arcSin = np.arcsin(squareRootResult)
    middleEarthRadiusKm = 6371.009
    middleEarthRadiusM = middleEarthRadiusKm * 1000
    distanceMeters = 2 * middleEarthRadiusM * arcSin
    
    return distanceMeters

def haversine_vectorized(coord0, coords):
    """ prepares data vor vectorized use of haversine-function
    """
    lat0, lon0 = coord0
    lats = coords[:,0] # 1. Spalte von coords (= array von latitudes)
    lons = coords[:,1] # 2. Spalte von coords (= array von longitudes)
    result_vectorized = haversine(lat0, lats, lon0, lons)
    
    return result_vectorized

if __name__ == "__main__":
    coords = get_airport_positions("data\\airports.dat", "Germany")
    coord0 = coords[0]
    coords = np.array(coords)
    coord0 = np.array(coord0)
    result_vectorized = haversine_vectorized(coord0, coords) / 1000
    np.set_printoptions(formatter={'float': '{: 0.1f}'.format})
    print(result_vectorized)

Ergebnis:
[ 0.0 142.6 176.1 111.5 68.3 180.4 175.7 81.8 263.9 159.7 146.4
370.2 198.5 239.4 222.4 35.4 299.1 77.4 29.6 148.7 52.9 249.9
439.0 482.7 409.6 161.6 517.0 539.6 371.8 308.4 160.8 571.1
470.2 173.9 360.8 442.2 436.7 530.5 467.0 407.9 470.2 470.0
283.0 454.0 410.9 460.3 387.9 568.1 482.6 460.9 527.5 421.6
558.1 412.0 537.6 480.8 401.2 484.3 336.7 447.6 417.1 292.3
299.9 493.0 533.2 310.2 244.5 405.7 243.7 285.6 212.7 310.4
500.9 565.9 479.7 552.4 561.1 596.7 509.0 408.6 70.8 300.1
357.7 333.8 478.6 552.5 534.2 508.5 593.1 567.2 525.9 472.2
584.1 568.8 528.5 365.4 393.3 345.5 505.8 365.1 468.0 514.6
380.8 327.8 352.9 246.9 407.7 291.0 345.6 489.1 451.1 592.5
482.1 519.7 338.7 551.2 494.6 533.5 377.0 261.0 466.0 455.0
433.9 537.2 368.8 406.3 113.3 348.6 425.5 363.6 437.0 471.2
584.2 532.6 582.6 565.9 576.2 299.6 501.5 555.3 363.6 476.1
203.6 425.3 417.3 76.4 532.8 362.3 532.8 531.6 553.3 8132.5
461.1 360.8 496.6 492.1 691.9 8003.1 928.8 8215.7 146.9 515.2
512.4 573.7 286.7 226.7 347.4 456.3 214.7 324.2 50.1 413.0
498.2 466.9 236.0 428.2 458.7 403.9 458.2 504.7 270.6 378.9
412.8 344.7 341.1 467.7 165.7 406.3 409.2 371.2 489.9 187.3
477.0 301.7 442.8 471.9 392.2 473.5 523.2 306.9 237.7 217.4
522.6 487.7 524.8 549.4 170.2 413.1 136.2 190.2 157.8 504.7
456.2 567.7 226.5 333.4 311.5 315.9 617.5 578.4 334.2 341.1
391.3 408.4 270.3 445.3 539.9 276.2 437.2 430.5 422.2 392.7
362.4 409.7 408.0 400.4 404.6 444.9 280.2 291.4 186.7]
Ich programmiere erfolglos, also bin ich nicht.
Sirius3
User
Beiträge: 17741
Registriert: Sonntag 21. Oktober 2012, 17:20

Du mußt halt die gesamte Matrix aller Kombinationen ausrechnen, also lat1/lon1 als Spalten- und lat2/lon2 als Zeilenvektor benutzen.
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

@Sirius3: Aber bedeutet das nicht, dass ich wieder zur Schleife zurückkehren muss? Bisher wurde erreicht, dass über einen ganzen Vektor berechnet werden kann.
Ich programmiere erfolglos, also bin ich nicht.
Sirius3
User
Beiträge: 17741
Registriert: Sonntag 21. Oktober 2012, 17:20

@Strawk: nein, keine Schleife, nur Matrizen statt Vektoren.
Benutzeravatar
Strawk
User
Beiträge: 229
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: Aachen
Kontaktdaten:

Hallo!

Code: Alles auswählen

print(len(coords))
for i in range(7184):
    coord0 = coords[i]
ergibt:
7184
IndexError: index 119 is out of bounds for axis 0 with size 44
Was ist da los?
Ich programmiere erfolglos, also bin ich nicht.
Antworten