Differenz in Bogenmaß aus Entfernung in Meter

[Steffen81]

Gegeben sei eine Entfernung in Meter auf der Erde. Ich hätte gern daraus den maximalen Differenzbetrag in Bogenmaß (WGS-84 Ellipsiod). Umgekehrt habe ich eine Funktion (die krieg ich grad nur nicht mal so eben invertiert). Das Ergebnis muss nicht exakt sein (im Zweifel kann die Bogenmaßdifferenz ein wenig zu groß sein). Kennt da zufällig jemand eine Lösung?

[Christopher Spies]

Umgekehrt habe ich eine Funktion (die krieg ich grad nur nicht mal so eben invertiert).

Das liegt vermutlich daran, dass es einen Unterschied macht, wie die Strecke liegt (z.B. zwischen Nord-Süd- und Ost-West-Richtung).
Wenn Du nur eine obere Schranke für die Winkeldifferenz zwischen zwei Orten mit gegebener Entfernung voneinander brauchst, könnte folgendes funktionieren...

ACHTUNG: Das nachfolgend Geschriebene ist ohne Gewähr!

WGS84 ist ein Rotationsellipsoid, dessen kleinere Halbachse 6356,752315 km lang ist (Quelle: Wikipedia).

Im Extremfall befinden wir uns am Nordpol (oder besser gesagt genau in der Rotationsachse der Erde), dort hat der Ellipsoid den minimalen Radius.
Um eine Entfernung in eine Winkeldifferenz umzurechnen, dividiert man diese durch den Radius des Ellipsoids. Da der Radius an keinem Punkt kleiner ist als besagte 6356,75 km, kann die Winkeldifferenz den Kehrwert des Radius (1,57313e-7 rad/m) nicht überschreiten.
Die untere Grenze für die Winkeldifferenz ist dann durch den Kehrwert des größten Radius gegeben (6378,137 km => 1,56786e-7 rad/m). Dieser liegt am Äquator vor.

Eine Strecke von 1000 m schließt also nirgendwo auf der Erde einen Winkel von mehr als 1,57313e-7 rad/m * 1000 m = 1,57313e-4 rad = 0,00901337 ° = 32,4481 ",
und nirgendwo auf der Erde einen Winkel von weniger als 1,56786e-7 rad/m * 1000 m = 1,56786e-4 rad = 0,00898315 ° = 32,3394 " ein.

Wenn Deine Aufgabenstellung nur eine bestimmte Region der Erde (z.B. Deutschland) betrifft, könntest Du durch Bildung der partiellen Ableitungen Deiner Umrechenformel nach den Koordinaten in einem von Dir gewählten Arbeitspunkt linearisieren. Du erhältst dann eine bilineare Gleichung (statt wie hier eine lineare) mit größerer Genauigkeit.

ACHTUNG: Das vorstehend Geschriebene ist ohne Gewähr!

Man korrigiere mich, wenn ich falschliege.

Gruß,
- Christopher

[Steffen81]

Wenn Du nur eine obere Schranke für die Winkeldifferenz zwischen zwei Orten mit gegebener Entfernung voneinander brauchst, könnte folgendes funktionieren...

Genau dies suche ich
...aus Performancegründen. Eigentlich möchte ich alle Punkte die innerhalb einer bestimmten Distanz zu einem Referenzpunkt liegen. Allerdings wird dies bei sehr vielen Punkten sehr rechenintensiv. Daher möchte ich die zu vergleichenden Punkte vorher schon eingrenzen.

Ich habe es jetzt durch probieren abgeschätzt. Jedoch traue ich dem Ergebnis irgendwie nicht recht über den Weg. Wie du schon sagtest, hängt die Bogenmaßdistanz ja vom Ort auf der Erde ab. So ergab durch probieren:

Code: Alles auswählen

A (lat: 0.000, long: 0.000)
B (lat: 0.005, long: 0.000)
d = 554m

A (lat: 89.905, long: 0.000)
B (lat: 89.900, long: 0.000)
d = 558m

Müsste der Unterschied zwischen Äquator und Nordpol nicht viel größer sein?

Die Entfernung wird nach Vincenty berechnet.

[Christopher Spies]

Hallo Steffen,

Ich habe es jetzt durch probieren abgeschätzt. Jedoch traue ich dem Ergebnis irgendwie nicht recht über den Weg. Wie du schon sagtest, hängt die Bogenmaßdistanz ja vom Ort auf der Erde ab. So ergab durch probieren:

Code: Alles auswählen

A (lat: 0.000, long: 0.000)
B (lat: 0.005, long: 0.000)
d = 554m

