4
$\begingroup$

Do these conversion formulas look correct? I'm writing my own local navigation script in Python, and I want to ensure that my understanding of the underlying math is solid.

Reference location
Conversion from GPS to ECEF:

$$ \begin{pmatrix} X_{r}\\ Y_{r}\\ Z_{r} \end{pmatrix} = \begin{pmatrix} \left(N_{r}+h_{r}\right)\,\cos\left(\phi_{r}\right)\,\cos\left(\lambda_{r}\right)\\ \left(N_{r}+h_{r}\right)\,\cos\left(\phi_{r}\right)\,\sin\left(\lambda_{r}\right)\\ \left(N_{r}\left(1-e^{2}\right)+h_{r}\right)\,\sin\left(\phi_{r}\right) \end{pmatrix} $$

Aircraft location
Conversion from GPS to ECEF:

$$ \begin{pmatrix} X_{p}\\ Y_{p}\\ Z_{p} \end{pmatrix} = \begin{pmatrix} \left(N_{p}+h_{p}\right)\,\cos\left(\phi_{p}\right)\,cos\left(\lambda_{p}\right)\\ \left(N_{p}+h_{p}\right)\,\cos\left(\phi_{p}\right)\,\sin\left(\lambda_{p}\right)\\ \left(N_{p}\left(1-e^{2}\right)+h_{p}\right)\,\sin\left(\phi_{p}\right) \end{pmatrix} $$

Aircraft-reference displacement vector:

$$ \Delta_{r_{ECEF}} = \begin{pmatrix} \Delta X\\ \Delta Y\\ \Delta Z \end{pmatrix} = \begin{pmatrix} X_{p}-X_{r}\\ Y_{p}-Y_{r}\\ Z_{p}-Z_{r} \end{pmatrix} $$

Direction-Cosine-Matrix (DCM):

$$ \mathbf{R}\left(\phi_{r},\,\lambda_{r}\right) = \begin{pmatrix} -\sin\left(\phi_{r}\right)\,\cos\left(\lambda_{r}\right) & -\sin\left(\phi_{r}\right)\,\sin\left(\lambda_{r}\right) & \cos\left(\phi_{r}\right)\\ \cos\left(\lambda_{r}\right) & -\sin\left(\lambda_{r}\right) & 0\\ -\cos\left(\phi_{r}\right)\,\cos\left(\lambda_{r}\right) & -\cos\left(\phi_{r}\right)\,\sin\left(\lambda_{r}\right) & -\sin\left(\phi_{r}\right) \end{pmatrix} $$

North-East-Down (NED) conversion:

$$ \begin{pmatrix} N\\ E\\ D \end{pmatrix} = \begin{pmatrix} -\sin\left(\phi_{r}\right)\,\cos\left(\lambda_{r}\right) & -\sin\left(\phi_{r}\right)\,\sin\left(\lambda_{r}\right) & \cos\left(\phi_{r}\right)\\ \cos\left(\lambda_{r}\right) & -\sin\left(\lambda_{r}\right) & 0\\ -\cos\left(\phi_{r}\right)\,\cos\left(\lambda_{r}\right) & -\cos\left(\phi_{r}\right)\,\sin\left(\lambda_{r}\right) & -\sin\left(\phi_{r}\right) \end{pmatrix} \begin{pmatrix} \Delta X\\ \Delta Y\\ \Delta Z \end{pmatrix} $$

$\endgroup$
3
  • 3
    $\begingroup$ Don't think this is an Aviation topic. Sure, planes use GPS, but so do cars, and hikers, and cell phones, and plenty of other things. Maybe more on-topic at edit: GIS. $\endgroup$
    – Ralph J
    Commented Nov 7, 2023 at 17:02
  • 1
    $\begingroup$ @RalphJ I have a different perspective on that; these coordinate transformations are typical for aviation system developers, whether you are working on FMS navigation implementations, radar data fusion and tracking, flight simulation (both for ATC and for flight training) and probably a few more cases. I agree it may have no direct relation to piloting an aircraft, but our scope is wider than that. Certainly it does not belong at Electrical Engineering, perhaps at GIS it would get an answer. For now I would leave it here, and if nobody answers it, I may have a dig at it myself. $\endgroup$
    – DeltaLima
    Commented Nov 7, 2023 at 18:47
  • $\begingroup$ @DeltaLima Yeah, GIS not EE, that seems like definitely the better fit. But "how do I implement these formula in Python" seems pretty far outside of Aviation. The amount of coding that gets posted here in answers is almost zero. That said, I'll leave it be & let the community weigh in to close or leave open. $\endgroup$
    – Ralph J
    Commented Nov 7, 2023 at 20:38

