Distanz zwischen zwei Wegpunkten berechnen (Luftlinie)

[webmantz]

Hallo zusammen,

ich bin gerade dabei, eine Datenbank für das maps4free-Projekt zu entwickeln.
Hierbei stellt sich u.A. die Aufgabe, die Entfernung zwischen zwei angegebenen Koordinaten zu errechnen, wobei die Luftlinie ohne Berücksichtigung des Höhenprofils ausreicht. Die Koordinaten sind alle im WGS84 Datum.

Kennt jemand eine Formel, wie ich aus zwei Koordinaten die Entfernung berechnen kann?

Danke für eure Hilfe.

Gruß, André

[RaGo]

Hallo André,
hier die Formeln, habe ich mir von einer Seite kopiert, weiss leider nicht mehr die Webadresse:
Entfernungsberechnung

Fuer unsere Zwecke hier ist die gaengige Formel zwischen A und B wohl

arccos(sin(B_lat)*sin(A_lat)+cos(B_lat)*cos(A_lat)*cos(B_lon - A_lon))

... multipliziert mit dem Erdradius. Aber wie gross ist der? Unsere Erde ist ja keine perfekte Kugel mit z.B. Radius 6371,110 km oder Umfang = 40000 km sondern eher etwas flachgedrueckt. Unterschiedliche Quellen nennen zwischen 6336 und 6399 km. Fuer unsere Breitengrade sollte 6380 km ganz gut hinkommen.


1) Pythagoras

d = sqrt((X2 - X1)^2 + (Y2 - Y1)^2)

kartesische Koordinaten fuer X/Y; Genauigkeit angeblich ausreichend bis ca. 20 km: Breitengrad / Meter Abweichung

70 / 30
50 / 20
30 / 9

2) Haversine Formel

dlon = B_lon - A_lon
dlat = B_lat - A_lat
a = (sin(dlat/2))^2 + cos(A_lat) * cos(B_lat) * (sin(dlon/2))^2
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c

R = Erdradius

Bei zwei Punkte auf gegenueberliegenden Kugelseiten betraegt der Fehler etwa 2 km, verursacht durch Rundungsfehler in der Sinus-Berechnung bei kleinen Winkeln. Ansonsten ist die Formel gut geeignet. Fehler entstehen insgesamt aus der rechnerinternen Rundung der Zahlenwerte, wie auch aus der begrenzten Laenge von Pi.

Fuer weit auseinander liegende Punkte bietet sich z.B. eine modifizierte Variante an:

dlon = B_lon - A_lon
dlat = B_lat - A_lat
a = (sin(dlat/2))^2 + cos(A_lat) * cos(B_lat) * (sin(dlon/2))^2
c = 2 * atan2( sqrt(a), sqrt(1-a) )
d = R * c

3) nahe der Polarkappen

a = pi/2 - A_lat
b = pi/2 - B_lat
c = sqrt( a^2 + b^2 - 2 * a * b * cos(B_lon - A_lon) )
d = R * c

4) arccos

... diese Formel ist mathematisch korrekt, aber ungeeignet fuer geringe Entfernungen!

a = sin(A_lat) * sin(B_lat)
b = cos(A_lat) * cos(B_lat) * cos(B_lon - A_lon)
c = arccos(a + b)
d = R * c

cos (5 degree) = 0.996194698
cos (1 degree) = 0.999847695
cos (1 minute) = 0.9999999577
cos (1 second) = 0.9999999999882
cos (0.05 sec) = 0.999999999999971

Bei nur sieben signifikanten Stellen reicht die Rechenaufloesung bis max. 1 Minute. Tatsaechlich duerften wir mit z.B. 15 signifikanten Stellen ganz gut im Bereich der maximalen Aufloesung von Sekunden liegen. Mit dem Taschenrechner ist da aber nicht mehr viel zu wollen.

All das relativiert sich mit der Genauigkeit des genutzten Erdumfangs. Ich wusste bisher nicht, dass die Nautische Meile / Seemeile mit 1,852 km als Minute vom Aequatorumfang ableitet.