A (lat: 89.905, long: 0.000)
B (lat: 89.900, long: 0.000)
d = 558m

Das passt doch gut zu den von mir angegebenen Umrechenfaktoren.
Eine Strecke von 554 m schließt immer einen Winkel von 0.00499340 ° oder weniger ein
und eine Strecke von 558m einen Winkel von 0.00501260 ° oder mehr.
Das stimmt doch genau mit denen von Dir durch "ausprobieren" gefundenen Werten überein (im Rahmen der Messgenauigkeit).

Müsste der Unterschied zwischen Äquator und Nordpol nicht viel größer sein?

Nein, warum?
So krass verbiegt sich ein Felsklumpen jetzt auch wieder nicht, wenn er mit einer Umdrehung pro Tag rotiert...
Und 22 km (Höhendifferenz, wohlgemerkt!) sind eigentlich schon ganz schön viel!

Gruß,
- Christopher

[Steffen81]

Das passt doch gut zu den von mir angegebenen Umrechenfaktoren...

Das find ich ja schon mal gut - spricht ja für deine Lösung...

Müsste der Unterschied zwischen Äquator und Nordpol nicht viel größer sein?

Nein, warum?
So krass verbiegt sich ein Felsklumpen jetzt auch wieder nicht, wenn er mit einer Umdrehung pro Tag rotiert...
Und 22 km (Höhendifferenz, wohlgemerkt!) sind eigentlich schon ganz schön viel!

Ok, das mag noch stimmen für die Breitengraddifferenz. Für die Längengraddifferenz habe ich allerdings die falschen Beispiele getestet.
Am Äquator entspricht 1/200 Längengraddifferenz (Ost-West) ca. 557m. Am Nordpol geht die Distanz beim gleichen Winkel gegen Null.

Code: Alles auswählen

A (lat: 89.9999, long: 0.0001)
B (lat: 89.9999, long: 0.0051)
d(NP) = 0.001m

A (lat: 0.0001, long: 0.0001)
B (lat: 0.0001, long: 0.0051)
d(Äq) = 557m

Das ist auch das, was ich erwartet hatte und daher traute ich meiner ersten Lösung nicht.

Deine Lösung beschreibt immer den einen Winkel zwischen zwei Punkten. Mir liegen aber für jeden Punkt jeweils nur Längen- und Breitengrad vor und das sind auch die Werte, auf die ich aus Performancegründen eingrenzen möchte. Wie gesagt, die exakte Distanzberechnung ist kein Problem.

Ich möchte also praktisch auf ein "Quadrat" (ist natürlich kein Echtes) auf der Erde eingrenzen mit der Kantenlänge 2*d (d = Distanz, nicht Durchmesser). Und dieses "Quadrat" wäre in Nordeuropa schon etwas kleiner als in Südeuropa. Allein Deutschland reicht mir nicht, wenn auch der Nordpol nicht ganz so wichtig ist.

Man könnte vom Referenzpunkt um die Distanz nach Norden (bzw. in Richtung des nächsten Pols) gehen und für diese Breite die Längengraddifferenz berechnen. Die Breitengraddifferenz haben wir ja schon, aber da die sich kaum unterscheidet muss man hier für die obere Schranke auch nicht wirklich rechnen.

[Christopher Spies]

Ok, das mag noch stimmen für die Breitengraddifferenz. Für die Längengraddifferenz habe ich allerdings die falschen Beispiele getestet. [...] Am Nordpol geht die Distanz beim gleichen Winkel gegen Null. [...] Das ist auch das, was ich erwartet hatte und daher traute ich meiner ersten Lösung nicht.

Ich habe natürlich nur entlang der Meridiane gedacht . Die an den Polen verschwindende Länge der Breitengrade habe ich ignoriert (peinlich!).

Ich möchte also praktisch auf ein "Quadrat" (ist natürlich kein Echtes) auf der Erde eingrenzen mit der Kantenlänge 2*d (d = Distanz, nicht Durchmesser). Und dieses "Quadrat" wäre in Nordeuropa schon etwas kleiner als in Südeuropa. Allein Deutschland reicht mir nicht, wenn auch der Nordpol nicht ganz so wichtig ist.

Mir ist immer noch nicht ganz klar, was Du brauchst.
Brauchst Du nur eine obere Schranke oder eine schnelle Näherung, um schnell einige in Frage kommende Punkte ausfindig zu machen (wobei einige falsch positive dabei sein dürfen und diese mit der exakten Berechnung ausgesiebt werden), oder brauchst Du eine "exakte" Formel?

