¿Cómo se calcula el radio de la Tierra en una latitud geodésica dada?

19

(Veo que hay una ecuación en wikipedia que hace exactamente lo que estoy pidiendo, pero no hay referencias. ¡No tengo forma de confirmar la validez de esta ecuación!)

Ya entiendo la diferencia entre la latitud geocéntrica y la latitud geodésica.

Suponiendo que se dan radios semi mayores ay semi menores conocidos b. ¿Cómo se calcula el radio en una latitud geodésica dada?

Necesito algún tipo de confirmación experta (derivación, enlace a derivación, confirmación del experto, explicación, etc.).

Trevor Boyd Smith
fuente

Respuestas:

25

Esta pregunta supone un modelo elipsoidal de la tierra. Su superficie de referencia se obtiene girando una elipse alrededor de su eje menor (trazada verticalmente por convención). Tal elipse es solo un círculo que se ha extendido horizontalmente por un factor de a y verticalmente por un factor de b . Usando la parametrización estándar del círculo unitario,

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

(que define coseno y seno), obtenemos una parametrización

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

(Los dos componentes de esta parametrización describen un viaje alrededor de la curva: especifican, en coordenadas cartesianas, nuestra ubicación en el "tiempo" t ).

La latitud geodésica , f , de cualquier punto es el ángulo que "hacia arriba" forma el plano ecuatorial. Cuando a difiere de b , el valor de f difiere del de t (excepto a lo largo del ecuador y en los polos).

Figura

En esta imagen, la curva azul es un cuadrante de dicha elipse (muy exagerada en comparación con la excentricidad de la tierra). El punto rojo en la esquina inferior izquierda es su centro. La línea discontinua designa el radio a un punto en la superficie. Su dirección "hacia arriba" se muestra con un segmento negro: es, por definición, perpendicular a la elipse en ese punto. Debido a la excentricidad exagerada, es fácil ver que "arriba" no es paralelo al radio.

En nuestra terminología, t está relacionado con el ángulo formado por el radio con respecto a la horizontal yf es el ángulo formado por ese segmento negro. (Tenga en cuenta que cualquier punto de la superficie se puede ver desde esta perspectiva. Esto nos permite limitar que t y f se encuentren entre 0 y 90 grados; sus cosenos y senos serán positivos, por lo que no tenemos que preocuparnos por los negativos raíces cuadradas en las fórmulas.)

El truco consiste en convertir de la parametrización t a una en términos de f , porque en términos de t el radio R es fácil de calcular (a través del teorema de Pitágoras). Su cuadrado es la suma de cuadrados de los componentes del punto,

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

Para realizar esta conversión, necesitamos relacionar la dirección "hacia arriba" f con el parámetro t . Esta dirección es perpendicular a la tangente de la elipse. Por definición, se obtiene una tangente a una curva (expresada como un vector) diferenciando su parametrización:

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

(La diferenciación calcula la tasa de cambio. La tasa de cambio de nuestra posición a medida que viajamos alrededor de la curva es, por supuesto, nuestra velocidad , y eso siempre apunta a lo largo de la curva).

Gire esto en el sentido de las agujas del reloj 90 grados para obtener el perpendicular, llamado vector "normal":

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

La pendiente de este vector normal, igual a (a sin (t)) / (b cos (t)) ("aumento sobre la carrera"), es también la tangente del ángulo que hace a la horizontal, de donde

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

Equivalentemente

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

(Si tiene una buena visión de la geometría euclidiana, podría obtener esta relación directamente de la definición de una elipse sin pasar por ningún trigonométrico o cálculo, simplemente reconociendo que las expansiones horizontales y verticales combinadas por a y b respectivamente tienen el efecto de cambiar todas las pendientes por este factor b / a .)

Mire nuevamente la fórmula para R (t) ^ 2: sabemos a y b - ellos determinan la forma y el tamaño de la elipse - por lo que solo necesitamos encontrar cos (t) ^ 2 y sin (t) ^ 2 en términos de f , que la ecuación anterior nos permite hacer fácilmente:

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

(Cuando tan (f) es infinito, estamos en el polo, así que solo establece f = t en ese caso).

Esta es la conexión que necesitamos. Sustituya estos valores por cos (t) ^ 2 y sin (t) ^ 2 en la expresión para R (t) ^ 2 y simplifique para obtener

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

Una transformación simple muestra que esta ecuación es la misma que la encontrada en Wikipedia. Porque a ^ 2 b ^ 2 = (ab) ^ 2 y (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
fuente
+1 ... excepto que creo que la fórmula final tiene un par fuera de lugar ... ¿no (b^4 sin(f))^2debería cambiarse a (b^4 sin(f)^2)?
Kirk Kuykendall
Realmente me alegro de que haya algunos expertos en este tema =).
Trevor Boyd Smith
¿Se puede publicar un archivo Geogebra (html) en este sitio? Tengo un radio de la vertical principal que podría demostrar visualmente lo que está sucediendo.
Puede exportar el original en formato .png, @Dan: use el cuadro de diálogo Archivo | Exportar. Recomiendo usar fuentes grandes (16 o 18 puntos parecen funcionar bien) y acercar la imagen lo más que pueda primero.
whuber
Asumo que la interactividad se perderá entonces. La demostración demuestra cómo la variación de los radios y la latitud de interés cambia las propiedades.
3

Es interesante descubrir que mi solución analfabeta matemática funcionó con 5 minutos de reflexión y codificación, ¿no debería considerarse el factor de aplanamiento en lugar de un modelo elíptico perfecto?

        double pRad = 6356.7523142;
        double EqRad = 6378.137;                      
        return pRad + (90 - Math.Abs(siteLatitude)) / 90 * (EqRad - pRad); 
Howard Grover
fuente
1
Donde pRad es el radio polar, y EqRad es el radio ecuatorial.
Stefan Steiger
Esta es la única respuesta que pude leer. Parece funcionar para mí
Sean Bradley
1
Veo que estás haciendo una interpolación lineal de radio entre el polo y el ecuador. Si bien no hay ninguna razón para creer que la interpolación lineal es precisa , usaré esto como "lo suficientemente bueno" para la Tierra, dado su leve factor de aplanamiento. Por cierto, creo que es un poco más fácil leer el equivalente: return E + (P - E) * Abs(Lat) / 90así que no es necesario tenerlo 90 - ...en la fórmula.
ToolmakerSteve
2

ingrese la descripción de la imagen aquí

Al menos esa es la fórmula que encontré en el Centro de Análisis y Evaluación de Datos de EE. UU. (DAAC) para el wiki del Programa de Modernización de la Computación de Alto Rendimiento (HPCMP) del Departamento de Defensa (DoD) . Dice que tomaron prestado mucho de la entrada de Wikipedia . Aún así, el hecho de que conservaron esa fórmula debería contar para algo.

RK
fuente
¿Puedes proporcionar un enlace al contenido?
Trevor Boyd Smith
donde φ es la latitud geodésica, y a (eje semi-mayor) yb (eje semi-menor) son, respectivamente, el radio ecuatorial y el radio 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