Der tatsaechliche Radius ist wohl ganz gut angenaehert durch die Formel

R' = a * (1 - e^2) / (1 - e^2 * (sin(lat))^2)^(3/2)

a: Erdumfang
b: Polarradius
e: Exxentrizitaet sqrt(1 - b^2/a^2)

bzw.

a = 6378 km (3963 mi) Equatorial radius (surface to center distance)
b = 6357 km (3950 mi) Polar radius (surface to center distance)
e = 0.081082 Eccentricity

gives the following table of values for the

Radii of Curvature:

latitude...........r...................R'..................N..........

00 degrees . 6357 km (3950 mi) . 6336 km (3937 mi) . 6378 km (3963 mi)

15 degrees . 6360 km (3952 mi) . 6340 km (3940 mi) . 6379 km (3964 mi)

30 degrees . 6367 km (3957 mi) . 6352 km (3947 mi) . 6383 km (3966 mi)

45 degrees . 6378 km (3963 mi) . 6367 km (3957 mi) . 6389 km (3970 mi)

60 degrees . 6388 km (3970 mi) . 6383 km (3966 mi) . 6394 km (3973 mi)

75 degrees . 6396 km (3974 mi) . 6395 km (3974 mi) . 6398 km (3975 mi)

90 degrees . 6399 km (3976 mi) . 6399 km (3976 mi) . 6399 km (3976 mi)

.... aber interessant ist das wohl eher fuer Entfernungen auf See.

Denn bei uns auf dem Festland wird das 'Normal-Null' ohnehin willkuerlich festgelegt und liegt z.B. um 20 cm tiefer(?) als in der Schweiz.


Rainer

[webmantz]

Moin Rainer,

erst mal danke für die ausführlichen Formeln.

Allerdings scheint da was nicht zu stimmen. ich habe mal wahllos 2 Strecken rausgesucht und diese berechnen lassen:

erster Versuch:
9,4782611 Ost 52,675321 N => 9,48305922 O 52,6626576 N
Entfernung laut Magic Maps 3D: 1,45 km
Entfernung laut o.g. Formel: 83,906 km

zweiter Versuch:
10,38356389 Ost 53,27534084 N => 8,98390816 O 53,04155291 N
Entfernung laut Magic Maps 3D: 97,2 km
Entfernung laut o.g. Formel: 8706,866901 km

Ich habe die Berechung in Excel vorgenommen. Die Excel-Formel lautet:

Code: Alles auswählen

=ARCCOS( SIN(B_LAT)*SIN(A_LAT) + COS(B_LAT)*COS(A_LAT)*COS(B_LON-A_LON) ) * 6380

wobei A_LAT, A_LON, B_LAT, B_LON benannte Zellen mit den entsprechenden Werten sind.

Was mache ich falsch?

Gruß, André

[Jacob]

Nimm doch eine Strecke, deren genaue Länge du kennst, z.B. zwischen deinem Computer und dem Kühlschrank und spiele solange mit der Formul und eventuellen Fehlern rum, bis diss Ergebins hinkommt...

[webmantz]

Hallo Jacob,

Nimm doch eine Strecke, deren genaue Länge du kennst, z.B. zwischen deinem Computer und dem Kühlschrank und spiele solange mit der Formul und eventuellen Fehlern rum, bis diss Ergebins hinkommt...

danke für diesen sehr wissenschaftlichen Vorschlag *kopfschüttel*

Gruß, André

[macnetz]

Ich habe die Berechung in Excel vorgenommen....
Was mache ich falsch?

Hallo André,

soweit ich weiss, klappt das in Excel nicht mit der direkten Grad-Eingabe bei Winkelfunktionen. Du musst AFAIK stattdessen den Kreis-Bogen der Gradangabe nehmen . . .

Grüsse - Anton

[Roland]

Hallo Anton!

5 Minuten zu früh ...
Andrè: Wandle Deine Altgrad mit PI/180° um.

Grüße Roland

[webmantz]

Moin,

