Zum Beispiel habe ich Koordinaten für drei Basispunkte an einer Küste und ich muss die Koordinaten des Punktes vor der Küste finden, der von allen drei gleich weit entfernt ist. Es ist eine einfache Übung in Geometrie, aber alle Messungen müssen die Geodäsie berücksichtigen.
Wenn ich mich diesem Punkt auf euklidische Weise näherte, könnte ich die geodätischen Pfade messen, die die Basispunkte verbinden, die Mittelpunkte der Seiten des resultierenden Dreiecks finden und senkrechte Orthodrome zu jedem dieser Pfade erstellen. Die drei Loxodrome würden vermutlich am äquidistanten Punkt zusammenlaufen. Wenn dies die richtige Methode ist, muss es in Arc einen einfacheren Weg geben.
Antworten:
Diese Antwort ist in mehrere Abschnitte unterteilt:
Analyse und Reduzierung des Problems , um zu zeigen, wie man den gewünschten Punkt mit "eingemachten" Routinen findet.
Abbildung: Ein funktionierender Prototyp mit Arbeitscode.
Beispiel , zeigt Beispiele für die Lösungen.
Fallstricke , mögliche Probleme besprechen und wie man mit ihnen umgeht .
ArcGIS-Implementierung , Kommentare zum Erstellen eines benutzerdefinierten ArcGIS-Tools und zum Abrufen der erforderlichen Routinen.
Analyse und Reduzierung des Problems
Beginnen wir mit der Beobachtung , dass in der (vollkommen rund) sphärischen Modell gibt es immer eine Lösung sein --in der Tat, genau zwei Lösungen. Bei gegebenen Basispunkten A, B und C bestimmt jedes Paar seine "senkrechte Winkelhalbierende", die die Menge von Punkten ist, die von den beiden gegebenen Punkten gleich weit entfernt sind. Diese Halbierende ist eine Geodätische (Großkreis). Die sphärische Geometrie ist elliptisch : Zwei beliebige Geodäten kreuzen sich (in zwei eindeutigen Punkten). Somit sind die Schnittpunkte der AB-Halbierungslinie und der BC-Halbierungslinie per Definition äquidistant von A, B und C, wodurch das Problem gelöst wird. (Siehe die erste Abbildung unten.)
Auf einem Ellipsoid sieht es komplizierter aus, aber da es sich um eine kleine Störung der Kugel handelt, können wir ein ähnliches Verhalten erwarten. (Die Analyse würde uns zu weit führen.) Die komplizierten Formeln, die (innerhalb eines GIS) verwendet werden, um genaue Entfernungen auf einem Ellipsoid zu berechnen, sind jedoch keine konzeptionelle Komplikation: Das Problem ist im Grunde dasselbe. Um zu sehen, wie einfach das Problem wirklich ist, geben wir es etwas abstrakter an. In dieser Aussage bezieht sich "d (U, V)" auf den wahren, vollständig genauen Abstand zwischen den Punkten U und V.
Diese drei Entfernungen hängen alle vom unbekannten X ab . So werden die Unterschiede in den Abständen u (X) = d (X, A) - D (X, B) und V (x) = d (X, B) - D (X, C) sind Reellwertige Funktionen von X. Nochmals, etwas abstrakt, können wir diese Unterschiede zu einem geordneten Paar zusammenfassen. Wir werden auch (lat, lon) als Koordinaten für X verwenden, so dass wir es auch als geordnetes Paar betrachten können, sagen wir X = (phi, lambda). In diesem Setup ist die Funktion
ist eine Funktion aus einem Teil eines zweidimensionalen Raums, die Werte im zweidimensionalen Raum annimmt, und unser Problem reduziert sich auf
Hier zahlt sich die Abstraktion aus: Es gibt eine Menge großartiger Software, um dieses (rein numerische, mehrdimensionale Root-Finding-) Problem zu lösen. Die Art und Weise, wie es funktioniert, ist, dass Sie eine Routine schreiben, um F zu berechnen , und diese zusammen mit allen Informationen über Einschränkungen der Eingabe an die Software weitergeben ( phi muss zwischen -90 und 90 Grad liegen und Lambda muss zwischen -180 und 180 liegen Grad). Es dreht für den Bruchteil einer Sekunde ab und gibt (normalerweise) nur einen Wert von ( phi , lambda ) zurück, wenn es einen finden kann.
Es gibt Details zu handhaben, weil das eine Kunst ist: Es gibt verschiedene Lösungsmethoden zur Auswahl, abhängig davon, wie sich F "verhält"; es hilft, die Software zu "steuern", indem es einen vernünftigen Ausgangspunkt für ihre Suche gibt (dies ist eine Möglichkeit, die nächstgelegene Lösung zu erhalten, anstatt irgendeine andere); In der Regel müssen Sie angeben, wie genau die Lösung sein soll (damit Sie wissen, wann Sie die Suche beenden müssen). (Weitere Informationen darüber, was GIS-Analysten über solche Details wissen müssen, die häufig bei GIS-Problemen auftreten, finden Sie unter Empfehlen von Themen, die in einen Kurs für Informatik für Geoinformatik aufgenommen werden sollen, und im Abschnitt "Verschiedenes" am Ende. )
Abbildung: ein funktionierender Prototyp
Die Analyse zeigt, dass wir zwei Dinge programmieren müssen: eine grobe Anfangsschätzung der Lösung und die Berechnung von F selbst.
Die anfängliche Schätzung kann durch einen "sphärischen Durchschnitt" der drei Basispunkte erfolgen. Dies wird erreicht, indem sie in geozentrischen kartesischen (x, y, z) Koordinaten dargestellt werden, diese Koordinaten gemittelt werden und dieser Durchschnitt zurück auf die Kugel projiziert und in Breiten- und Längengraden erneut ausgedrückt wird. Die Größe der Kugel spielt keine Rolle und die Berechnungen sind daher unkompliziert: Da dies nur ein Ausgangspunkt ist, benötigen wir keine Ellipsoidberechnungen.
Für diesen funktionierenden Prototyp habe ich Mathematica 8 verwendet.
(Die endgültige
IfBedingung prüft, ob der Durchschnitt möglicherweise nicht eindeutig eine Länge angibt. In diesem Fall wird auf ein gerades arithmetisches Mittel der Breiten- und Längengrade seiner Eingabe zurückgegriffen - möglicherweise keine gute Wahl, aber zumindest eine gültige. Wenn Sie diesen Code als Implementierungsanleitung verwenden, beachten Sie, dass die Argumente von Mathematica imArcTanVergleich zu den meisten anderen Implementierungen umgekehrt sind: Das erste Argument ist die x-Koordinate, das zweite ist die y-Koordinate und es wird der vom Vektor gemachte Winkel zurückgegeben ( x, y).)Da Mathematica - wie ArcGIS und fast alle anderen GIS - Code zur Berechnung genauer Abstände auf dem Ellipsoid enthält, gibt es im zweiten Teil fast nichts zu schreiben. Wir nennen nur die Wurzelfindungsroutine:
Der bemerkenswerteste Aspekt dieser Implementierung ist, wie sie die Notwendigkeit vermeidet, den Breitengrad (
f) und den Längengrad (q) einzuschränken, indem sie immer modulo 180 bzw. 360 Grad berechnen. Dies vermeidet, dass das Problem eingeschränkt werden muss (was häufig zu Komplikationen führt). Die SteuerparameterMaxIterationsusw. wurden so angepasst, dass dieser Code die größtmögliche Genauigkeit bietet.Um es in Aktion zu sehen, wenden wir es auf die drei Basispunkte an, die in einer verwandten Frage angegeben sind :
Die berechneten Abstände zwischen dieser Lösung und den drei Punkten sind
(Dies sind Meter). Sie stimmen durch die elfte signifikante Ziffer überein (was eigentlich zu genau ist, da Abstände selten besser als ein Millimeter oder so genau sind). Hier ist ein Bild dieser drei Punkte (schwarz), ihrer drei gegenseitigen Winkelhalbierenden und der Lösung (rot):
Beispiel
Um diese Implementierung zu testen und ein besseres Verständnis für das Verhalten des Problems zu erhalten, finden Sie hier eine Konturdarstellung der mittleren quadratischen Abweichung in Abständen für drei weit voneinander entfernte Basispunkte. (Die RMS-Diskrepanz wird erhalten, indem alle drei Differenzen d (X, A) -d (X, B), d (X, B) -d (X, C) und d (X, C) -d (X) berechnet werden , A), deren Quadrate gemittelt werden und die Quadratwurzel gezogen wird. Wenn X das Problem löst und ansonsten zunimmt, wenn X von einer Lösung abweicht, entspricht dies Null und misst, wie nahe wir einer Lösung an einem beliebigen Ort sind. )
Die Basispunkte (60, -120), (10, -40) und (45,10) sind in dieser Plate-Carree-Projektion rot dargestellt. Die Lösung (49.2644488, -49.9052992), für deren Berechnung 0,03 Sekunden benötigt wurden, ist gelb. Die RMS-Diskrepanz beträgt weniger als drei Nanometer , obwohl alle relevanten Entfernungen Tausende von Kilometern betragen. Die dunklen Bereiche zeigen kleine Werte des Effektivwerts und die hellen Bereiche zeigen hohe Werte.
Diese Karte zeigt deutlich, dass eine andere Lösung in der Nähe von (-49.2018206, 130.0297177) liegt (berechnet auf einen Effektivwert von zwei Nanometern, indem der anfängliche Suchwert diametral gegenüber der ersten Lösung eingestellt wird).
Tücken
Numerische Instabilität
Wenn die Basispunkte nahezu kollinear und nahe beieinander liegen, sind alle Lösungen fast eine halbe Welt entfernt und äußerst schwer genau zu bestimmen. Der Grund dafür ist, dass kleine Änderungen an einem Ort auf der ganzen Welt - in Richtung der Basispunkte oder von diesen weg - nur unglaublich kleine Änderungen der Entfernungsunterschiede hervorrufen. Es ist einfach nicht genug Genauigkeit und Präzision in die übliche Berechnung von geodätischen Entfernungen eingebaut, um das Ergebnis zu bestimmen.
Wenn Sie beispielsweise mit den Basispunkten (45.001, 0), (45, 0) und (44.999,0) beginnen, die entlang des Nullmeridians nur 111 Meter voneinander entfernt sind, erhalten Sie
tridie Lösung (11.8213, 77.745) ). Die Abstände zu den Basispunkten betragen 8.127.964.998 77; 8.127.964.998 41; bzw. 8.127.964.998 65 Meter. Sie stimmen auf den Millimeter genau überein! Ich bin mir nicht sicher, wie genau dieses Ergebnis sein mag, wäre aber nicht im Geringsten überrascht, wenn andere Implementierungen Standorte zurückgeben würden, die weit von diesem entfernt sind und eine fast so gute Gleichheit der drei Entfernungen aufweisen.Rechenzeit
Diese Berechnungen sind nicht schnell, da sie eine beträchtliche Suche unter Verwendung komplizierter Entfernungsberechnungen erfordern, und erfordern gewöhnlich einen merklichen Sekundenbruchteil. Echtzeitanwendungen müssen dies berücksichtigen.
ArcGIS-Implementierung
Python ist die bevorzugte Skriptumgebung für ArcGIS (ab Version 9). Das Paket scipy.optimize enthält einen multivariaten Rootfinder,
rootder die FunktionenFindRootdes Mathematica- Codes ausführen soll. Natürlich bietet ArcGIS selbst genaue Ellipsoid-Abstandsberechnungen. Der Rest sind dann alle Implementierungsdetails: Entscheiden Sie, wie die Basispunktkoordinaten erhalten werden (von einer Ebene, die der Benutzer eingibt, aus einer Textdatei, von der Maus) und wie die Ausgabe dargestellt wird (als Koordinaten) Schreiben Sie diese Schnittstelle, portieren Sie den hier gezeigten Mathematica- Code (unkompliziert), und schon sind Sie fertig.quelle
Wie Sie bemerken, tritt dieses Problem bei der Bestimmung der Seegrenzen auf. Es wird oft als "Tri-Point" -Problem bezeichnet, und Sie können dies googeln und mehrere Artikel finden, die sich damit befassen. Eines dieser Papiere ist von mir (!) Und ich biete eine genaue und schnell konvergente Lösung. Siehe Abschnitt 14 von http://arxiv.org/abs/1102.1215
Die Methode besteht aus folgenden Schritten:
Die notwendige Formel für die Dreipunktlösung in der Projektion ist in der Arbeit angegeben. Solange Sie eine genaue azimutale äquidistante Projektion verwenden, ist die Antwort genau. Konvergenz ist quadratisch, was bedeutet, dass nur wenige Iterationen erforderlich sind. Dies wird mit ziemlicher Sicherheit die von @whuber vorgeschlagenen allgemeinen Methoden zum Auffinden von Wurzeln übertreffen.
Ich kann Ihnen nicht direkt mit ArcGIS helfen. Sie können mein Python-Paket für geodätische Berechnungen unter https://pypi.python.org/pypi/geographiclib herunterladen und die darauf basierende Projektion ganz einfach codieren.
Bearbeiten
Das Problem, den Dreipunkt in @ whubers entartetem Fall (45 + eps, 0) (45,0) (45-eps, 0) zu finden, wurde von Cayley in Auf den geodätischen Linien auf einem abgeflachten Sphäroid , Phil, betrachtet. Mag. (1870), http://books.google.com/books?id=4XGIOoCMYYAC&pg=PA15
In diesem Fall erhält man den Dreipunkt, indem man einer Geodät von (45,0) mit dem Azimut 90 folgt und den Punkt findet, an dem die geodätische Skala verschwindet. Für das WGS84-Ellipsoid ist dieser Punkt (-0,10690908732248, 89,89291072793167). Die Entfernung von diesem Punkt zu jedem der Punkte (45,001,0), (45,0) und (44,999) beträgt 10010287,665788943 m (innerhalb eines Nanometers oder so). Dies sind ungefähr 1882 km mehr als Whubers Schätzung (was nur zeigt, wie instabil dieser Fall ist). Für eine kugelförmige Erde wäre der Dreipunkt natürlich (0,90) oder (0, -90).
ADDENDUM: Hier ist eine Implementierung der azimutalen äquidistanten Methode mit Matlab
Ich teste das mit Octave
als Tri-Point für New York, Tokio und Wellington.
Diese Methode ist für benachbarte kolineare Punkte ungenau, z. B. [45.001,0], [45,0], [44.999,0]. In diesem Fall sollten Sie für eine Geodät aus [45,0] bei Azimut 90 nach M 12 = 0 auflösen.
Für das Beispiel ergibt dies:
quelle
Ich war gespannt, wie schnell sich @ cffks Ansatz einer Lösung annähert, und schrieb einen Test mit ArcObjects, der diese Ausgabe ergab. Entfernungen sind in Metern:
Hier ist der Quellcode. (Bearbeiten) Das FindCircleCenter wurde geändert, um Schnittpunkte (Mittelpunkt) zu behandeln, die vom Rand der Azimutalprojektion abfallen:
Es gibt auch einen alternativen Ansatz in der Juni 2013 Ausgabe des MSDN Magazins, Amoeba Method Optimization using C # .
Bearbeiten
In einigen Fällen wurde der zuvor veröffentlichte Code zum Antipode konvergiert. Ich habe den Code so geändert, dass er diese Ausgabe für @ cffks Testpunkte erzeugt.
Hier ist die Ausgabe, die es jetzt erzeugt:
Hier ist der überarbeitete Code:
Bearbeiten
Hier sind die Ergebnisse, die ich mit esriSRProjCS_WGS1984N_PoleAziEqui erhalte
quelle