Jak mogę użyć ArcGIS 10.1, aby znaleźć punkt geodezyjny w równej odległości zdefiniowany przez trzy punkty?

Na przykład mam współrzędne trzech punktów bazowych na linii brzegowej i muszę znaleźć współrzędne punktu u wybrzeży, który jest w równej odległości od wszystkich trzech. Jest to proste ćwiczenie z geometrii, ale wszystkie pomiary muszą uwzględniać geodezję.

Gdybym podchodził do tego w sposób euklidesowy, mógłbym zmierzyć ścieżki geodezyjne łączące punkty bazowe, znaleźć punkty środkowe boków powstałego trójkąta i stworzyć prostopadłe ortodromy do każdej z tych ścieżek. Te trzy loksodromy prawdopodobnie zbiegną się w równo odległym punkcie. Jeśli jest to poprawna metoda, musi istnieć łatwiejszy sposób, aby to zrobić w Arc.

Ta odpowiedź jest podzielona na wiele sekcji:

• Analiza i redukcja problemu , pokazująca, jak znaleźć pożądany punkt za pomocą procedur „w puszkach”.

• Ilustracja: działający prototyp , podający działający kod.

• Przykład pokazujący przykłady rozwiązań.

• Pułapki , omawianie potencjalnych problemów i sposoby radzenia sobie z nimi.

• Implementacja ArcGIS , komentarze na temat tworzenia niestandardowego narzędzia ArcGIS i gdzie uzyskać potrzebne procedury.

Analiza i redukcja problemu

Zacznijmy od obserwacji, że w (idealnie okrągły) Model sferycznej będzie zawsze być rozwiązaniem --W rzeczywistości, dokładnie dwa rozwiązania. Biorąc pod uwagę punkty bazowe A, B i C, każda para określa swój „prostopadły dwusieczny”, który jest zbiorem punktów w równej odległości od dwóch podanych punktów. Ten dwusieczny jest geodezyjny (wielki okrąg). Geometria sferyczna jest eliptyczna : przecinają się dowolne dwie geodezy (w dwóch unikalnych punktach). Zatem punkty przecięcia dwusiecznego AB i dwusiecznego BC są - z definicji - w równej odległości od A, B i C, rozwiązując w ten sposób problem. (Zobacz pierwszy rysunek poniżej.)

Elipsoida wygląda na bardziej skomplikowaną, ale ponieważ jest to niewielkie zaburzenie sfery, możemy spodziewać się podobnego zachowania. (Analiza tego zabrałaby nas za daleko). Skomplikowane formuły stosowane (wewnętrznie w GIS) do obliczania dokładnych odległości na elipsoidzie nie są jednak komplikacją pojęciową: problem jest w zasadzie taki sam. Aby zobaczyć, jak prosty jest naprawdę problem, określmy go nieco abstrakcyjnie. W tym stwierdzeniu „d (U, V)” odnosi się do prawdziwej, w pełni dokładnej odległości między punktami U i V.

Biorąc pod uwagę trzy punkty A, B, C (jako pary lat-lon) na elipsoidzie, znajdź punkt X, dla którego (1) d (X, A) = d (X, B) = d (X, C) i ( 2) ta wspólna odległość jest jak najmniejsza.

Te trzy odstępy wszystko zależy od nieznanego X . Zatem różnice w odległościach u (X) = d (X, A) - d (X, B) i v (X) = d (X, B) - d (X, C) są funkcjami X w wartościach rzeczywistych. Ponownie, nieco abstrakcyjnie, możemy połączyć te różnice w uporządkowaną parę. Użyjemy również (lat, lon) jako współrzędnych dla X, co pozwoli nam również rozważyć to jako parę uporządkowaną, powiedzmy X = (phi, lambda). W tej konfiguracji funkcja

F (phi, lambda) = (u (X), v (X))

jest funkcją z części dwuwymiarowej przestrzeni przyjmującej wartości w przestrzeni dwuwymiarowej, a nasz problem zmniejsza się do

Znajdź wszystkie możliwe (phi, lambda), dla których F (phi, lambda) = (0,0).

