Każde oprogramowanie, które może dokładnie rzutować współrzędne, może obliczyć dokładne wskaźniki Tissot .
Dobrym źródłem formuł jest Snyder, John, Map Projections - A Working Manual , głównie na s. 20–26. (Nie będę ich tutaj reprodukować, ponieważ ta strona nie ma odpowiednich narzędzi do przekazywania wzorów matematycznych). Wymagają one wszystkich czterech pierwszych pochodnych rzutowanych współrzędnych (x, y) w odniesieniu do współrzędnych sferycznych (lat, lon) = (phi, lambda):
dx / d(phi), dx / d(lambda);
dy / d(phi), dy / d(lambda).
Cała reszta w TI jest obliczana w kategoriach tych (przy użyciu niektórych funkcji arytmetycznych i trygonometrycznych: cosinus, sinus główny odwrotny i główna styczna odwrotna). Obliczenia wymagają opisu kształtu ziemi. Aby uzyskać największą dokładność, należy zastosować elipsoidalny układ odniesienia z półosiową osią a i mimośrodowością e. (Będą one znane oprogramowaniu.)
Książka Snydera zawiera instrukcje, jak obliczyć wszystko oprócz tych pochodnych. Zrób to numerycznie. Miałem doskonałe wyniki, stosując szacunki centralnej różnicy skończonych pierwszego rzędu w odległości h = 10 ^ (- 5,2) radianów (zwykle około 50 metrów): jest to dobry kompromis między próbą zbliżenia się do nieskończenie bliskiego i utratą zbyt dużej precyzji od zaokrąglenie zmiennoprzecinkowe (przy założeniu podwójnej precyzji), ponieważ popełniony błąd jest proporcjonalny do (10 ^ (- 5.2)) ^ 2 = 10 ^ (- 10.4), a 10 ^ (- 5.2) jest równy 10 ^ 10.4 razy dokładność podwójnej precyzji IEEE z 10 ^ (- 15,6) i wciąż jest znacznie większa niż typowa precyzja w rzutach, które zwykle wynoszą od 10 ^ (- 10) do około 10 ^ (- 14).
Jak zatem obliczyć szacunki różnic skończonych? Ta część jest zaskakująco łatwa. Aby uzyskać dx / d (phi) w punkcie (phi, lambda), poproś GIS o rzutowanie punktów
(phi - h/2, lambda) --> (x0,y0),
(phi + h/2, lambda) --> (x1,y1).
Użyj szacunków
dx / d(phi) = (x1 - x0)/h,
dy / d(phi) = (y1 - y0)/h.
Podobnie rzutuj punkty
(phi, lambda - h/2) --> (x2,y2),
(phi, lambda + h/2) --> (x3,y3)
i użyj szacunków
dx / d(lambda) = (x3 - x2)/h,
dy / d(lambda) = (y3 - y2)/h.
To wymaga czterech projekcji i odrobiny arytmetyki. (Możesz zredukować go do trzech, stosując różnice niecentralne, ale dokładność nieco się obniża. Mądrze jest dążyć do wysokiej dokładności, nie pozwalając, by h stał się zbyt mały, chyba że masz pewność, że twój GIS używa oceny ankietowej (milimetr) dokładność w formułach projekcyjnych.)
Z tych pochodnych, wraz ze wzorami Snydera (zwracając uwagę na modyfikacje opisane w 4-19 i 4-21), można uzyskać długości osi wskaźnika Tissot w (phi, lambda) i jego orientacji. Na mapach na skalę światową TI będzie tak mały, że będzie niewidoczny, więc ostateczną rzeczą do zrobienia jest decyzja, ile chcesz przeskalować każdego TI. Współczynnik skali określam, sprawdzając, jak duża będzie mapa, znajdując rozmiary typowych wskaźników TI na mapie i skalując, aby te wskaźniki miały około 6% szerokości tak jak mapa. W każdym razie to dobry początek; Pozwalam użytkownikowi dostosować rozmiar TI. Oczywiście przeskalujesz wszystkie TI o tę samą wartość, aby można je było porównać, a każdy z nich zostanie przeskalowany wokół własnego centrum (uzyskanego w piątej projekcji (phi, lambda) -> (x, y) ).
Miłym dodatkiem do eliptycznego przedstawienia TI jest pokazanie kierunków lokalnego południka i równoległości: wtedy, na pierwszy rzut oka, możesz ocenić zbieżność siatki . Pokazuję również standardowy okrąg (reprezentujący brak zniekształceń) koncentryczny z każdym TI, ponieważ poprawia on zdolność czytelnika do oceny stopnia zniekształcenia reprezentowanego przez każdą elipsę.
Godne uwagi w tej projekcji Mollweide jest ekstremalne TI w pobliżu bieguna południowego. Jest to nadal idealna elipsa i dokładnie opisuje zniekształcenie mapy.
