Haversine-Formel liefert falsche Werte

[Strawk]

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]

[Sirius3]

Du mußt halt die gesamte Matrix aller Kombinationen ausrechnen, also lat1/lon1 als Spalten- und lat2/lon2 als Zeilenvektor benutzen.

[Strawk]

@Sirius3: Aber bedeutet das nicht, dass ich wieder zur Schleife zurückkehren muss? Bisher wurde erreicht, dass über einen ganzen Vektor berechnet werden kann.

[Sirius3]

@Strawk: nein, keine Schleife, nur Matrizen statt Vektoren.

[Strawk]

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?

[__blackjack__]

@Strawk: Wo kommt denn die zweite Zeile her? Oder wo sind die Zeilen dazwischen geblieben?

[Strawk]

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.

[ThomasL]

mach mal print(coords.shape)
in der Annahme, dass das ein numpy ndarray ist.

[Strawk]

Hallo ThomasL!
Wäre ja nicht das erste Mal, dass du mir entscheidend hilfst.

(7184, 2)

Grüße
Strawk

[Strawk]

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.

[__blackjack__]

@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.

[ThomasL]

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?

[Strawk]

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))

[Sirius3]

@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))

[Strawk]

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

[Sirius3]

Bis auf den Punkt, dass ich lat und lon ein bißchen durcheinander gebracht habe, rechnet ›calc_dist_all2all_airports‹ genau das was du willst.

[Strawk]

Hallo!

@Sirius3: Habe ich denn etwas nicht gleichlautendes behauptet?

Grüße, Strawk