Oto, gdzie abstrakcja się opłaca: istnieje mnóstwo świetnego oprogramowania do rozwiązania tego (czysto numerycznego wielowymiarowego problemu znajdowania korzeni). Działa to tak, że piszesz procedurę do obliczenia F , a następnie przekazujesz ją do oprogramowania wraz z wszelkimi informacjami o ograniczeniach na jego wejściu ( phi musi leżeć między -90 a 90 stopni, a lambda musi leżeć między -180 a 180 stopni). Rozwija się na ułamek sekundy i zwraca (zazwyczaj) tylko jedną wartość ( phi , lambda ), jeśli ją znajdzie.

Są szczegóły do ​​rozwiązania, ponieważ jest w tym sztuka: istnieją różne metody rozwiązania do wyboru, w zależności od tego, jak F „zachowuje się”; pomaga „sterować” oprogramowaniem, zapewniając mu odpowiedni punkt wyjścia do wyszukiwania (jest to jeden ze sposobów uzyskania najbliższego rozwiązania, a nie inne); i zazwyczaj musisz określić, jak dokładne ma być rozwiązanie (aby wiedział, kiedy zatrzymać wyszukiwanie). (Aby dowiedzieć się więcej na temat tego, co analitycy GIS powinni wiedzieć o takich szczegółach, które często pojawiają się w problemach z GIS, odwiedź Zalecane tematy, które należy uwzględnić w kursie informatyki dla technologii geoprzestrzennych i zajrzyj do sekcji „Różne” pod koniec. )

Ilustracja: działający prototyp

Analiza pokazuje, że musimy zaprogramować dwie rzeczy: przybliżone wstępne oszacowanie rozwiązania i obliczenie samego F.

Wstępnej oceny można dokonać na podstawie „średniej sferycznej” trzech punktów bazowych. Uzyskuje się to poprzez przedstawienie ich w geocentrycznych współrzędnych kartezjańskich (x, y, z), uśrednienie tych współrzędnych i rzutowanie tej średniej z powrotem na sferę i ponowne wyrażenie jej w szerokości i długości geograficznej. Rozmiar kuli jest nieistotny, dlatego obliczenia są proste: ponieważ jest to tylko punkt początkowy, nie potrzebujemy obliczeń elipsoidalnych.

Do tego działającego prototypu użyłem Mathematica 8.

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
  ]

( IfWarunek końcowy sprawdza, czy średnia może nie wskazywać wyraźnie długości geograficznej; jeśli tak, to wraca do prostej arytmetycznej szerokości i długości geograficznej jej danych wejściowych - może nie jest to doskonały wybór, ale przynajmniej poprawny. Dla tych, którzy używają tego kodu do wskazówek dotyczących implementacji, zauważ, że argumenty Mathematiki ArcTan są odwrócone w porównaniu z większością innych implementacji: jej pierwszy argument to współrzędna x, drugi to współrzędna y i zwraca kąt wykonany przez wektor ( x, y).)

Jeśli chodzi o drugą część, ponieważ Mathematica - podobnie jak ArcGIS i prawie wszystkie inne GIS - zawiera kod do obliczania dokładnych odległości na elipsoidzie, prawie nie ma co pisać. Po prostu nazywamy procedurę znajdowania roota:

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
   ];

Najbardziej godnym uwagi aspektem tej implementacji jest to, że unika ona potrzeby ograniczania szerokości fi szerokości geograficznej ( q), zawsze obliczając je odpowiednio modulo 180 i 360 stopni. Pozwala to uniknąć konieczności ograniczania problemu (co często powoduje komplikacje). Parametry kontrolne MaxIterationsitp. Są modyfikowane, aby kod ten zapewniał najwyższą możliwą dokładność.

Aby zobaczyć to w akcji, zastosujmy ją do trzech punktów bazowych podanych w powiązanym pytaniu :

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

{-6.29692, 106.907}

Obliczone odległości między tym rozwiązaniem a trzema punktami wynoszą

{1450.23206979, 1450.23206979, 1450.23206978}

(to są metry). Zgadzają się co do jedenastej cyfry znaczącej (która jest zbyt precyzyjna, ponieważ odległości rzadko są dokładne z dokładnością lepszą niż milimetr). Oto zdjęcie tych trzech punktów (czarny), ich trzech wzajemnych dwusiecznych i rozwiązania (czerwony):

