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: 160
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: RegBez Köln
Kontaktdaten:

Dienstag 7. Mai 2019, 18:28

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: 9896
Registriert: Sonntag 21. Oktober 2012, 17:20

Dienstag 7. Mai 2019, 21:10

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: 160
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: RegBez Köln
Kontaktdaten:

Donnerstag 9. Mai 2019, 11:50

@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: 9896
Registriert: Sonntag 21. Oktober 2012, 17:20

Donnerstag 9. Mai 2019, 12:12

@Strawk: nein, keine Schleife, nur Matrizen statt Vektoren.
Benutzeravatar
Strawk
User
Beiträge: 160
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: RegBez Köln
Kontaktdaten:

Donnerstag 9. Mai 2019, 12:48

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.
Benutzeravatar
__blackjack__
User
Beiträge: 3353
Registriert: Samstag 2. Juni 2018, 10:21

Donnerstag 9. Mai 2019, 12:50

@Strawk: Wo kommt denn die zweite Zeile her? Oder wo sind die Zeilen dazwischen geblieben?
“I am Dyslexic of Borg, Your Ass will be Laminated” -- unknown
Benutzeravatar
Strawk
User
Beiträge: 160
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: RegBez Köln
Kontaktdaten:

Donnerstag 9. Mai 2019, 13:12

Code: Alles auswählen

coords = get_airport_positions("data\\airports.dat", "All")
counter = 0
# print(len(coords))
for i in range(7184):
    coord0 = coords[i]
    coords = np.array(coords[i+1:]) # slicing von i+1 bis Ende Array
ergibt:
IndexError: index 119 is out of bounds for axis 0 with size 44
Mehr ist da nicht, alles andere ist auskommentiert.
Ich programmiere erfolglos, also bin ich nicht.
Benutzeravatar
ThomasL
User
Beiträge: 686
Registriert: Montag 14. Mai 2018, 14:44
Wohnort: Kreis Unna NRW

Donnerstag 9. Mai 2019, 13:16

mach mal print(coords.shape)
in der Annahme, dass das ein numpy ndarray ist.
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: 160
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: RegBez Köln
Kontaktdaten:

Donnerstag 9. Mai 2019, 13:21

Hallo ThomasL! :-)
Wäre ja nicht das erste Mal, dass du mir entscheidend hilfst.
(7184, 2)
Grüße
Strawk
Ich programmiere erfolglos, also bin ich nicht.
Benutzeravatar
Strawk
User
Beiträge: 160
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: RegBez Köln
Kontaktdaten:

Donnerstag 9. Mai 2019, 13:38

Code: Alles auswählen

coords = get_airport_positions("data\\airports.dat", "All")
coords = np.array(coords)
counter = 0
# print(len(coords))
print(coords.shape)
for i in range(10):
    coord0 = coords[i]
    coords = coords[i+1:] # slicing von i+1 bis Ende Array
    print(len(coords))
(7184, 2)
7183
7181
7178
7174
7169
7163
7156
7148
7139
7129
Verkürzt sich in immer größeren Schritten, statt 1er-Schritten. Ich habe wohl ein größeres Brett vorm Kopf.
Ich programmiere erfolglos, also bin ich nicht.
Benutzeravatar
__blackjack__
User
Beiträge: 3353
Registriert: Samstag 2. Juni 2018, 10:21

Donnerstag 9. Mai 2019, 14:01

@Strawk: Du veränderst/verkürzt die Länge von `coords` ja bei jedem Schleifendurchlauf *selbst* und wunderst Dich dann, dass die kürzer wird? Und genau das hattest Du beim ersten mal als Du das Problem beschrieben hast, auch gar nicht gezeigt, dass Du `coords` immer wieder an kürzere Arrays bindest.
“I am Dyslexic of Borg, Your Ass will be Laminated” -- unknown
Benutzeravatar
ThomasL
User
Beiträge: 686
Registriert: Montag 14. Mai 2018, 14:44
Wohnort: Kreis Unna NRW

Donnerstag 9. Mai 2019, 16:37

Länge coords = 7184
i läuft von 0 bis 9 also
i=0 coords = coords[1:] Länge um 1 gekürzt 7183
i=1 coords = coords[2;] Länge nochmals um 2 gekürzt also 7181
i=2 coords = coords[3;] Länge nochmals um 3 gekürzt also 7178
etc
Was willst du denn mit dieser Zuweisung bewirken?
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: 160
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: RegBez Köln
Kontaktdaten:

Donnerstag 9. Mai 2019, 20:37

Hi!

So! :)

Bitte kritisieren. Die Konvention Funktionen/Variablen kann ich nicht ändern; mein Lehrer will es so.

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
import time