danke! Genau das wars
Ich hätte ja auch mal einen Blick in die Excelhilfe werfen können, aber das kommt davon wenn man glaubt, ein Programm gut genug zu kennen.

Gruß, André

[Horst Frank]

Tach!

Leider habe ich das aber noch nicht verstanden.
Was muss ich mit den werten machen, mit PI/180 multiplizieren?

[Roland]

Hallo Horst,

ja Winkeleinheiten machen auch unseren Lehrlingen Probleme- mir natürlich nie-nie-nie ... Auch auf den Taschenrechnern muss man einstellen ob man in Bogenmass, Altgrad oder Neugrad rechnen will.
Einmal hat der gestreckte Winkel 3,14159, einmal 180° oder (Vermessung) 200 gon.

Excel braucht z.B. 'Bogenmass'.
D.h. ein sin(90°) muss umgerechnet werden in sin(90°/180°xPI)= sin (PI/2).
Probier's aus, da muss 1 rauskommen (oder 0,9999999999 ... )

Grüße Roland

[Horst Frank]

Ich hatte das problem grade in einem perlscript .... nun stimmt alles ...

DANKE!!!

[macci]

Hallo zusammen
hier noch ein C-script:
-------------------------------------------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
float PI=3.14159;


// Calulate the distance between to co-ordinates in degrees
float check_distance(float lat1, float lon1, float lat2, float lon2)
{
// returns a value in kilometres
return 1.852 * 2 * asin( sqrt( pow(sin((lat1/180*PI-lat2/180*PI)/2),2) + cos(lat1/180*PI)*cos(lat2/180*PI)* pow(sin((lon1/180*PI-lon2/180*PI)/2),2)))*180*60/PI;
}


// Calulate the distance between to co-ordinates in degrees
float check_distance1(float lat1, float lon1, float lat2, float lon2)
{
// returns a value in kilometres
return 100*(acos(sin(lat2)*sin(lat1)+cos(lat2)*cos(lat1)*cos(lon2 - lon1)));
}


// Calulate the distance between to co-ordinates in degrees
float check_distance2(float lat1, float lon1, float lat2, float lon2)
{
// returns a value in kilometres
return 100*(sqrt(((lat2 - lat1)*(lat2 - lat1)) + ((lon2 - lon1)*(lon2 - lon1))));
}

int main (int argc, const char * argv[]) {

float lat1=52.675321;
float lon1=9.4782611;
float lat2=52.6626576;
float lon2=9.48305922;

printf("distance: %f \n",check_distance(lat1,lon1,lat2,lon2));
printf("distance1: %f \n",check_distance1(lat1,lon1,lat2,lon2));
printf("distance2: %f \n",check_distance2(lat1,lon1,lat2,lon2));

//controll it on http://gpsvisualizer.com/calculators

return 0;
}

[Roland]

Hallo macci,

das hat aber gedauert

Meinst Du, dass fünf Nachkommastellen für Pi ausreichen. Theoretisch sind das Meter, aber wenn Differenzen zwischen nahezu gleichgroßen Winkeln gebildet werden ?
Ich nehme immer mehr mit.

Grüße Roland

[macci]

Naja, bei mir hat's so gute Resultate gegeben.
Am genauesten war "distance" dann "distance1" und am ungenauesten war "distance2". Für mich und meinen 8-bit µC aber auch noch gut genug.

Gruss Martin

[krabat]

Die obige Formel geht davon aus, dass die Erde eine perfekte Kugel ist, was natürlich nicht ganz stimmt. Immerhin ist der Äquatorialradius ca. 20 km größer als der Polare.

Für noch genauere Ergebnisse gibt es den Algorithmus von Vincenty, der vom WGS-84-Ellipsoid ausgeht. Eine geschlossene Formel ist hier nicht möglich, aber der Algorithmus konvergiert sehr schnell.

http://www.movable-type.co.uk/scripts/latlong-vincenty.html (mit Javascript Code)

Zur Person: http://www.gfy.ku.dk/~iag/newslett/news7606.htm oder natürlich Wikipedia.