¿Entender los términos en la fórmula de la longitud del grado?

13

Las calculadoras en línea como http://www.csgnetwork.com/degreelenllavcalc.html (ver fuente de la página) usan las fórmulas a continuación para obtener metros por grado. En general, entiendo cómo la distancia por grado varía según la ubicación de la latitud, pero no entiendo cómo se traduce eso a continuación. Más específicamente, ¿de dónde provienen las constantes, los 3 términos "cos" en cada fórmula y los coeficientes (2, 4, 6; 3 y 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));
convert latitude-longitude gis-principle length Brent
fuente
3
En un círculo, los términos de la forma cos (m * x) para m = 0, 1, 2, ... juegan el mismo papel que los monomios 1, x, x ^ 2, x ^ 3, ... para Taylor serie en la línea. Cuando vea una expansión de este tipo, puede pensar en ella de la misma manera: cada término da una aproximación de orden superior a una función. Por lo general, tales series trigonométricas son infinitas; pero en el uso práctico pueden truncarse tan pronto como el error de aproximación sea aceptable. Parte de esta tecnología se encuentra bajo el capó de cada SIG porque muchas proyecciones esferoidales se calculan utilizando tales series.
whuber
Esto es muy útil para calcular distancias donde la distancia entre líneas de latitud varía, también útil para ayudar a determinar dónde trazar puntos en un mapa de mercator si tiene una cuadrícula x, y como superposición
Consejo: no se olvide de usar radianes para lat(a pesar de las variables resultantes latleny longlenson en metros por grado, no metros por radián). Si usa grados para lat, incluso puede terminar con un valor negativo para longlen.
Luke Hutchison el

Respuestas:

23

El radio principal del esferoide WGS84 es a = 6378137 metros y su aplanamiento inverso es f = 298.257223563, de donde la excentricidad cuadrada es

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

El radio de curvatura meridional en la latitud phi es

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

y el radio de curvatura a lo largo del paralelo es

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

Además, el radio del paralelo es 

r = N cos(phi)

Estas son correcciones multiplicativas de los valores esféricos de M y N , que son iguales al radio esférico a , que es a lo que se reducen cuando e2 = 0.

Figura

En el punto amarillo a 45 grados de latitud norte, el disco azul de radio M es el círculo de osculación ("beso") en la dirección del meridiano y el disco rojo de radio N es el círculo de osculación en la dirección del paralelo: ambos los discos contienen la dirección "hacia abajo" en este punto. Esta cifra exagera el aplanamiento de la tierra en dos órdenes de magnitud.

Los radios de curvatura determinar las longitudes de grados: cuando un círculo tiene un radio de R , su perímetro de longitud 2 pi R cubiertas 360 grados, de donde la longitud de un grado es pi * R / 180. Sustituyendo M y r para R - - es decir, multiplicando M y r por pi / 180 - da simples exactas fórmulas para las longitudes de grado.

Estas fórmulas, que se basan únicamente en los valores dados de a y f (que se pueden encontrar en muchos lugares ) y la descripción del esferoide como un elipsoide de rotación, concuerdan con los cálculos en la pregunta dentro de 0.6 partes por millones (unos pocos centímetros), que es aproximadamente el mismo orden de magnitud de los coeficientes más pequeños en la pregunta, lo que indica que están de acuerdo. (La aproximación siempre es un poco baja.) En el gráfico, el error relativo en la longitud de un grado de latitud es negro y el de la longitud está punteado en rojo:

Figura

En consecuencia, podemos entender que los cálculos en la pregunta son aproximaciones (a través de series trigonométricas truncadas) a las fórmulas dadas anteriormente.

Los coeficientes se pueden calcular a partir de la serie del coseno de Fourier para M y r como funciones de la latitud. Se dan en términos de funciones elípticas de e2, lo que sería demasiado complicado para reproducir aquí. Para el esferoide WGS84, mis cálculos dan

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

(Puede adivinar cómo p4ingresa la fórmula. :) La cercanía de estos valores a los parámetros en el código atestigua la exactitud de esta interpretación. Esta aproximación mejorada es precisa mucho mejor que una parte por mil millones en todas partes.

Para probar esta respuesta, ejecuté un Rcódigo para llevar a cabo ambos 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")

El cálculo exacto con radiise puede utilizar para imprimir tablas de longitudes de grados, como en

zapsmall(radii(phi) * pi / 180)

La salida está en metros y tiene este aspecto (con algunas líneas eliminadas):

          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

Referencias

LM Bugayevskiy y JP Snyder, Proyecciones de mapas: un manual de referencia. Taylor y Francis, 1995. (Apéndice 2 y Apéndice 4)

JP Snyder, Proyecciones de mapas: un manual de trabajo. USGS Professional Paper 1395, 1987. (Capítulo 3)

whuber
fuente
No sé por qué se usaría una aproximación tan complicada a un simple par de fórmulas ...
whuber
¡Qué minuciosa, excelente respuesta! Parece correcto; ahora solo necesito repasar esta matemática para entenderla. :)
Brent
@Brent He agregado una figura para ayudarte a entender las matemáticas.
whuber
0

Esa es la fórmula de Haversine , aunque expresada de manera extraña.

tmcw
fuente
¡Claramente no es la fórmula de Haversine! Esto es (relacionado con) una perturbación del mismo utilizado para el esferoide. Ni siquiera encuentra distancias entre pares de puntos arbitrarios, que es para lo que se usa la fórmula de Haversine (en la esfera).
whuber
1
En otras palabras, la fórmula de Haversine calcula la distancia del gran círculo, y esta fórmula es una perturbación de la misma que calcula la distancia elipsoide más precisa.
Brent