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.
