Zrozumienie terminów we wzorze Długość stopnia?

13

Kalkulatory online, takie jak http://www.csgnetwork.com/degreelenllavcalc.html (zobacz źródło strony), wykorzystują poniższe formuły, aby uzyskać metry na stopień. Rozumiem ogólnie, jak odległość na stopień różni się w zależności od położenia szerokości geograficznej, ale nie rozumiem, jak to przekłada się na poniższe. Mówiąc dokładniej, skąd biorą się stałe, terminy 3 „cos” w każdej formule oraz współczynniki (2, 4, 6; 3 i 5) dla „lat”?

    // Set up "Constants"
    m1 = 111132.92;     // latitude calculation term 1
    m2 = -559.82;       // latitude calculation term 2
    m3 = 1.175;         // latitude calculation term 3
    m4 = -0.0023;       // latitude calculation term 4
    p1 = 111412.84;     // longitude calculation term 1
    p2 = -93.5;         // longitude calculation term 2
    p3 = 0.118;         // longitude calculation term 3

    // Calculate the length of a degree of latitude and longitude in meters
    latlen = m1 + (m2 * Math.cos(2 * lat)) + (m3 * Math.cos(4 * lat)) +
            (m4 * Math.cos(6 * lat));
    longlen = (p1 * Math.cos(lat)) + (p2 * Math.cos(3 * lat)) +
                (p3 * Math.cos(5 * lat));
Brent
źródło
3
Na kole warunki postaci cos (m * x) dla m = 0, 1, 2, ... grają tę samą rolę, co monomiały 1, x, x ^ 2, x ^ 3, ... działają dla Taylora seria na linii. Kiedy widzisz tego rodzaju rozwinięcie, możesz myśleć o tym w ten sam sposób: każdy termin daje przybliżenie funkcji wyższego rzędu. Zwykle takie serie trygonometryczne są nieskończone; ale w praktyce można je obciąć, gdy tylko błąd przybliżenia będzie akceptowalny. Niektóre takie technologie leżą pod maską każdego GIS, ponieważ wiele rzutów sferoidalnych jest obliczanych przy użyciu takich serii.
whuber
Jest to bardzo przydatne do obliczania odległości, w których odległość między liniami szerokości jest różna, a także przydatne, aby pomóc określić, gdzie
Wskazówka: Nie zapomnij użyć radianów na lat(choć powstałych zmiennych latleni longlensą w metrach na stopniu nie metrów na radian). Jeśli użyjesz stopni dla lat, możesz nawet otrzymać ujemną wartość dla longlen.
Luke Hutchison

Odpowiedzi:

23

Główny promień sferoidy WGS84 wynosi a = 6378137 metrów, a jej odwrotne spłaszczenie wynosi f = 298.257223563, skąd kwadratowa mimośrodowość wynosi

e2 = (2 - 1/f)/f = 0.0066943799901413165.

Promień południkowy krzywizny na szerokości geograficznej phi wynosi

M = a(1 - e2) / (1 - e2 sin(phi)^2)^(3/2)

a promień krzywizny wzdłuż równoległości wynosi

N = a / (1 - e2 sin(phi)^2)^(1/2)

Ponadto promień równoległości wynosi

r = N cos(phi)

Są to multiplikatywne korekty sferycznych wartości M i N , z których oba są równe promieniu sferycznemu a , do czego redukują się, gdy e2 = 0.

Postać

W żółtym punkcie na 45 stopniach szerokości geograficznej północnej niebieski dysk o promieniu M jest oscylującym („całującym”) okręgiem w kierunku południka, a czerwony dysk o promieniu N jest oscylującym okręgiem w kierunku równoległym: oba dyski zawierają w tym momencie kierunek „w dół”. Liczba ta przesadza spłaszczanie ziemi o dwa rzędy wielkości.

Promienie krzywizny określenia długości stopni: jeżeli koło ma promień R , jego obwód długości 2 pi R pokrywy 360 stopni, z którego długość jednego stopnia jest pi * R / 180 Podstawiając M i R do R - - to znaczy pomnożenie M i r przez pi / 180 - daje proste dokładne wzory na długości stopni.

Wzory te - oparte wyłącznie na podanych wartościach a i f (które można znaleźć w wielu miejscach ) i opisie sferoidy jako elipsoidy obrotu - zgadzają się z obliczeniami w pytaniu z dokładnością do 0,6 części na milion (kilka centymetrów), co jest w przybliżeniu tym samym rzędem wielkości najmniejszych współczynników w pytaniu, co oznacza, że ​​się zgadzają. (Przybliżenie jest zawsze trochę niskie.) Na wykresie błąd względny długości stopnia szerokości geograficznej jest czarny, a błąd długości geograficznej jest przerywany na czerwono:

Postać

W związku z tym możemy zrozumieć, że obliczenia w pytaniu są przybliżeniami (poprzez skrócone szeregi trygonometryczne) do podanych powyżej wzorów.


Współczynniki można obliczyć z szeregu cosinus Fouriera dla M i r jako funkcji szerokości geograficznej. Podano je w kategoriach eliptycznych funkcji e2, które byłyby zbyt nieporządne, aby się tutaj odtworzyć. Dla sferoidy WGS84 moje obliczenia podają

  m1 = 111132.95255
  m2 = -559.84957
  m3 = 1.17514
  m4 = -0.00230
  p1 = 111412.87733
  p2 = -93.50412
  p3 = 0.11774
  p4 = -0.000165

(Można się domyślać, jak p4wchodzi do formuły. :) Bliskość tych wartości do parametrów w kodzie świadczy o poprawności tej interpretacji. To ulepszone przybliżenie jest dokładne wszędzie o wiele lepiej niż jedna część na miliard.


Aby przetestować tę odpowiedź, wykonałem Rkod, aby wykonać oba obliczenia:

#
# Radii of meridians and parallels on a spheroid.  Defaults to WGS84 meters.
# Input is latitude (in degrees).
#
radii <- function(phi, a=6378137, e2=0.0066943799901413165) {
  u <- 1 - e2 * sin(phi)^2
  return(cbind(M=(1-e2)/u, r=cos(phi)) * (a / sqrt(u))) 
}
#
# Approximate calculation.  Same interface (but no options).
#
m.per.deg <- function(lat) {
  m1 = 111132.92;     # latitude calculation term 1
  m2 = -559.82;       # latitude calculation term 2
  m3 = 1.175;         # latitude calculation term 3
  m4 = -0.0023;       # latitude calculation term 4
  p1 = 111412.84;     # longitude calculation term 1
  p2 = -93.5;         # longitude calculation term 2
  p3 = 0.118;         # longitude calculation term 3

  latlen = m1 + m2 * cos(2 * lat) + m3 * cos(4 * lat) + m4 * cos(6 * lat);
  longlen = p1 * cos(lat) + p2 * cos(3 * lat) + p3 * cos(5 * lat);
  return(cbind(M.approx=latlen, r.approx=longlen))
}
#
# Compute the error of the approximation `m.per.deg` compared to the 
# correct formula and plot it as a function of latitude.
#
phi <- pi / 180 * seq(0, 90, 10)
names(phi) <- phi * 180 / pi
matplot(phi * 180 / pi, 10^6 * ((m.per.deg(phi) - radii(phi) * pi / 180) / 
       (radii(phi) * pi / 180)),
        xlab="Latitude (degrees)", ylab="Relative error * 10^6",lwd=2, type="l")

Dokładne obliczenia za pomocą radiimogą być używane do drukowania tabel o długości stopni, jak w

zapsmall(radii(phi) * pi / 180)

Dane wyjściowe są w metrach i wyglądają tak (z usuniętymi niektórymi liniami):

          M         r
0  110574.3 111319.49
10 110607.8 109639.36
20 110704.3 104647.09
...
80 111659.9  19393.49
90 111694.0      0.00

Bibliografia

LM Bugayevskiy i JP Snyder, Projekcje map - Podręcznik referencyjny. Taylor i Francis, 1995. (Załącznik 2 i Załącznik 4)

JP Snyder, Projekcje map - Podręcznik roboczy. USGS Professional Paper 1395, 1987. (Rozdział 3)

Whuber
źródło
Nie wiem, dlaczego kiedykolwiek miałoby być tak skomplikowane zbliżenie do prostej pary wzorów ...
whuber
Cóż za dokładna, doskonała odpowiedź! Wydaje się poprawne; teraz muszę tylko odświeżyć matematykę, aby ją zrozumieć. :)
Brent,
@Brent Dodałem rysunek, który pomoże ci zrozumieć matematykę.
whuber
0

To formuła Haversine , choć wyrażona w dziwny sposób.

tmcw
źródło
To wyraźnie nie jest formuła Haversine! Jest to (związane z) zaburzeniem stosowanym przez sferoidę. Nie znajduje nawet odległości między dowolnymi parami punktów, do czego używana jest formuła Haversine (na kuli).
whuber
1
Innymi słowy, wzór Haversine oblicza odległość wielkiego koła, a ta formuła jest zaburzeniem, które oblicza dokładniejszą odległość elipsoidy?
Brent,