Przykład

Aby przetestować tę implementację i lepiej zrozumieć, jak zachowuje się problem, oto wykres konturowy średniej kwadratowej rozbieżności w odległościach dla trzech szeroko rozmieszczonych punktów bazowych. (Rozbieżność RMS uzyskuje się przez obliczenie wszystkich trzech różnic d (X, A) -d (X, B), d (X, B) -d (X, C) i d (X, C) -d (X , A), uśredniając ich kwadraty i przyjmując pierwiastek kwadratowy. Jest równy zero, gdy X rozwiązuje problem, a w przeciwnym razie wzrasta, gdy X odchodzi od rozwiązania, a zatem mierzy, jak „jesteśmy blisko” rozwiązania w dowolnym miejscu. )

Punkty bazowe (60, -120), (10, -40) i (45,10) są pokazane na czerwono w tym rzucie płyty Carree; rozwiązanie (49,2644488, -49,9052992) - którego obliczenie wymagało 0,03 sekundy - jest w kolorze żółtym. Jego rozbieżność RMS wynosi mniej niż trzy nanometry , mimo że wszystkie istotne odległości wynoszą tysiące kilometrów. Ciemne obszary pokazują małe wartości RMS, a jasne obszary pokazują wysokie wartości.

Ta mapa wyraźnie pokazuje, że inne rozwiązanie leży w pobliżu (-49.2018206, 130.0297177) (obliczone na RMS dwóch nanometrów poprzez ustawienie początkowej wartości wyszukiwania diametralnie przeciwnej do pierwszego rozwiązania).

Pułapki

Niestabilność numeryczna

Kiedy punkty bazowe są prawie współliniowe i zbliżają się do siebie, wszystkie rozwiązania będą blisko pół świata i niezwykle trudno będzie je dokładnie określić. Powodem jest to, że niewielkie zmiany w lokalizacji na całym świecie - przenoszenie jej w kierunku punktów bazowych lub od nich - powodują tylko niewiarygodnie małe zmiany w różnicach odległości. Po prostu nie ma wystarczającej dokładności i precyzji wbudowanych w zwykłe obliczenia odległości geodezyjnych, aby dokładnie określić wynik.

Na przykład, zaczynając od punktów bazowych w (45,001, 0), (45, 0) i (44,999,0), które są oddzielone wzdłuż Prime Meridian tylko 111 metrów między każdą parą, triuzyskuje rozwiązanie (11.8213, 77.745 ). Odległości od punktów bazowych wynoszą 8 127 964 998 77; 8 127 964 998 41; i 8 127 964 998, odpowiednio 65 metrów. Zgadzają się z dokładnością do milimetra! Nie jestem pewien, jak dokładny może być ten wynik, ale nie byłbym w najmniejszym stopniu zaskoczony, gdyby inne implementacje zwróciły lokalizacje daleko od tego, pokazując prawie równie dobrą równość trzech odległości.

Czas obliczeń

Obliczenia te, ponieważ wymagają znacznego wyszukiwania przy użyciu skomplikowanych obliczeń odległości, nie są szybkie, zwykle wymagają zauważalnego ułamka sekundy. Aplikacje w czasie rzeczywistym muszą o tym wiedzieć.

Implementacja ArcGIS

Python jest preferowanym środowiskiem skryptowym dla ArcGIS (od wersji 9). Pakiet scipy.optimize ma wieloczynnikowej rootfinder rootktóry powinien robić to, co FindRootrobi w Mathematica kodu. Oczywiście sam ArcGIS oferuje dokładne obliczenia elipsoidalne odległości. Pozostała część to wszystkie szczegóły implementacji: zdecyduj, w jaki sposób zostaną uzyskane współrzędne punktu bazowego (z warstwy? Wpisanej przez użytkownika? Z pliku tekstowego? Z myszy?) Oraz w jaki sposób zostaną przedstawione dane wyjściowe (jako współrzędne wyświetlany na ekranie? jako punkt graficzny? jako nowy obiekt punktowy w warstwie?), napisz ten interfejs, przenieś pokazany tutaj kod Mathematica (proste), a wszystko będzie gotowe.

Jak zauważasz, problem ten powstaje przy określaniu granic morskich; jest często określany jako problem „trzypunktowy” i możesz to zrobić w Google i znaleźć kilka artykułów na ten temat. Jeden z tych artykułów jest przeze mnie (!) I oferuję dokładne i szybko zbieżne rozwiązanie. Patrz sekcja 14 http://arxiv.org/abs/1102.1215

Metoda składa się z następujących kroków:

1. zgadnij trzypunktowy O
2. użyj O jako środka azymutalnego równoodległego rzutu
3. projekt A, B, C, do tej projekcji
4. znajdź punkt w tej projekcji, O '
5. użyj O 'jako nowego centrum projekcji
6. powtarzaj, aż O 'i O się pokrywają

Niezbędny wzór na rozwiązanie trójpunktowe w rzucie podano w artykule. Tak długo, jak używasz dokładnej projekcji azymutalnej w jednakowej odległości, odpowiedź będzie dokładna. Konwergencja jest kwadratowa, co oznacza, że ​​potrzeba tylko kilku iteracji. To prawie na pewno przewyższy ogólne metody wyszukiwania rootów sugerowane przez @whuber.

Nie mogę bezpośrednio pomóc w ArcGIS. Możesz pobrać mój pakiet Pythona do wykonywania obliczeń geodezyjnych z https://pypi.python.org/pypi/geographiclib i kodowanie projekcji na tej podstawie jest proste.

Edytować

Problem znalezienia trójki w zdegenerowanym przypadku @ Whubera (45 + eps, 0) (45,0) (45-eps, 0) był rozważany przez Cayley'a w: O liniach geodezyjnych na sferoidie spłaszczonej , Phil. Mag. (1870), http://books.google.com/books?id=4XGIOoCMYYAC&pg=PA15

W tym przypadku tri-punkt jest uzyskiwany przez podążanie za geodezją z (45,0) z azymutem 90 i znalezienie punktu, w którym zanika skala geodezyjna. W przypadku elipsoidy WGS84 ten punkt to (-0.10690908732248, 89,89291072793167). Odległość od tego punktu do każdego z (45,001,0), (45,0), (44,999) wynosi 10010287.665788943 m (w granicach około nanometra). To o 1882 km więcej niż szacunek Whubera (który po prostu pokazuje, jak niestabilny jest ten przypadek). W przypadku kulistej ziemi trzypunkt wynosiłby oczywiście (0,90) lub (0, -90).

DODATEK: Oto implementacja metody azymutalnej w jednakowej odległości przy użyciu Matlaba

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

Testuję to za pomocą Octave

oktawa: 1> format długi
oktawa: 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

jako punkt triumfalny w Nowym Jorku, Tokio i Wellington.

Ta metoda jest niedokładna dla sąsiadujących punktów kolinearnych, np. [45,001,0], [45,0], [44,999,0]. W takim przypadku powinieneś rozwiązać dla M ₁₂ = 0 na geodezji pochodzącej z [45,0] na azymucie 90. Następująca funkcja wykonuje lewę (używając metody Newtona):

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

Na przykład daje to:

[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

Byłem ciekawy, jak szybko podejście @ cffk zbiega się w rozwiązaniu, więc napisałem test przy użyciu arcobjectów, które dały ten wynik. Odległości w metrach:

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

Oto kod źródłowy. (Edytuj) Zmieniono FindCircleCenter, aby obsługiwał przecięcia (punkty środkowe), które spadają z krawędzi projekcji azymutalnej:

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;
}

Istnieje również alternatywne podejście do wydania MSDN Magazine z czerwca 2013 r., Amoeba Method Optimization przy użyciu C # .

Edytować

W niektórych przypadkach poprzednio opublikowany kod był zbieżny z antypodem. Zmieniłem kod, aby generował ten wynik dla punktów testowych @ cffk.

Oto dane wyjściowe, które teraz wytwarza:

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

Oto poprawiony kod:

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;
    }
}

Edytować

Oto wyniki, które otrzymuję dzięki esriSRProjCS_WGS1984N_PoleAziEqui

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