def haversine(lat1, lat2, lon1, lon2):
    """ calculates the distance in meters between geographical points
    """
    # in: 2 pairs of coordinates
    # out: distance in meters of them
    
    # 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
    distancesMeters = 2 * middleEarthRadiusM * arcSin
    
    return distancesMeters

def haversine_vectorized(coord0, coords):
    """ prepares data vor vectorized use of haversine-function
    """
    # in: one tupel of first 2 coordinates, list of all remainung coordinates
    # out: vectorized result, which is list of kilometers
    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

def calc_dist_all2all_airports(airport_file, country):
    """ calculates the distances of each airport to each and adds them
    """
    # in: the file where the airport data are stored, the desired country
    # out: added kilometers of all the distances, number of calculations done
    coords = get_airport_positions(airport_file, country)
    coords = np.array(coords)
    coordsOneValue = coords
    multipleKilometers = []
    counter = 0
    for i in range(len(coords)):
        coord0 = coordsOneValue[i]
        coords = coords[i+1:] # slicing from i+1 until end array
        result_vectorized = haversine_vectorized(coord0, coords) / 1000
        multipleKilometers.extend(result_vectorized)
        counter += 1
    sumKilometers = sum(multipleKilometers)
    nCalculations = counter
    
    return sumKilometers, nCalculations

# number of airports: 7184
if __name__ == "__main__":
    t0 = time.clock()
    sumKilometers, nCalculations = calc_dist_all2all_airports("data\\airports.dat", country="All")
    print("sum kilometers: %10.1f" % (sumKilometers))
    print("number of calculations done: %10.f" % (nCalculations))
    print("Computation of distances took %6.2f secs" % (time.clock()-t0))

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

Donnerstag 9. Mai 2019, 21:34

@Strawk: Konvention hin oder her, Variablen schreibt man klein_mit_unterstrich. Warum nennst Du nCalculations mal counter? und zählst die Länge von coords, das ist doch unnötig. Und warum n? number_of_calculations (richtig geschrieben) sind das ja auch nicht, sondern nur die Anzahl der Flughäfen (number_of_airports).
Was soll dann noch coords_one_value? Das ist das selbe wie coords, also überflüssig.

Mit ein bißchen numy wird das ganze dann deutlich kürzer.

Code: Alles auswählen

EARTH_RADIUS_IN_M = 63710090

def haversine(lat1, lat2, lon1, lon2):
    """ calculates the distance in meters between geographical points
    in: 2 pairs of coordinates in radians
    out: distance in meters of them
    """
    mean_lat = (lat2 - lat1) * 0.5
    mean_lon = (lon2 - lon1) * 0.5
    a = np.sin(mean_lat)**2 + np.cos(lat1)*np.cos(lat2)* np.sin(mean_lon)**2
    return 2 * EARTH_RADIUS_IN_M * np.arcsin(a**0.5)

def calc_dist_all2all_airports(airport_file, country):
    """ calculates the distances of each airport to each and adds them
    in: the file where the airport data are stored, the desired country
    out: added kilometers of all the distances, number of calculations done
    """
    coords = get_airport_positions(airport_file, country)
    lat, lon = np.radians(coords).T
    distances = haversine(lat[:, None], lon[:, None], lat[None, :], lon[None, :])
    return distances.sum()/2, len(coords)

# number of airports: 7184
if __name__ == "__main__":
    t0 = time.clock()
    total_distance, number_of_airports = calc_dist_all2all_airports("data\\airports.dat", country="All")
    print("sum of kilometers: %10.1f" % total_distance)
    print("number of airports: %10d" % number_of_airports)
    print("Computation of distances took %6.2f secs" % (time.clock()-t0))
Benutzeravatar
Strawk
User
Beiträge: 160
Registriert: Mittwoch 15. Februar 2017, 11:42
Wohnort: RegBez Köln
Kontaktdaten:

Sonntag 12. Mai 2019, 15:25

Hallo!

Aufgabe ist aktuell, eine Matrix distance_airports(i)(j) zu schaffen, in der jeweils die Abstände in Kilometern des i-ten Flughafens zum j-ten Flughafen enthalten sind. Daran beiße ich mir seit Dienstag die Zähne aus. Problem ist, dass die Daten ja bisher mittels Vektor generiert werden. Das kriege ich mit einer Struktur (i)(j) nicht zusammen gebacken. Benötige dazu mal einen heißen Tipp. Ich habe runde Klammern benutzt, um die Forumssoftware nicht zu verwirren (italics).

Strawk :)
Ich programmiere erfolglos, also bin ich nicht.
Antworten