Wie kann ich mit ArcGIS 10.1 einen geodätischen Punkt mit gleichem Abstand finden, der durch drei Punkte definiert ist?

12

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.

Ich muss O finden

Fett
quelle
Gibt es Einschränkungen bezüglich der relativen Positionen der 3 Punkte? Bild Ostküste, Mittelpunkt ist am weitesten östlich. Ihre Lösung würde nicht funktionieren, da die Lotrechte nicht auf See konvergieren würden. Ich bin sicher, wir können uns andere schlimme Fälle einfallen lassen!
mkennedy
Ich frage mich, ob Sie eine entfernungserhaltende Projektion verwenden und die Berechnung von dort aus ausführen könnten. progonos.com/furuti/MapProj/Normal/CartProp/DistPres/… Ich bin nicht sicher, welchen Algorithmus ich verwenden soll, es muss einen geben ... Vielleicht ist es der Schwerpunkt: en.wikipedia.org/wiki/Barycentric_coordinate_system
Alex Leith
Suchen Sie auf unserer Website nach Lösungen für ein eng verwandtes Problem nach "Trilateration" . Auch gis.stackexchange.com/questions/10332/… ist ein Duplikat, hat aber keine angemessenen Antworten (höchstwahrscheinlich, weil die Frage verwirrt gestellt wurde).
whuber
@mkennedy Es gibt grundsätzlich keine schlechten Fälle, nur numerisch instabile. Diese treten auf, wenn die drei Basispunkte kollinear sind; Die beiden Lösungen (auf einem sphärischen Modell) treten an den beiden Polen der gemeinsamen Geodäsie auf. In einem ellipsoidalen Modell treten sie in der Nähe der erwarteten Pole auf.
whuber
Die Verwendung von Loxodromen wäre hier nicht richtig: Sie sind nicht die senkrechten Halbierungslinien. Auf der Kugel sind diese Linien Teile großer Kreise (Geodäten), auf dem Ellipsoid weichen sie jedoch geringfügig von den geodätischen Linien ab.
whuber

Antworten:

10

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.

Finden Sie für drei Punkte A, B, C (als Lat-Lon-Paare) auf einem Ellipsoid einen Punkt X, für den (1) d (X, A) = d (X, B) = d (X, C) und ( 2) dieser gemeinsame Abstand ist so gering wie möglich.

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

F (Phi, Lambda) = (u (X), v (X))

ist eine Funktion aus einem Teil eines zweidimensionalen Raums, die Werte im zweidimensionalen Raum annimmt, und unser Problem reduziert sich auf

Finde alle möglichen (phi, lambda) für die F (phi, lambda) = (0,0) ist.

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.

sphericalMean[points_] := Module[{sToC, cToS, cMean},
  sToC[{f_, l_}] := {Cos[f] Cos[l], Cos[f] Sin[l], Sin[f]};
  cToS[{x_, y_, z_}] := {ArcTan[x, y], ArcTan[Norm[{x, y}], z]};
  cMean = Mean[sToC /@ (points Degree)];
  If[Norm[Most@cMean] < 10^(-8), Mean[points], cToS[cMean]] / Degree
  ]

(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 im ArcTan Vergleich 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:

tri[a_, b_, c_] := Block[{d = sphericalMean[{a, b, c}], sol, f, q},
   sol = FindRoot[{GeoDistance[{Mod[f, 180, -90], Mod[q, 360, -180]}, a] == 
                   GeoDistance[{Mod[f, 180, -90], Mod[q, 360, -180]}, b] ==
                   GeoDistance[{Mod[f, 180, -90], Mod[q, 360, -180]}, c]}, 
           {{f, d[[1]]}, {q, d[[2]]}}, 
           MaxIterations -> 1000, AccuracyGoal -> Infinity, PrecisionGoal -> 8];
   {Mod[f, 180, -90], Mod[q, 360, -180]} /. sol
   ];

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 Steuerparameter MaxIterationsusw. 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 :

sol = tri @@ (bases = {{-6.28530175, 106.9004975375}, {-6.28955287, 106.89573839}, {-6.28388865789474, 106.908087643421}})

{-6.29692, 106.907}

Die berechneten Abstände zwischen dieser Lösung und den drei Punkten sind

{1450.23206979, 1450.23206979, 1450.23206978}

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

Abbildung 1


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

Figur 2

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 Funktionen FindRootdes 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.

whuber
quelle
3
+1 Sehr gründlich. Ich denke, Sie müssen möglicherweise anfangen, dafür Gebühren zu erheben, @whuber.
Radar
2
@ Radar Danke. Ich hoffe, die Leute werden mein Buch kaufen, wenn es (jemals) irgendwann erscheint :-).
whuber
1
Will do Bill ... Sende einen Entwurf !!!
Ausgezeichnet! Dennoch scheint eine analytische Lösung möglich zu sein. Durch Wiederholung des Problems in den kartesischen 3D-Raum mit 3 Punkten A, B, C und E, wobei E der Mittelpunkt der Erde ist. Als nächstes finden Sie zwei Ebenen Plane1 und Plane2. Ebene1 ist die Ebene, die normal zu EbeneABE ist und durch E, Mittelpunkt (A, B) verläuft. Ebenso wäre Ebene2 die Ebene, die normal zu Ebene A ist und durch E, Mittelpunkt (C, E) verläuft. Die Linie O, die durch Schnitt von Ebene 1 und Ebene 2 gebildet wird, stellt Punkte dar, die gleich weit von den drei Punkten entfernt sind. Der nähere der beiden Punkte zu A (oder B oder C), wo die Linie O die Kugel schneidet, ist Punkt O.
Kirk Kuykendall
Diese analytische Lösung, @Kirk, gilt nur für die Kugel. (Überschneidungen von Ebenen mit dem Ellipsoid sind mit Ausnahme einiger offensichtlicher Ausnahmefälle, bei denen es sich um Meridiane oder den Äquator handelt, keine senkrechten Winkelhalbierenden in der Metrik des Ellipsoids .)
whuber
3

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:

  1. rate einen Dreipunkt O
  2. Verwenden Sie O als Mittelpunkt einer azimutalen äquidistanten Projektion
  3. projizieren Sie A, B, C zu dieser Projektion
  4. finde den Dreipunkt in dieser Projektion, O '
  5. Verwenden Sie O 'als neues Projektionszentrum
  6. Wiederholen, bis O 'und O zusammenfallen

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

function [lat, lon] = tripoint(lat1, lon1, lat2, lon2, lat3, lon3)
% Compute point equidistant from arguments
% Requires:
%   http://www.mathworks.com/matlabcentral/fileexchange/39108
%   http://www.mathworks.com/matlabcentral/fileexchange/39366
  lats = [lat1, lat2, lat3];
  lons = [lon1, lon2, lon3];
  lat0 = lat1;  lon0 = lon1; % feeble guess for tri point
  for i = 1:6
    [x, y] = eqdazim_fwd(lat0, lon0, lats, lons);
    a = [x(1), y(1), 0];
    b = [x(2), y(2), 0];
    c = [x(3), y(3), 0];
    z = [0, 0, 1];
    % Eq. (97) of http://arxiv.org/abs/1102.1215
    o = cross((a*a') * (b - c) + (b*b') * (c - a) + (c*c') * (a - b), z) ...
        / (2 * dot(cross(a - b, b - c), z));
    [lat0, lon0] = eqdazim_inv(lat0, lon0, o(1), o(2))
  end
  % optional check
  s12 = geoddistance(lat0, lon0, lats, lons); ds12 = max(s12) - min(s12)
  lat = lat0; lon = lon0;
end

Ich teste das mit Octave

Oktave: 1> Format lang
Oktave: 2> [lat0, lon0] = Tripoint (41, -74, 36, 140, -41, 175)
lat0 = 15.4151378380375
lon0 = -162.479314381144
lat0 = 15.9969703299812
lon0 = -147.046790722192
lat0 = 16.2232960167545
lon0 = -147.157646039471
lat0 = 16.2233394851560
lon0 = -147.157748279290
lat0 = 16.2233394851809
lon0 = -147.157748279312
lat0 = 16.2233394851809
lon0 = -147.157748279312
ds12 = 3,72529029846191e-09
lat0 = 16.2233394851809
lon0 = -147.157748279312

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.

function [lat2,lon2] = semiconj(lat1, lon1, azi1)
% Find the point where neighboring parallel geodesics emanating from
% close to [lat1, lon1] with azimuth azi1 intersect.

  % First guess is 90 deg on aux sphere
  [lat2, lon2, ~, ~, m12, M12, M21, s12] = ...
      geodreckon(lat1, lon1, 90, azi1, defaultellipsoid(), true);
  M12
  % dM12/ds2 = - (1 - M12*M21/m12)
  for i = 1:3
    s12 = s12 - M12 / ( -(1 - M12*M21)/m12 ); % Newton
    [lat2, lon2, ~, ~, m12, M12, M21] = geodreckon(lat1, lon1, s12, azi1);
    M12
  end
end

Für das Beispiel ergibt dies:

[lat2, lon2] = semiconj (45, 0, 90)
M12 = 0,00262997817649321
M12 = -6.08402492665097e-09
M12 = 4.38017677684144e-17
M12 = 4.38017677684144e-17
lat2 = -0,106909087322479
lon2 = 89.8929107279317
cffk
quelle
+1. Es ist jedoch unklar, ob ein allgemeiner Root-Finder eine schlechtere Leistung erbringt: Die Funktion verhält sich in der Nähe ihres Optimums und Newtons Methode wird beispielsweise auch quadratisch konvergieren. (Die Konvergenz von Mathematica dauert in der Regel ungefähr vier Schritte. Für jeden Schritt sind vier Auswertungen erforderlich, um den Jacobian zu schätzen.) Der wahre Vorteil Ihrer Methode besteht darin, dass sie ohne Verwendung eines Root-Finders problemlos in einem GIS skriptiert werden kann.
Whuber
Genau. Meine Methode ist äquivalent zu Newton und im Gegensatz zu Mathematics Wurzelfindungsmethode ist es nicht erforderlich, den Gradienten anhand von Differenzen abzuschätzen.
cffk
Richtig - aber Sie müssen jedes Mal die Neuprojektion durchführen, was so aussieht, als wäre es ungefähr gleich viel Arbeit. Ich schätze jedoch die Einfachheit und Eleganz Ihres Ansatzes: Es ist sofort klar, dass es funktionieren sollte und schnell konvergieren wird.
Whuber
Ich habe in meiner Antwort Ergebnisse für dieselben Testpunkte veröffentlicht.
Kirk Kuykendall
3

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:

0 longitude: 0 latitude: 90
    Distances: 3134.05443974188 2844.67237777542 3234.33025754997
    Diffs: 289.382061966458 -389.657879774548 -100.27581780809
1 longitude: 106.906152157596 latitude: -6.31307123035178
    Distances: 1450.23208989615 1450.23208089398 1450.23209429293
    Diffs: 9.00216559784894E-06 -1.33989510686661E-05 -4.39678547081712E-06
2 longitude: 106.906583669013 latitude: -6.29691590176649
    Distances: 1450.23206976414 1450.23206976408 1450.23206976433
    Diffs: 6.18456397205591E-11 -2.47382558882236E-10 -1.85536919161677E-10
3 longitude: 106.906583669041 latitude: -6.29691590154641
    Distances: 1450.23206976438 1450.23206976423 1450.23206976459
    Diffs: 1.47565515362658E-10 -3.61751517630182E-10 -2.14186002267525E-10
4 longitude: 106.906583669041 latitude: -6.29691590154641
    Distances: 1450.23206976438 1450.23206976423 1450.23206976459
    Diffs: 1.47565515362658E-10 -3.61751517630182E-10 -2.14186002267525E-10

Hier ist der Quellcode. (Bearbeiten) Das FindCircleCenter wurde geändert, um Schnittpunkte (Mittelpunkt) zu behandeln, die vom Rand der Azimutalprojektion abfallen:

public static void Test()
{
    var t = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
    var srf = Activator.CreateInstance(t) as ISpatialReferenceFactory2;
    var pcs = srf.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_WGS1984N_PoleAziEqui)
        as IProjectedCoordinateSystem2;

    var pntA = MakePoint(106.9004975375, -6.28530175, pcs.GeographicCoordinateSystem);
    var pntB = MakePoint(106.89573839, -6.28955287, pcs.GeographicCoordinateSystem);
    var pntC = MakePoint(106.908087643421, -6.28388865789474, pcs.GeographicCoordinateSystem);

    int maxIter = 5;
    for (int i = 0; i < maxIter; i++)
    {
        var msg = string.Format("{0} longitude: {1} latitude: {2}", i, pcs.get_CentralMeridian(true), pcs.LatitudeOfOrigin);
        Debug.Print(msg);
        var newCenter = FindCircleCenter(ProjectClone(pntA, pcs), ProjectClone(pntB, pcs), ProjectClone(pntC, pcs));
        newCenter.Project(pcs.GeographicCoordinateSystem); // unproject
        var distA = GetGeodesicDistance(newCenter, pntA);
        var distB = GetGeodesicDistance(newCenter, pntB);
        var distC = GetGeodesicDistance(newCenter, pntC);
        Debug.Print("\tDistances: {0} {1} {2}", distA, distB, distC);
        var diffAB = distA - distB;
        var diffBC = distB - distC;
        var diffAC = distA - distC;
        Debug.Print("\tDiffs: {0} {1} {2}", diffAB, diffBC, diffAC);

        pcs.set_CentralMeridian(true, newCenter.X);
        pcs.LatitudeOfOrigin = newCenter.Y;
    }
}
public static IPoint FindCircleCenter(IPoint a, IPoint b, IPoint c)
{
    // from http://blog.csharphelper.com/2011/11/08/draw-a-circle-through-three-points-in-c.aspx
    // Get the perpendicular bisector of (x1, y1) and (x2, y2).
    var x1 = (b.X + a.X) / 2;
    var y1 = (b.Y + a.Y) / 2;
    var dy1 = b.X - a.X;
    var dx1 = -(b.Y - a.Y);

    // Get the perpendicular bisector of (x2, y2) and (x3, y3).
    var x2 = (c.X + b.X) / 2;
    var y2 = (c.Y + b.Y) / 2;
    var dy2 = c.X - b.X;
    var dx2 = -(c.Y - b.Y);

    // See where the lines intersect.
    var cx = (y1 * dx1 * dx2 + x2 * dx1 * dy2 - x1 * dy1 * dx2 - y2 * dx1 * dx2)
        / (dx1 * dy2 - dy1 * dx2);
    var cy = (cx - x1) * dy1 / dx1 + y1;

    // make sure the intersection point falls
    // within the projection.
    var earthRadius = ((IProjectedCoordinateSystem)a.SpatialReference).GeographicCoordinateSystem.Datum.Spheroid.SemiMinorAxis;

    // distance is from center of projection
    var dist = Math.Sqrt((cx * cx) + (cy * cy));
    double factor = 1.0;
    if (dist > earthRadius * Math.PI)
    {
        // apply a factor so we don't fall off the edge
        // of the projection
        factor = earthRadius / dist;
    }
    var outPoint = new PointClass() as IPoint;
    outPoint.PutCoords(cx * factor, cy* factor);
    outPoint.SpatialReference = a.SpatialReference;
    return outPoint;
}

public static double GetGeodesicDistance(IPoint pnt1, IPoint pnt2)
{
    var pc = new PolylineClass() as IPointCollection;
    var gcs = pnt1.SpatialReference as IGeographicCoordinateSystem;
    if (gcs == null)
        throw new Exception("point does not have a gcs");
    ((IGeometry)pc).SpatialReference = gcs;
    pc.AddPoint(pnt1);
    pc.AddPoint(pnt2);
    var t = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
    var srf = Activator.CreateInstance(t) as ISpatialReferenceFactory2;
    var unit = srf.CreateUnit((int)esriSRUnitType.esriSRUnit_Meter) as ILinearUnit;
    var pcGeodetic = pc as IPolycurveGeodetic;
    return pcGeodetic.get_LengthGeodetic(esriGeodeticType.esriGeodeticTypeGeodesic, unit);
}

public static IPoint ProjectClone(IPoint pnt, ISpatialReference sr)
{
    var clone = ((IClone)pnt).Clone() as IPoint;
    clone.Project(sr);
    return clone;
}

public static IPoint MakePoint(double longitude, double latitude, ISpatialReference sr)
{
    var pnt = new PointClass() as IPoint;
    pnt.PutCoords(longitude, latitude);
    pnt.SpatialReference = sr;
    return pnt;
}

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:

0 0
0 longitude: 0 latitude: 0
    MaxDiff: 1859074.90170379 Distances: 13541157.6493561 11682082.7476523 11863320.2116807
1 longitude: 43.5318402621384 latitude: -17.1167429904981
    MaxDiff: 21796.9793742411 Distances: 12584188.7592282 12588146.4851222 12566349.505748
2 longitude: 32.8331167578493 latitude: -16.2707976739314
    MaxDiff: 6.05585224926472 Distances: 12577536.3369782 12577541.3560203 12577542.3928305
3 longitude: 32.8623898057665 latitude: -16.1374156408507
    MaxDiff: 5.58793544769287E-07 Distances: 12577539.6118671 12577539.6118666 12577539.6118669
4 longitude: -147.137582018133 latitude: 16.1374288796667
    MaxDiff: 1.12284109462053 Distances: 7441375.08265703 7441376.12671342 7441376.20549812
5 longitude: -147.157742373074 latitude: 16.2233413614432
    MaxDiff: 7.45058059692383E-09 Distances: 7441375.70752843 7441375.70752842 7441375.70752842
5 longitude: -147.157742373074 latitude: 16.2233413614432 Distance 7441375.70752843
iterations: 5

Hier ist der überarbeitete Code:

class Program
{
    private static LicenseInitializer m_AOLicenseInitializer = new tripoint.LicenseInitializer();

    [STAThread()]
    static void Main(string[] args)
    {
        //ESRI License Initializer generated code.
        m_AOLicenseInitializer.InitializeApplication(new esriLicenseProductCode[] { esriLicenseProductCode.esriLicenseProductCodeStandard },
        new esriLicenseExtensionCode[] { });
        try
        {
            var t = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
            var srf = Activator.CreateInstance(t) as ISpatialReferenceFactory2;
            var pcs = srf.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_World_AzimuthalEquidistant)
                as IProjectedCoordinateSystem2;
            Debug.Print("{0} {1}", pcs.get_CentralMeridian(true), pcs.LatitudeOfOrigin);
            int max = int.MinValue;
            for (int i = 0; i < 1; i++)
            {
                var iterations = Test(pcs);
                max = Math.Max(max, iterations);
                Debug.Print("iterations: {0}", iterations);
            }
            Debug.Print("max number of iterations: {0}", max);
        }
        catch (Exception ex)
        {
            Debug.Print(ex.Message);
            Debug.Print(ex.StackTrace);
        }
        //ESRI License Initializer generated code.
        //Do not make any call to ArcObjects after ShutDownApplication()
        m_AOLicenseInitializer.ShutdownApplication();
    }
    public static int Test(IProjectedCoordinateSystem2 pcs)
    {
        var pntA = MakePoint(-74.0, 41.0, pcs.GeographicCoordinateSystem);
        var pntB = MakePoint(140.0, 36.0, pcs.GeographicCoordinateSystem);
        var pntC = MakePoint(175.0, -41.0, pcs.GeographicCoordinateSystem);


        //var r = new Random();
        //var pntA = MakeRandomPoint(r, pcs.GeographicCoordinateSystem);
        //var pntB = MakeRandomPoint(r, pcs.GeographicCoordinateSystem);
        //var pntC = MakeRandomPoint(r, pcs.GeographicCoordinateSystem);

        int maxIterations = 100;
        for (int i = 0; i < maxIterations; i++)
        {
            var msg = string.Format("{0} longitude: {1} latitude: {2}", i, pcs.get_CentralMeridian(true), pcs.LatitudeOfOrigin);
            Debug.Print(msg);
            var newCenter = FindCircleCenter(ProjectClone(pntA, pcs), ProjectClone(pntB, pcs), ProjectClone(pntC, pcs));
            var c = ((IClone)newCenter).Clone() as IPoint;
            newCenter.Project(pcs.GeographicCoordinateSystem); // unproject
            //newCenter = MakePoint(-147.1577482, 16.2233394, pcs.GeographicCoordinateSystem);
            var distA = GetGeodesicDistance(newCenter, pntA);
            var distB = GetGeodesicDistance(newCenter, pntB);
            var distC = GetGeodesicDistance(newCenter, pntC);
            var diffAB = Math.Abs(distA - distB);
            var diffBC = Math.Abs(distB - distC);
            var diffAC = Math.Abs(distA - distC);
            var maxDiff = GetMax(new double[] {diffAB,diffAC,diffBC});
            Debug.Print("\tMaxDiff: {0} Distances: {1} {2} {3}",maxDiff, distA, distB, distC);
            if (maxDiff < 0.000001)
            {
                var earthRadius = pcs.GeographicCoordinateSystem.Datum.Spheroid.SemiMinorAxis;
                if (distA > earthRadius * Math.PI / 2.0)
                {
                    newCenter = AntiPode(newCenter);
                }
                else
                {
                    Debug.Print("{0} longitude: {1} latitude: {2} Distance {3}", i, pcs.get_CentralMeridian(true), pcs.LatitudeOfOrigin, distA);
                    return i;
                }
            }
            //Debug.Print("\tDiffs: {0} {1} {2}", diffAB, diffBC, diffAC);

            pcs.set_CentralMeridian(true, newCenter.X);
            pcs.LatitudeOfOrigin = newCenter.Y;
        }
        return maxIterations;
    }

    public static IPoint FindCircleCenter(IPoint a, IPoint b, IPoint c)
    {
        // from http://blog.csharphelper.com/2011/11/08/draw-a-circle-through-three-points-in-c.aspx
        // Get the perpendicular bisector of (x1, y1) and (x2, y2).
        var x1 = (b.X + a.X) / 2;
        var y1 = (b.Y + a.Y) / 2;
        var dy1 = b.X - a.X;
        var dx1 = -(b.Y - a.Y);

        // Get the perpendicular bisector of (x2, y2) and (x3, y3).
        var x2 = (c.X + b.X) / 2;
        var y2 = (c.Y + b.Y) / 2;
        var dy2 = c.X - b.X;
        var dx2 = -(c.Y - b.Y);

        // See where the lines intersect.
        var cx = (y1 * dx1 * dx2 + x2 * dx1 * dy2 - x1 * dy1 * dx2 - y2 * dx1 * dx2)
            / (dx1 * dy2 - dy1 * dx2);
        var cy = (cx - x1) * dy1 / dx1 + y1;

        // make sure the intersection point falls
        // within the projection.
        var earthRadius = ((IProjectedCoordinateSystem)a.SpatialReference).GeographicCoordinateSystem.Datum.Spheroid.SemiMinorAxis;

        // distance is from center of projection
        var dist = Math.Sqrt((cx * cx) + (cy * cy));
        double factor = 1.0;
        if (dist > earthRadius * Math.PI)
        {
            // apply a factor so we don't fall off the edge
            // of the projection
            factor = earthRadius / dist;
        }
        var outPoint = new PointClass() as IPoint;
        outPoint.PutCoords(cx * factor, cy* factor);
        outPoint.SpatialReference = a.SpatialReference;
        return outPoint;
    }

    public static IPoint AntiPode(IPoint pnt)
    {
        if (!(pnt.SpatialReference is IGeographicCoordinateSystem))
            throw new Exception("antipode of non-gcs projection not supported");
        var outPnt = new PointClass() as IPoint;
        outPnt.SpatialReference = pnt.SpatialReference;
        if (pnt.X < 0.0)
            outPnt.X = 180.0 + pnt.X;
        else
            outPnt.X = pnt.X - 180.0;
        outPnt.Y = -pnt.Y;
        return outPnt;
    }

    public static IPoint MakeRandomPoint(Random r, IGeographicCoordinateSystem gcs)
    {
        var latitude = (r.NextDouble() - 0.5) * 180.0;
        var longitude = (r.NextDouble() - 0.5) * 360.0;
        //Debug.Print("{0} {1}", latitude, longitude);
        return MakePoint(longitude, latitude, gcs);
    }
    public static double GetMax(double[] dbls)
    {
        var max = double.MinValue;
        foreach (var d in dbls)
        {
            if (d > max)
                max = d;
        }
        return max;
    }
    public static IPoint MakePoint(IPoint[] pnts)
    {
        double sumx = 0.0;
        double sumy = 0.0;
        foreach (var pnt in pnts)
        {
            sumx += pnt.X;
            sumy += pnt.Y;
        }
        return MakePoint(sumx / pnts.Length, sumy / pnts.Length, pnts[0].SpatialReference);
    }
    public static double GetGeodesicDistance(IPoint pnt1, IPoint pnt2)
    {
        var pc = new PolylineClass() as IPointCollection;
        var gcs = pnt1.SpatialReference as IGeographicCoordinateSystem;
        if (gcs == null)
            throw new Exception("point does not have a gcs");
        ((IGeometry)pc).SpatialReference = gcs;
        pc.AddPoint(pnt1);
        pc.AddPoint(pnt2);
        var t = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
        var srf = Activator.CreateInstance(t) as ISpatialReferenceFactory2;
        var unit = srf.CreateUnit((int)esriSRUnitType.esriSRUnit_Meter) as ILinearUnit;
        var pcGeodetic = pc as IPolycurveGeodetic;
        return pcGeodetic.get_LengthGeodetic(esriGeodeticType.esriGeodeticTypeGeodesic, unit);
    }

    public static IPoint ProjectClone(IPoint pnt, ISpatialReference sr)
    {
        var clone = ((IClone)pnt).Clone() as IPoint;
        clone.Project(sr);
        return clone;
    }

    public static IPoint MakePoint(double longitude, double latitude, ISpatialReference sr)
    {
        var pnt = new PointClass() as IPoint;
        pnt.PutCoords(longitude, latitude);
        pnt.SpatialReference = sr;
        return pnt;
    }
}