1 Answer 1

5
$\begingroup$

If I understand it correctly, you would like to calculate the position of the aircraft in a North East Down (NED) reference frame. The origin of that reference frame is at the reference location.

The standard approach to do this uses 3 steps:

  1. convert both positions into the same carthesian frame (usually Earth Centered, Earth Fixed, ECEF)
  2. move the origin of the first frame to the reference point
  3. transform to the coordinates by rotating the reference frame to the local NED orientation.

You are performing these steps in your formula, let's take a look at them one-by-one.


Your first step converts the positions of the aircraft and the reference location from geodetic latitude, longitude, height above ellipsoid ($\phi ,\lambda, h$) into the ECEF frame.

Your conversion seems correct, but for completeness sake, I will add here that you have to calculate $N_r$ and $N_p$ based on their geodetic latitude $\phi$:

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

With:

$\begin{eqnarray} a &=& \text{length of Earth's semi-major axis}\\ b &=& \text{length of Earth's semi-minor axis}\\ e &=& \sqrt{1-{\frac {b^{2}}{a^{2}}}} \text{the first numerical eccentricity of the ellipsoid}. \\ \end {eqnarray} $

The second step is a simple subtraction, which seems to be correct in your approach.

The third step is the rotation of the frame to the NED attitude, which you do using the direction cosine matrix. It seems to me that there is a mistake there.

I used the following verification test:

  • Set the reference point at $(0,0,0)$, meaning it is the intersection of the prime meridian and the equator, at the ellipsoid surface.
  • In ECEF that would be $(a,0,0)$.
  • Suppose our aircraft is 100 m east of that location, in ECEF: $(a,100,0)$,
  • the difference vector being $(0,100,0)$.
  • calculating the east coordinate, according to your DCM would be $0\cdot\cos(\lambda_r) -100\cdot\sin(\lambda_r) = 0$, but the expected result would be 100.

Therefore, my conclusion is that your DCM is wrong.


I have this DCM matrix for conversion to the ENU (east north up) frame:

$ \begin{pmatrix} E \\ N \\ U \end{pmatrix} = \begin{bmatrix} -\sin(\lambda_r) & \cos(\lambda_r) & 0 \\ -\sin(\phi_r)\cos(\lambda_r)& -\sin(\phi_{r})\sin(\lambda_{r}) & \cos(\phi_r) \\ \cos(\phi_r)\cos(\lambda_{r}) & \cos(\phi_{r})\sin(\lambda_{r}) & \sin(\phi_r) \end{bmatrix} \begin{bmatrix} X_{p}-X_{r}\\ Y_{p}-Y_{r}\\ Z_{p}-Z_{r} \end{bmatrix} $

Rearanging the first two rows and negating the third, would yield the NED DCM matrix:

$ \begin{pmatrix} N \\ E \\ D \end{pmatrix} = \begin{bmatrix} -\sin(\phi_r)\cos(\lambda_r)& -\sin(\phi_{r})\sin(\lambda_{r}) & \cos(\phi_r) \\ -\sin(\lambda_r) & \cos(\lambda_r) & 0 \\ -\cos(\phi_r)\cos(\lambda_{r}) & -\cos(\phi_{r})\sin(\lambda_{r}) & -\sin(\phi_r) \end{bmatrix} \begin{bmatrix} X_{p}-X_{r}\\ Y_{p}-Y_{r}\\ Z_{p}-Z_{r} \end{bmatrix} $

Performing the same quick verification test yields $-0\cdot\sin(\lambda_r) + 100\cdot\cos(\lambda_r) = 100$, and that is the expected result.

This is only proof that the result is not wrong for this simple input; I recommend doing further verification to ensure general correctness of both the formula and your implementation.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .