Como você calcula o raio da Terra em uma dada latitude geodésica?

19

(Vejo que há uma equação na Wikipedia que faz exatamente o que estou pedindo, mas não há referências. Não tenho como confirmar a validade dessa equação!)

Eu já entendo a diferença entre Latitude Geocêntrica vs Latitude Geodésica.

Supondo que sejam conhecidos os raios semi-maior ae semi-menor conhecidos b. Como você calcula o raio em uma determinada latitude geodésica?

Preciso de algum tipo de confirmação de especialista (derivação, link para derivação, confirmação de especialista, explicação etc.).

Trevor Boyd Smith
fonte

Respostas:

25

Esta questão assume um modelo elipsoidal da terra. Sua superfície de referência é obtida girando uma elipse em torno de seu eixo menor (plotado verticalmente por convenção). Tal elipse é apenas um círculo que foi esticado horizontalmente por um fator de a e verticalmente por um fator de b . Usando a parametrização padrão do círculo unitário,

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

(que define cosseno e seno), obtemos uma parametrização

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

(Os dois componentes dessa parametrização descrevem uma viagem em torno da curva: eles especificam, em coordenadas cartesianas, nossa localização no "tempo" t .)

A latitude geodésica , f , de qualquer ponto é o ângulo que "para cima" faz com o plano equatorial. Quando a difere de b , o valor de f difere do valor de t (exceto ao longo do equador e nos polos).

Figura

Nesta figura, a curva azul é um quadrante dessa elipse (muito exagerado em comparação com a excentricidade da Terra). O ponto vermelho no canto inferior esquerdo é o centro. A linha tracejada designa o raio para um ponto na superfície. Sua direção "para cima" é mostrada com um segmento preto: é, por definição, perpendicular à elipse naquele ponto. Devido à excentricidade exagerada, é fácil ver que "para cima" não é paralelo ao raio.

Em nossa terminologia, t está relacionado ao ângulo feito pelo raio em relação à horizontal ef é o ângulo feito por esse segmento preto. (Note-se que qualquer . Ponto da superfície pode ser visto a partir desta perspectiva Isso nos permite limitar tanto t e f a mentira entre 0 e 90 graus, os seus co-senos e senos será positivo, de modo que não precisa se preocupar com negativa raízes quadradas nas fórmulas.)

O truque é converter da parametrização t para uma em termos de f , porque em termos de t é fácil calcular o raio R (através do teorema de Pitágoras). Seu quadrado é a soma dos quadrados dos componentes do ponto,

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

Para fazer essa conversão, precisamos relacionar a direção "para cima" f com o parâmetro t . Essa direção é perpendicular à tangente da elipse. Por definição, uma tangente a uma curva (expressa como vetor) é obtida diferenciando sua parametrização:

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

(A diferenciação calcula a taxa de mudança. A taxa de mudança de nossa posição à medida que percorremos a curva é, obviamente, a nossa velocidade , e que sempre aponta ao longo da curva.)

Gire no sentido horário 90 graus para obter a perpendicular, chamada de vetor "normal":

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

A inclinação desse vetor normal, igual a (a sen (t)) / (b cos (t)) ("subida sobre corrida"), também é a tangente do ângulo que ele faz em relação à horizontal, de onde

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

Equivalentemente,

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

(Se tiver boa perspectiva geometria euclidiana, pode-se obter este relacionamento directamente da definição de uma elipse, sem passar por qualquer trig ou cálculo, simplesmente através do reconhecimento de que as expansões horizontais e verticais combinados por um e b , respectivamente, têm o efeito de alterar todas as inclinações por esse fator b / a .)

Olhe novamente para a fórmula para R (t) ^ 2: sabemos a e b - eles determinam a forma eo tamanho da elipse - por isso só precisa encontrar cos (t) ^ 2 e sin (t) ^ 2 em termos de f , que a equação anterior nos permite fazer facilmente:

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).

(Quando tan (f) é infinito, estamos no polo, então basta definir f = t nesse caso.)

Essa é a conexão que precisamos. Substitua esses valores por cos (t) ^ 2 e sin (t) ^ 2 na expressão para R (t) ^ 2 e simplifique para obter

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

Uma transformação simples mostra que essa equação é a mesma encontrada na Wikipedia. Como a ^ 2 b ^ 2 = (ab) ^ 2 e (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
fonte
+1 .. exceto que eu acho que a fórmula final tem um ponto fora do lugar ... não deve (b^4 sin(f))^2ser alterado para (b^4 sin(f)^2)?
22412 Kirk Kuykendall
realmente feliz que haja alguns especialistas por aqui =).
Trevor Boyd Smith
Um arquivo Geogebra (html) pode ser publicado neste site? Eu tenho um raio da vertical principal que pode demonstrar visualmente o que está acontecendo.
Você pode exportar o original no formato .png, @Dan: use a caixa de diálogo Arquivo | Exportar. Eu recomendo o uso de fontes grandes (16 ou 18 pontos parecem funcionar bem) e aplicar o zoom o máximo possível na imagem.
whuber
Eu assumo que a interatividade será perdida então. A demonstração demonstra como a variação dos raios e da latitude de interesse altera as propriedades.
3

Interessante descobrir que minha solução analfabeta de matemática fez o trabalho com 5 minutos de reflexão e codificação, não seria necessário considerar o fator de achatamento em vez de um modelo elíptico perfeito?

        double pRad = 6356.7523142;
        double EqRad = 6378.137;                      
        return pRad + (90 - Math.Abs(siteLatitude)) / 90 * (EqRad - pRad); 
Howard Grover
fonte
1
Onde pRad é o raio polar e EqRad é o raio equatorial.
Stefan Steiger
esta é a única resposta que eu pude ler. Parece funcionar para mim.
Sean Bradley
1
Vejo que você está fazendo uma interpolação linear de raio entre o polo e o equador. Embora não haja razão para acreditar que a interpolação linear seja precisa , usarei isso como "bom o suficiente" para a Terra, dado seu leve fator de achatamento. BTW eu acho que é um pouco mais fácil ler o equivalente:, return E + (P - E) * Abs(Lat) / 90então não precisa ter 90 - ...na fórmula.
Home
2

insira a descrição da imagem aqui

Pelo menos essa é a fórmula encontrada no Centro de Análise e Avaliação de Dados dos EUA (DAAC) para o wiki do Programa de Modernização da Computação de Alto Desempenho (HPCMP) do Departamento de Defesa (DoD) . Diz que eles pegaram emprestado pesadamente da entrada da Wikipedia . Ainda assim, o fato de manterem essa fórmula deve valer alguma coisa.

RK
fonte
Você pode fornecer um link para o conteúdo?
Trevor Boyd Smith
onde φ é a latitude geodésica e a (eixo semi-maior) eb (eixo semi-menor) são, respectivamente, o raio equatorial e o raio polar. 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