¿Cómo se calcula la latitud y la longitud de un punto de la tierra utilizando el acimut y la elevación vistos desde un satélite geoestacionario?

nota: creo que las coordenadas terrestres deben escribirse con mayúsculas$(X, Y, Z)$pero las citas de la otra respuesta usan minúsculas$(x, y, z)$. Teniendo en cuenta que el término radical apenas encaja con las letras pequeñas, he mantenido las minúsculas en la mayoría de los lugares. Para esta respuesta, las mayúsculas y minúsculas representan las mismas coordenadas.

tl; dr: conecte la posición de su satélite ($x, y, z$)

$$\mathbf{x_0} = x\mathbf{\hat{x}} + y\mathbf{\hat{y}} + z\mathbf{\hat{z}}$$ y la normal del vector ($u, v, w$) $$\mathbf{\hat{n}} = u\mathbf{\hat{x}} + v\mathbf{\hat{y}} + w\mathbf{\hat{z}}$$ Llegar$t$la longitud del vector desde el satélite hasta el elipsoide, luego use $$\mathbf{x_{intersect}} = \mathbf{x_0} + t\mathbf{\hat{n}}$$ para obtener las coordenadas de la intersección con el elipsoide de referencia.

Luego use la solución analítica en la segunda sección a continuación, o una de las soluciones iterativas a las que se hace referencia allí para convertirlas a latitud y longitud.

---

Matemáticas de intersección para obtener (x, y, z)

De esta práctica respuesta en GIS SE:

El elipsoide de referencia WGS84 es un elipsoide biaxial (y achatado). Es más corto en los polos que el ecuador, y el ecuador es un círculo.

La ecuación para ello es:

$$\frac{x^2}{a^2} + \frac{y^2}{a^2} + \frac{z^2}{b^2} - 1 = 0$$

dónde$x, y, z$es un punto en la superficie del elipsoide, y

dónde$a$es el semieje mayor (6.378.137 metros) y$b$es el semieje menor del elipsoide WGS84 (6.356.752,3142 metros).

Continuaré usando las matemáticas de esta respuesta y formatearé en MathJax y Python posterior:

La posición de su satélite es$x,y, z$y la dirección normal es$u, v, w$. La longitud del vector desde el satélite hasta la primera y más cercana intersección (asumiendo, por supuesto, que su satélite está fuera de la Tierra ):

t = -(1/(b^2 (u^2 + v^2) +  a^2 w^2)) * (b^2 (u x + v y) + a^2 w z + 1/2 Sqrt[
     4 (b^2 (u x + v y) + a^2 w z)^2 - 
     4 (b^2 (u^2 + v^2) + a^2 w^2) (b^2 (-a^2 + x^2 + y^2) + a^2 z^2)])

En MathJax:

$$A = -\frac{1}{b^2(u^2+v^2) + a^2w^2}$$

$$B = b^2(ux + vy) + a^2wz$$

$$C = \sqrt{(b^2(ux + vy) + a^2wz)^2 - (b^2(u^2 + v^2) + a^2w^2) (b^2(-a^2 + x^2 + y^2) + a^2z^2)}$$

$$t = A(B+C)$$

que es más fácil de leer que la forma en que MathJax lo muestra todo en una sola línea:

$$t = -\frac{1}{b^2(u^2+v^2) + a^2w^2} \left( b^2(ux + vy) + a^2wz + \sqrt{(b^2(ux + vy) + a^2wz)^2 - (b^2(u^2 + v^2) + a^2w^2) (b^2(-a^2 + x^2 + y^2) + a^2z^2)} \right)$$

Implementado en Python:

import numpy as np
a, b = 0.9, 1.1
asq, bsq = a**2, b**2
x, y, z = 5.0, 0.0, 0.0
u, v, w = -np.sqrt(1 - 0.1**2 - 0.1**2), 0.1, 0.1
xsq, ysq, zsq = x**2, y**2, z**2
usq, vsq, wsq = u**2, v**2, w**2

A = -(1/(bsq*(usq + vsq) +  asq*wsq))
B = bsq*(u*x + v*y) + asq*w*z
C = 0.5*np.sqrt(4*(bsq*(u*x + v*y) + asq*w*z)**2 -
                      4*(bsq*(usq + vsq) + asq*wsq) *
                      (bsq*(-asq + xsq + ysq) + asq*zsq))
t = A * (B + C)
print "t: ", t
xyz, uvw   = np.array([x, y, z]), np.array([u, v, w])
xyzi       = xyz + t*uvw
xi, yi, zi = xyzi
print "point of intersection: ", xyzi
print "check, is it zero within roundoff? ", xi**2/asq + yi**2/asq + zi**2/bsq - 1

cede:

t:  4.33961998769
point of intersection:  [ 0.70399539  0.433962    0.433962  ]
check, is it zero?  -8.54871728961e-15

De hecho, la solución cae en el elipsoide.

Convertir (x, y, z) a latitud y longitud

Necesitamos invertir las ecuaciones (que se muestran en esta respuesta ) para obtener la latitud y la longitud en el elipsoide WSG84 correspondiente a estas coordenadas.

Según [Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates][1] de Wikipedia

Las coordenadas cartesianas 3D$X, Y, Z$en coordenadas centradas en la Tierra y fijas en la Tierra asumiendo una forma elipsoidal viene dada por:

$$X = \left(N(\phi) + h \right) \cos\phi \cos\lambda$$

$$Y = \left(N(\phi) + h \right) \cos\phi \sin\lambda$$

$$Z = \left(\frac{b^2}{a^2} N(\phi) + h \right) \sin\phi$$

dónde$\phi, \lambda, h$son latitud, longitud y altitud, y$a, b$son los radios ecuatorial y polar del elipsoide utilizado, y

$$N(\phi) = \frac{a^2}{\sqrt{a^2\cos^2\phi + b^2 \sin^2\phi}}.$$

La inversión se trata en Transformaciones de Datum de Posiciones GPS; Nota de aplicación que se encuentra en esta pregunta GIS , así como en este enlace en los comentarios.

De la nota de aplicación mostraré la solución analítica .

Definiciones:

$$a = 6378137$$ $$b = a(1-f) = 6356752.31424518$$ $$f = \frac{1}{298.257223563}$$ $$e = \sqrt{\frac{a^2-b^2}{a^2}}$$ $$e' = \sqrt{\frac{a^2-b^2}{b^2}}.$$

Resolver para la longitud es trivial. La nota de aplicación da:

$$\lambda = \arctan\left(\frac{y}{x}\right),$$

sin embargo, en la práctica, no puede usar arctan a ciegas porque no puede distinguir los cuatro cuadrantes. Así que en su lugar usa

$$\lambda = \arctan2(y, x).$$

Resolviendo analíticamente para la latitud, la nota de aplicación da

$$\phi= \arctan\left( \frac{z + e'^2b \sin^3(\theta)}{p-e^2a\cos^3(\theta)} \right)$$

dónde

$$p = \sqrt{x^2+y^2}$$ $$\theta = \arctan\left(\frac{za}{pb}\right).$$

Nuevamente, probablemente deberías usar$arctan2()$en lugar de cuidar los cuadrantes.

También puede encontrar soluciones iterativas en esta referencia y en la otra (ya mencionada anteriormente). Las soluciones iterativas probablemente fueron significativamente más rápidas en el pasado (y podrían serlo también ahora), lo que sería importante si estuviera haciendo un trazado de rayos para millones de puntos, por ejemplo, generando una imagen de la superficie de la Tierra desde el punto de vista de los satélites con una precisa representación.