Jak obliczyć promień ziemi na danej szerokości geograficznej?

19

(Widzę, że na wikipedii jest równanie, które robi dokładnie to, o co proszę, ale nie ma żadnych odniesień. Nie mam sposobu na potwierdzenie ważności tego równania!)

Rozumiem już różnicę między geocentryczną szerokością geograficzną a geodezyjną.

Zakładając, że znane są pół-duże ai pół-mniejsze, bpodano promienie. Jak obliczyć promień na danej szerokości geograficznej?

Potrzebuję jakiegoś potwierdzenia eksperta (wyprowadzenie, link do wyprowadzenia, potwierdzenie od eksperta, wyjaśnienie itp.).

Trevor Boyd Smith
źródło

Odpowiedzi:

25

To pytanie zakłada elipsoidalny model ziemi. Jego powierzchnię odniesienia uzyskuje się poprzez obrócenie elipsy wokół jej mniejszej osi (wykreślonej pionowo zgodnie z konwencją). Taka elipsa jest po prostu okręgiem, który został rozciągnięty poziomo o współczynnik a i pionowo o współczynnik b . Używając standardowej parametryzacji koła jednostki,

t --> (cos(t), sin(t))

(która definiuje cosinus i sinus), uzyskujemy parametryzację

t --> (a cos(t), b sin(t)).

(Dwa elementy tej parametryzacji opisują podróż wokół krzywej: określają, we współrzędnych kartezjańskich, naszą lokalizację w „czasie” t .)

Szerokości geodezyjnej , F , dla dowolnego punktu, to kąt, że „w górę” kieruje do płaszczyzny równikowej. Gdy a różni się od b , wartość f różni się od wartości t (z wyjątkiem wzdłuż równika i na biegunach).

Postać

Na tym zdjęciu niebieska krzywa jest jedną ćwiartką takiej elipsy (znacznie przesadzona w porównaniu z ekscentrycznością Ziemi). Czerwona kropka w lewym dolnym rogu to jej środek. Linia przerywana oznacza promień do jednego punktu na powierzchni. Jego „górny” kierunek jest pokazany za pomocą czarnego segmentu: z definicji jest on prostopadły do ​​elipsy w tym punkcie. Ze względu na przesadną mimośrodowość łatwo zauważyć, że „góra” nie jest równoległa do promienia.

W naszej terminologii t jest związany z kątem wykonanym przez promień do poziomu, a f jest kątem wykonanym przez ten czarny segment. (Zauważ, że z tej perspektywy można patrzeć na dowolny punkt na powierzchni. To pozwala nam ograniczyć zarówno t, jak i f, do wartości od 0 do 90 stopni; ich cosinus i sinus będą dodatnie, więc nie musimy się martwić o negatywne pierwiastki kwadratowe we wzorach).

Sztuką jest konwersja z parametryzacji t na jeden w kategoriach f , ponieważ w kategoriach t promień R jest łatwy do obliczenia (za pomocą twierdzenia Pitagorasa). Jego kwadrat jest sumą kwadratów składników punktu,

R(t)^2 = a^2 cos(t)^2 + b^2 sin(t)^2.

Aby dokonać tej konwersji, musimy powiązać kierunek „w górę” f z parametrem t . Ten kierunek jest prostopadły do ​​stycznej elipsy. Z definicji styczną do krzywej (wyrażoną jako wektor) uzyskuje się poprzez różnicowanie jej parametryzacji:

Tangent(t) = d/dt (a cos(t), b sin(t)) = (-a sin(t), b cos(t)).

(Różnicowanie oblicza szybkość zmian. Szybkość zmiany naszej pozycji podczas podróży po krzywej jest oczywiście naszą prędkością i zawsze wskazuje na krzywą.)

Obróć to w prawo o 90 stopni, aby uzyskać prostopadłą, zwaną wektorem „normalnym”:

Normal(t) = (b cos(t), a sin(t)).

Nachylenie tego wektora normalnego, równe (sin (t)) / (b cos (t)) („wzniesienie po przebiegu”), jest również styczną kąta, jaki robi do poziomu, skąd

tan(f) = (a sin(t)) / (b cos(t)).

Równoważnie

(b/a) tan(f) = sin(t) / cos(t) = tan(t).

(Jeśli masz dobry wgląd w geometrię euklidesową, możesz uzyskać ten związek bezpośrednio z definicji elipsy bez przechodzenia przez dowolny wyzwalacz lub rachunek różniczkowy, po prostu rozpoznając, że połączone poziome i pionowe rozwinięcia odpowiednio przez a i b powodują zmianę wszystkie zbocza tego współczynnika b / a .)