Wie groß soll Dein "Quadrat" denn sein?
Wenn es klein genug ist, kann man im Arbeitspunkt linearisieren, sonst wird es kompliziert.

Gruß,
- Christopher

[Steffen81]

Ich versteh im Moment nicht ganz was du mit schneller Näherung meinst. Ich versuchs nochmal so:

Es sei ein Punkt Z und eine Distanz dmax (Radius) gegeben:
Z(latz [°];longz [°])
dmax [m] (z.B. 500m)

Zu Z möchte ich nun alle Punkte, die innerhalb einer Distanz dmax zu Z liegen. dmax ist für jedes Z konstant, kann aber global z.B. bei der nächsten Berechnung variieren.

Der einfachste Algorithmus ist, jeden Punkt P außer Z zu nehmen, nach Vincenty die Distanz d(ZPn) zu berechnen und wenn diese kleiner oder gleich dmax ist, den Punkt Pn zu merken. Nun ist die Menge der Punkte sehr groß und die Berechnung dauert entsprechend lange. Erste Maßnahme war, die Bestimmtung von d(ZPn) nicht nach Vincenty, sondern z.B. nach Haversine durchzuführen. Jedoch wäre es deutlich effizienter, einfach die Anzahl der zu vergleichenden Punkte zu reduzieren.

Wenn ich also nur diejenigen Punkte zum Vergleich nach obigem Algorithmus heranziehe, die in einem "Quadrat" mit Mittelpunkt Z und Kantenlänge 2*dmax liegen, dürfte das ganze recht schnell gehen. dmax liegt mir in Meter vor. Ich brauche diese aber in Grad, und zwar einmal als Längengraddifferenz und einmal als Breitengraddifferenz, sodass ich sagen kann:
Alle Punkte mit (latz - dlat) <= latz <= (latz + dlat) und (longz - dlong) <= longz <= (longz + dlong) sind vorläufig relevant. Und nur auf diese wird dann der obige Algorithmus angewandt.

Und eben die Werte dlat und dlong, die ich suche, sind abhängig von Z(latz;longz). Sie müssen nicht exakt sein, lieber etwas zu groß als zu klein. Aber Konstanten für Pi*Daumen-in-Deutschland reichen mir eben auch nicht ganz.

[Roland Ziegler]

Wenn man die Inverse einer Funktion nicht auf einfachem Wege ableiten kann, bietet sich immer noch die Methode des systematischen Probierens, sprich eine Iteration. Wesentliche Voraussetzung für solch einen iterativen Ansatz: Die Funktion muss konvergieren. Die einfachste Lösung, weil unabhängig von der Funktion, nennt man Bisektion. Die erreichbare Genauigkeit entspricht der der Originalfunktion. Die eigentliche Aufgabe des Anwenders besteht darin, vernünftige Startwerte für die Iteration zu bestimmen, damit diese mit ausreichender Wahrscheinlichkeit konvergiert. Zur Bestimmung der Startwerte kann eine sehr einfache und fast banale Näherungsfunktion schon ausreichen.

Bisektion wird übrigens im AR recht häufig angewandt, hat man es doch bei Klothoide und Artverwandten schnell mit nicht geschlossen darstellbaren Funktionen zu tun.



Nachtrag:
Der Ingenieur denkt ja bei Mathematik häufig pragmatisch. Welche Aufgabe ist zu lösen, und wie geht man da geschickt ran.

Du schreibst, Du willst wissen, welche Punkte zu einem Referenzpunkt innerhalb einer kreisförmigen Hüllkurve liegen. Und dies unangenehmerweise auf eine Ellipsoidenoberfläche. Warum verwendest Du nicht eine Kartenprojektion und rechnest mit zweidimensionaler Geometrie? Oder ist der Durchmesser Deines Kreises extrem groß? Wenn nicht: Ist die Kreisform ein Muss? Oder darf es auch ein ungefährer Kreis sein? Ist die absolute Entfernung wichtig? Daraus ergeben sich geeignete Projektionen, mit Winkel- oder Längen-/Flächentreue und wohlbekannten fertigen Formeln für Projektion und Inverse.

[Christopher Spies]

Hallo,

Ich versteh im Moment nicht ganz was du mit schneller Näherung meinst. [...] Sie müssen nicht exakt sein, lieber etwas zu groß als zu klein.

Genau das meinte ich mit "schnelle Näherung": Eine Faustformel, die Dir aus einer großen Menge von Punkten mit wenig Rechenaufwand mögliche Kandidaten aussucht. Es dürfen ein paar "falsch positive" dabei sein, diese werden dann mit der Formel nach Vincenty eliminiert.

