Compreendendo os termos da fórmula Comprimento do grau?

13

Calculadoras online, como http://www.csgnetwork.com/degreelenllavcalc.html (ver código-fonte da página), usam as fórmulas abaixo para obter medidores por grau. Eu entendo em geral como a distância por grau varia dependendo da localização da latitude, mas não entendo como isso se traduz no abaixo. Mais especificamente, de onde vêm as constantes, os 3 termos "cos" em cada fórmula e os coeficientes (2, 4, 6; 3 e 5) para "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
fonte
3
Em um círculo, os termos da forma cos (m * x) para m = 0, 1, 2, ... desempenham o mesmo papel que os monômios 1, x, x ^ 2, x ^ 3, ... desempenham por Taylor série na linha. Quando você vê uma expansão desse tipo, pode pensar da mesma maneira: cada termo fornece uma aproximação de ordem superior a uma função. Geralmente essas séries trigonométricas são infinitas; mas, na prática, eles podem ser truncados assim que o erro de aproximação for aceitável. Algumas dessas tecnologias estão ocultas em todos os SIGs porque muitas projeções esferoidais são calculadas usando essas séries.
whuber
Isto é muito útil para distâncias de cálculo onde a distância entre as linhas de latitude variam, também úteis para ajudar a determinar onde plotar pontos em um mapa mercator se você tiver um x, y grade como uma sobreposição
Dica: não se esqueça de usar radianos para lat(mesmo que as variáveis resultantes latlene longlenestão em metros por grau, não metros por radiano). Se você usa graus para lat, pode até acabar com um valor negativo para longlen.
Luke Hutchison

Respostas:

23

O raio principal do esferóide WGS84 é a = 6378137 metros e seu achatamento inverso é f = 298.257223563, de onde a excentricidade ao quadrado é

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

O raio meridional da curvatura na latitude phi é

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

e o raio de curvatura ao longo do paralelo é

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

Além disso, o raio do paralelo é

r = N cos(phi)

Essas são correções multiplicativas aos valores esféricos de M e N , ambos iguais ao raio esférico a , que é o que eles reduzem quando e2 = 0.

Figura

No ponto amarelo a 45 graus de latitude norte, o disco azul do raio M é o círculo osculante ("beijo") na direção do meridiano e o disco vermelho do raio N é o círculo osculante na direção do paralelo: ambos discos contêm a direção "para baixo" neste momento. Esta figura exagera o achatamento da Terra em duas ordens de magnitude.

Os raios de curvatura determinam os comprimentos dos graus: quando um círculo tem um raio de R , seu perímetro de comprimento 2 pi R cobre 360 ​​graus, onde o comprimento de um grau é pi * R / 180. Substituindo M e r por R - - isto é, multiplicar M e r por pi / 180 - fornece fórmulas exatas simples para os comprimentos de graus.

Estas fórmulas - os quais são baseados unicamente nos valores dados de uma e f (que pode ser encontrado em muitos lugares ) e a descrição do esferóide como um elipsóide de rotação - de acordo com os cálculos na pergunta para dentro de 0,6 partes por milhões (alguns centímetros), que é aproximadamente a mesma ordem de magnitude dos menores coeficientes da pergunta, indicando que eles concordam. (A aproximação é sempre um pouco baixa.) No gráfico, o erro relativo no comprimento de um grau de latitude é preto e o de longitude é tracejado em vermelho:

Figura

Dessa forma, podemos entender que os cálculos na questão são aproximações (através de séries trigonométricas truncadas) às fórmulas fornecidas acima.


Os coeficientes podem ser calculados a partir da série de cosseno de Fourier para M e r como funções da latitude. Eles são dados em termos de funções elípticas de e2, o que seria muito confuso para se reproduzir aqui. Para o esferóide WGS84, meus cálculos dão

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

(Você pode adivinhar como p4entra a fórmula. :) A proximidade desses valores com os parâmetros no código atesta a exatidão dessa interpretação. Essa aproximação aprimorada é precisa para muito melhor do que uma parte por bilhão em qualquer lugar.


Para testar esta resposta, executei o Rcódigo para realizar os dois cálculos:

#
# 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")

O cálculo exato com radiipode ser usado para imprimir tabelas com comprimentos de graus, como em

zapsmall(radii(phi) * pi / 180)

A saída está em metros e é assim (com algumas linhas removidas):

          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

Referências

LM Bugayevskiy e JP Snyder, Projeções de Mapas - Um Manual de Referência. Taylor e Francis, 1995. (Apêndice 2 e Apêndice 4)

JP Snyder, Projeções de Mapas - Um Manual de Trabalho. USGS Professional Paper 1395, 1987. (Capítulo 3)

whuber
fonte
Não sei por que uma aproximação tão complicada a um simples par de fórmulas seria usada ....
whuber
Que resposta excelente e completa! Parece correto; agora eu só preciso revisar essa matemática para entender. :)
Brent
@Brent Adicionei uma figura para ajudá-lo a entender a matemática.
whuber
0

Essa é a fórmula de Haversine , embora expressa de uma maneira estranha.

tmcw
fonte
Claramente não é a fórmula de Haversine! Isso é (relacionado a) uma perturbação usada para o esferóide. Ele nem encontra distâncias entre pares de pontos arbitrários, e é para isso que a fórmula de Haversine é usada (na esfera).
whuber
1
Em outras palavras, a fórmula de Haversine calcula a distância do grande círculo, e essa fórmula é uma perturbação dela que calcula a distância elipsóide mais precisa?
Brent