Abstandsberechnung zwischen 2 Städten

[schöni]

Zitat:

Zitat von jaenicke

Routenplaner berücksichtigen ja die exakte Streckenführung. Deshalb ist da natürlich der Wert nicht Luftlinie sondern entspricht der zurückgelegten Strecke.

Allgemein kann man das ja nicht berechnen ohne die tatsächlichen Daten zur Streckenführung zu haben, wie es bei Routenplanern der Fall ist.

Schade. Ist also doch aufwendiger.

Zitat:

Zitat von sx2008

Ein Routenplaner ist ein recht komplexes Stück Software.
Falls dich das interessiert, kannst du ja mal OpenStreetMap anschauen.

Danke für den Link! Werd ich mir gleich mal ansehen.

[BoolString]

Für Koordinaten gibt es grundsätzlich 2 Verfahren. Das erste geht von einer idealen Kugelgestalt der Erde aus. Dies ist sicherlich innerhalb Deutschlands und für deine Belange ausreichend genau. Hierbei geht man meist von einem mittleren Planetenradius von 6371 km aus (die Angaben variieren je nach verwendeter Referenz).

Delphi-Quellcode:

Function CalcSimpleGeoDistance (Lat1, Lon1, Lat2, Lon2, aRadius : Double) : Double;
// Calculates the distance between two positions, based on given
// Latitude and Longitude values on an ideal sqhere with a given
// radius; Dr. Jan Schulz, 2007
Var aPart1 : Double;
    aPart2 : Double;
    aPart3 : Double;
Begin
  Lat1 := Lat1 * Pi / 180;
  Lon1 := Lon1 * Pi / 180;
  Lat2 := Lat2 * Pi / 180;
  Lon2 := Lon2 * Pi / 180;

  aPart1 := Cos (Lat1) * Cos (Lon1) * Cos (Lat2) * Cos (Lon2);
  aPart2 := Cos (Lat1) * Sin (Lon1) * Cos (Lat2) * Sin (Lon2);
  aPart3 := Sin (Lat1) * Sin (Lat2);

  Result := ArcCos (aPart1 + aPart2 + aPart3) * aRadius;
end;

Leider ist die Erde aber ein >68-eckiger Klumpen. Für genauere Betrachtungen sollte man deshalb berücksichtigen, daß die Erde durch ihre Rotation zu den Polen hin abgeflacht ist. Besonders wenn man über mehrere 1000 km Distanz mit Koordinaten ausrechnen will. Hier gibt es den Algorithmus von Andoyer, der ihn - um 1949 herum - in einem französischen Ephemeridenkalender veröffentlicht hat. Hab es nicht nachgemessen, aber die Abweichung soll für die Distanz Paris-New York im Bereich der quadrierten Abplattung liegen (<50 m). Natürlich nur in Bezug auf ein Referenzellipsoid. Tatsächlich Weglängen durch Abweichungen im Geoid sind nicht berücksichtigt (Abplattung, Referenzellipsoid und Geoid siehe Wikipedia)...

Delphi-Quellcode:

Function CalcAndoyersGeoDistance (Lat1, Lon1, Lat2, Lon2, RadiusEQ, RadiusPol : Double) : Double;
// Calculates a more precise distance between two positions, based on given
// Latitude and Longitude values for a planet with respect to the flattening
// calculated from the polar and equatorial distance; Dr. Jan Schulz, 2007
Var Flattening : Double;
    AndoyerF   : Double;
    AndoyerG   : Double;
    AndoyerL   : Double;
    AndoyerC   : Double;
    AndoyerS   : Double;
    AndoyerO   : Double;
    AndoyerR   : Double;
    H1         : Double;
    H2         : Double;
Begin
  Lat1 := Lat1 * Pi / 180;
  Lon1 := Lon1 * Pi / 180;
  Lat2 := Lat2 * Pi / 180;
  Lon2 := Lon2 * Pi / 180;

  Flattening := CalcGeometricFlattening (RadiusEQ, RadiusPol);

  AndoyerF := (Lat1 + Lat2) / 2;
  AndoyerG := (Lat1 - Lat2) / 2;
  AndoyerL := (Lon1 - Lon2) / 2;

  AndoyerC := Sqr (Cos (AndoyerG)) * Sqr (Cos (AndoyerL)) +
              Sqr (Sin (AndoyerF)) * Sqr (Sin (AndoyerL));

  AndoyerS := Sqr (Sin (AndoyerG)) * Sqr (Cos (AndoyerL)) +
              Sqr (Cos (AndoyerF)) * Sqr (Sin (AndoyerL));

  AndoyerO := ArcTan ((Sqrt (AndoyerS / AndoyerC)));

  AndoyerR := Sqrt (AndoyerS * AndoyerC) / AndoyerO;

  H1 := (3 * AndoyerR - 1) / (2 * AndoyerC);
  H2 := (3 * AndoyerR + 1) / (2 * AndoyerS);

  Result := 2 * AndoyerO * RadiusEQ * ( 1 +
            Flattening * H1 * Sqr (Sin (AndoyerF)) * Sqr (Cos (AndoyerG)) -
            Flattening * H2 * Sqr (Cos (AndoyerF)) * Sqr (Sin (AndoyerG)));
end;


Function CalcGeometricFlattening (RadiusEQ, RadiusPol : Double) : Double;
// Calculates the geometric flattening of a planet as relative difference
// between the equatorial and polar radius; Dr. Jan Schulz, 2007
Begin
  If RadiusEQ > 0 THen Result := (RadiusEQ - RadiusPol) / RadiusEQ
                  Else Result := 0;
end;

Bei dieser Methode nimmt man für den Äquatorialradius normalerweise 6378.137 km und für den Polradius 6356.752 km. Je nachdem welches Ellipsoid hierbei zu Grunde gelegt wird, weichen die einzelnen Werte voneinander ab.

Liebe Grüße aus ><)))°> Town

Jan