Ponownie przyjrzeć wzoru na R (t) ^ 2: wiemy, jak i B - oni określają kształt i wielkość elipsy - więc musimy tylko znaleźć cos (t) ^ 2 i sin (t) ^ 2 w kategoriach f , które powyższe równanie pozwala nam łatwo zrobić:

cos(t)^2 = 1/(1 + tan(t)^2) 
         = 1 / (1 + (b/a)^2 tan(f)^2) 
         = a^2 / (a^2 + b^2 tan(f)^2);
sin(t)^2 = 1 - cos(t)^2 
         = b^2 tan(f)^2 / (a^2 + b^2 tan(f)^2).

(Gdy tan (f) jest nieskończony, jesteśmy na biegunie, więc w takim przypadku po prostu ustaw f = t .)

To jest połączenie, którego potrzebujemy. Zamień te wartości na cos (t) ^ 2 i sin (t) ^ 2 w wyrażenie na R (t) ^ 2 i uproś, aby uzyskać

R(f)^2 = ( a^4 cos(f)^2 + b^4 sin(f)^2 ) / ( a^2 cos(f)^2 + b^2 sin(f)^2 ).

Prosta transformacja pokazuje, że to równanie jest takie samo jak w Wikipedii. Ponieważ a ^ 2 b ^ 2 = (ab) ^ 2 i (a ^ 2) ^ 2 = a ^ 4,

R(f)^2 = ( (a^2 cos(f))^2 + (b^2 sin(f))^2 ) / ( (a cos(f))^2 + (b sin(f))^2 )
Whuber
źródło
+1 .. poza tym, że sądzę, że ostateczna formuła ma paren nie na miejscu ... nie należy (b^4 sin(f))^2go zmieniać na (b^4 sin(f)^2)?
Kirk Kuykendall
naprawdę się cieszę, że są wokół nas eksperci =).
Trevor Boyd Smith
Czy plik Geogebra (HTML) może zostać opublikowany na tej stronie? Mam promień pionowy, który może wizualnie pokazać, co się dzieje.
Możesz wyeksportować oryginał w formacie .png, @Dan: użyj okna dialogowego Plik | Eksportuj. Zalecam używanie dużych czcionek (16 lub 18 punktów wydaje się działać dobrze) i powiększanie obrazu tak daleko, jak to możliwe.
whuber
Zakładam, że wtedy interaktywność zostanie utracona. Demo pokazuje, jak zróżnicowanie promieni i szerokości geograficznej zainteresowania zmienia właściwości.
3

Ciekawe, że moje niepiśmienne rozwiązanie matematyczne wykonało zadanie po 5 minutach namysłu i kodowania, czy nie należy brać pod uwagę czynnika spłaszczania zamiast idealnego modelu eliptycznego?

        double pRad = 6356.7523142;
        double EqRad = 6378.137;                      
        return pRad + (90 - Math.Abs(siteLatitude)) / 90 * (EqRad - pRad); 
Howard Grover
źródło
1
Gdzie pRad jest promieniem biegunowym, a EqRad jest promieniem równikowym.
Stefan Steiger,
to jedyna odpowiedź, którą mogłem przeczytać. Wydaje mi się, że to działa.
Sean Bradley,
1
Widzę, że wykonujesz liniową interpolację promienia między biegunem a równikiem. Chociaż nie ma powodu, aby sądzić, że interpolacja liniowa jest dokładna , użyję tego jako „wystarczająco dobrego” dla Ziemi, biorąc pod uwagę jej łagodny współczynnik spłaszczenia. BTW Myślę, że nieco łatwiej jest odczytać odpowiednik:, return E + (P - E) * Abs(Lat) / 90więc nie musisz mieć tego 90 - ...w formule.
ToolmakerSteve
2

wprowadź opis zdjęcia tutaj

Przynajmniej taką formułę znalazłem w amerykańskim Centrum analizy i oceny danych (DAAC) dla wiki Departamentu Obrony (DoD) High Performance Computing Modernization Program (HPCMP) wiki . Mówi się, że mocno pożyczyli od wpisu Wikipedii . Mimo to fakt, że zachowali tę formułę, powinien na coś liczyć.

RK
źródło
Czy możesz podać link do treści?
Trevor Boyd Smith
gdzie φ jest szerokością geodezyjną, a a (pół-główna oś) i b (pół-mniejsza oś) to odpowiednio promień równikowy i promień biegunowy. var a = 6378137; // m var b = 6356752.3142; // m en.wikipedia.org/wiki/Semi-major_and_semi-minor_axes en.wikipedia.org/wiki/World_Geodetic_System
Stefan Steiger