Du schreibst, Du willst wissen, welche Punkte zu einem Referenzpunkt innerhalb einer kreisförmigen Hüllkurve liegen. Und dies unangenehmerweise auf eine Ellipsoidenoberfläche. [...] ist der Durchmesser Deines Kreises extrem groß? Wenn nicht: Ist die Kreisform ein Muss? Oder darf es auch ein ungefährer Kreis sein? Ist die absolute Entfernung wichtig?

Steffen sucht wohl die Punkte, die innerhalb eines Quadrats liegen, dessen Mittelpunkt der Referenzpunkt ist und dessen Kanten parallel zu Äquator bzw. Meridian verlaufen. Es ist also wohl eine winkeltreue Projektion gewünscht. Falsch positive Ergebnisse (Punkte, die sich eigentlich nicht im "Quadrat" befinden, aber dennoch dort lokalisiert werden) sind OK, falsch negative Ergebnisse (das Gegenteil) hingegen unerwünscht.

Ich versuche mich mal an einer Lösung:

ACHTUNG: Das nachfolgend Geschriebene ist ohne Gewähr!

WGS84 ist ein Rotationsellipsoid, dessen Halbachsen

Code: Alles auswählen

a = 6378,137 km
b = 6356,752315 km

lang sind (Quelle: Wikipedia).
Der Punkt Z mit dem Breitengrad latz (0° = Äquator, +-90° = Pole) und dem Längengrad longz hat vom Mittelpunkt des Ellipsoids den Abstand

Code: Alles auswählen

R_Lat := a * b / sqrt(a^2 * cos(lat_z)^2 + b^2 * sin(lat_z)^2)

(vor Anwendung dieser Formel muss man Längen- und Breitengrad gegebenenfalls ins Bogenmaß umrechnen).
Dieser Abstand stellt den Radius des Längenkreises dar, auf dem Z liegt.
Der Radius des Breitenkreises ist dann

Code: Alles auswählen

R_Lon := R_Lat * cos(lat_z)

.
Der Kehrwert dieser Radien ist der Umrechenfaktor zwischen Distanz und Winkeldifferenz.
Da der Ellipsoid rotationssymmetrisch ist, hängen die Radien nicht vom Längengrad ab.

Beispiel: Jena (52,9272° Nord) liegt gemäß obiger Formeln etwa 6364532 m vom Erdmittelpunkt entfernt. Der Radius des Breitenkreises beträgt 3836726 m. Die Umrechenfaktoren sind dann 0,0324085 "/m entlang des Längenkreises und 0,0537606 "/m entlang des Breitenkreises, eine Quadrat der Kantenlänge 1000m würde also 32,4085 " in Nord-Süd-Richtung und 53,7606 " in Ost-West-Richtung einschließen.
Am Äquator beträgt der Radius des Längenkreises ebenso wie der des Breitenkreises 6378137 m. Ein Quadrat der Kantenlänge 1000 m würde also eine Winkeldifferenz von 32,3394 " (entlang beider Achsen) einschließen.
Bei 89,9° Breite beträgt der Radius des Längenkreises 6356752 m und der des Breitenkreises 11095 m. Ein Quadrat der Kantenlänge würde also 32,4481 " in Nord-Süd-Richtung, aber 5,164 ° in Ost-West-Richtung einschließen.

ACHTUNG: Das vorstehend Geschriebene ist ohne Gewähr!

Man korrigiere mich wiederum, wenn ich falschliege.

Gruß,
- Christopher

[Steffen81]

...die Methode des systematischen Probierens, sprich eine Iteration. Wesentliche Voraussetzung für solch einen iterativen Ansatz: Die Funktion muss konvergieren.

Mein Problem ist weniger, ob sie konvergiert, sondern wann. Mein eigentliches Problem ist ja gelöst, es dauert mir nur zu lange.

Der Ingenieur denkt ja bei Mathematik häufig pragmatisch. Welche Aufgabe ist zu lösen, und wie geht man da geschickt ran.

Genau darum geht's mir.

Die Koordinaten liegen nun mal so (auf dem Ellipsoiden) vor und eine Projektion bedeutete auch wieder (Rechen)aufwand.
Der Durchmesser des (Um)Kreises ist typischerweise nicht sehr groß, daher reichen mir auch ungefähre Winkeldifferenzen für Länge und Breite (wie gesagt, lieber etwas zu groß als zu klein).

Ich habe eine einfach Formel gefunden, die mir die Winkeldifferenz zumindest nährungsweise ermittelt - und das sehr einfach. Da hätte ich auch gleich drauf kommen können.

Ausgehen von der Erde als Kugel sind alle Breitengrade ca. 111,12km voneinander entfernt. Die Längengrade sind cos(Breitengrad) * 111,12km voneinander entfernt. Daher ergibt sich (nährungsweise):

dlat = dmax / 111120.0
dlong = dmax / cos(latz) * 111120.0

[Christopher Spies]

Der Durchmesser des (Um)Kreises ist typischerweise nicht sehr groß, daher reichen mir auch ungefähre Winkeldifferenzen für Länge und Breite (wie gesagt, lieber etwas zu groß als zu klein).

Ich habe eine einfach Formel gefunden, die mir die Winkeldifferenz zumindest nährungsweise ermittelt - und das sehr einfach. [...] Ausgehen von der Erde als Kugel sind alle Breitengrade ca. 111,12km voneinander entfernt.

So kann man´s auch machen -- das ist sogar noch einfacher als meine Lösung .

Aber nimm für den Abstand der Breitengrade lieber 110,946 km an (ergibt sich aus dem minimalen Radius von 6356,75 km). Dann kannst Du sicher nein, dass die von Dir geschätzte Winkeldifferenz immer etwas zu groß ist. Ansonsten unterschätzt Du die Winkeldifferenz in Polnähe möglicherweise.

Gruß,
- Christopher

[Roland Ziegler]

Aus reiner Neugierde würde mich interessieren, welche Datenmengen da innerhalb welcher Antwortzeiten auf was für einer Implementierungsplattform verarbeitet werden sollen.

[Steffen81]

So kann man´s auch machen -- das ist sogar noch einfacher als meine Lösung

Richtig. Wie gesagt, hätte man auch gleich drauf kommen können.

Aber nimm für den Abstand der Breitengrade lieber 110,946 km an (ergibt sich aus dem minimalen Radius von 6356,75 km). Dann kannst Du sicher nein, dass die von Dir geschätzte Winkeldifferenz immer etwas zu groß ist. Ansonsten unterschätzt Du die Winkeldifferenz in Polnähe möglicherweise.

Guter Einwand. Damit ist sogar der Ellipsoid berücksichtigt.
Korrekterweise sollte es also folgendermaßen lauten:

dlat = dmax / 110946.0
dlong = dmax / cos(rad(latz+dLat)) * 111120.0

@Roland:
Es handelt sich im Moment um folgende Dimensionen - wenn alles funktioniert, können diese sich jedoch noch multiplizieren:
Punkte: < 10.000 in einer Datenbank (im Moment noch MySQL)
Die Anwendung läuft im Moment über PHP5, später soll es Java sein.
Die Algorithmen, für die ich obiges Problem benötige, laufen aber im Großen und Ganzen nur einmal - bzw. in größeren zeitlichen Abständen. Zur Laufzeit sind die Berechnungen nicht ganz so umfangreich.

[Roland Ziegler]

@Roland:
Es handelt sich im Moment um folgende Dimensionen - wenn alles funktioniert, können diese sich jedoch noch multiplizieren:
Punkte: < 10.000 in einer Datenbank (im Moment noch MySQL)
Die Anwendung läuft im Moment über PHP5, später soll es Java sein.
Die Algorithmen, für die ich obiges Problem benötige, laufen aber im Großen und Ganzen nur einmal - bzw. in größeren zeitlichen Abständen. Zur Laufzeit sind die Berechnungen nicht ganz so umfangreich.

Dann mach Dir mal keine Sorgen um die Geschwindigkeit durch die Mathematik darin. Das langsamste wird der Datenbankzugriff sein. Und wenn das in Java ausgeführt wird, wirst Du für die reinen Rechenroutinen auch keine großen Einbußen gegenüber nativem Code erleben. 10000 Punkte, das sind beispielsweise drei bis zehn Zeilen einer typischen Bitmap, wie sie in TransDEM unter C++ verarbeitet wird. Beim Georeferenzieren wird da je Bildpunkt das volle Programm durchgezogen, Hin- und Rückprojektion, Datumstransformation. Auf einer Mehrkernmachine eilt hier der mathematische Thread dem Visualisierungsthread immer kräftig voraus (Visualisierung ohne DX). Damit wird klar, was dauert und was schnell geht und wo man ggf. in Optimierung investieren kann oder soll.

[Steffen81]

Das langsamste wird der Datenbankzugriff sein.

Exakt. Und den gilt es mit der hier berechneten Einschränkung zu minimieren.