Bearbeiten

Hier sind die Ergebnisse, die ich mit esriSRProjCS_WGS1984N_PoleAziEqui erhalte

0 90
0 longitude: 0 latitude: 90
    MaxDiff: 1275775.91880553 Distances: 8003451.67666723 7797996.2370572 6727675.7578617
1 longitude: -148.003774863594 latitude: 9.20238223616225
    MaxDiff: 14487.6784785809 Distances: 7439006.46128994 7432752.45732905 7447240.13580763
2 longitude: -147.197808459106 latitude: 16.3073233548167
    MaxDiff: 2.32572609744966 Distances: 7441374.94409209 7441377.26981819 7441375.90768183
3 longitude: -147.157734641831 latitude: 16.2233338760411
    MaxDiff: 7.72997736930847E-08 Distances: 7441375.70752842 7441375.70752848 7441375.7075284
3 longitude: -147.157734641831 latitude: 16.2233338760411 Distance 7441375.70752842
Kirk Kuykendall
quelle
Das ist beeindruckend schnelle Konvergenz! (+1)
whuber
Sie sollten eine ehrliche azimutale äquidistante Projektion verwenden, die auf newCenter ausgerichtet ist. Stattdessen verwenden Sie die Projektion, die auf dem N-Pol zentriert ist, und verschieben den Ursprung auf newCenter. Es kann daher ein Zufall sein, dass Sie in diesem Fall eine anständige Lösung erhalten (vielleicht, weil die Punkte nahe beieinander liegen?). Es wäre gut, es mit 3 Punkten zu versuchen, die Tausende von Kilometern voneinander entfernt sind. Eine Implementierung der azimutalen äquidistanten Projektion finden Sie in mathworks.com/matlabcentral/fileexchange/…
cffk am
@cffk Die einzige andere Möglichkeit, eine azimutale äquidistante Projektion zu erstellen, die auf einem bestimmten Punkt zentriert ist, besteht darin, dieselbe Methode zu verwenden, jedoch mit esriSRProjCS_World_AzimuthalEquidistant anstelle von esriSRProjCS_WGS1984N_PoleAziEqui (oder esriSRProjCS_S_S_S_S_S_S_S_ Der einzige Unterschied besteht darin, dass es auf 0,0 anstatt auf 0,90 (oder 0, -90) zentriert ist. Können Sie mich anleiten, einen Test mit der Mathematik durchzuführen, um zu sehen, ob dies zu unterschiedlichen Ergebnissen aus einer Projektion von "Ehrlichkeit zu Güte" führt?
Kirk Kuykendall
@ KirkKuykendall: siehe den Anhang zu meiner ersten Antwort.
cffk
1
@KirkKuykendall Vielleicht hat ESRI eine "ehrliche" Projektion implementiert? Die Schlüsseleigenschaft, die erforderlich ist, damit dieser Algorithmus funktioniert, ist, dass die vom "Mittelpunkt" gemessenen Abstände wahr sind. Und diese Eigenschaft ist leicht zu überprüfen. (Die azimutale Eigenschaft in Bezug auf den Mittelpunkt ist zweitrangig und kann nur die Konvergenzrate beeinflussen.)
cffk