The n-vector page

Position calculations – simple and exact solutions by means of n-vector

n-vector is used to replace latitude and longitude, making the calculations simple and non-singular. Full accuracy is achieved for any global position (and for any distance).

This web page contains solutions to ten common geographical position calculations.

Background

Written by Kenneth Gade (PhD, Principal Scientist in the Navigation Group at FFI, Norwegian Defence Research Establishment).

For references, please cite:

Gade, K. (2010). A Non-singular Horizontal Position Representation, The Journal of Navigation, Volume 63, Issue 03, pp 395-417, July 2010.

Used worldwide for critical applications:

n-vector is used worldwide, e.g. for navigation, traffic control, and surveillance, of aircraft, cars, ships and submarines/submersibles. One example is Thales, who is world-leading at air traffic management (ATM). Their newest generation of ATM system is based on n-vector (replacing previous alternatives).

Useful pages

Solving ten common problems

We will now illustrate how ten common geographical position calculations can be solved using n-vector. The solutions are also available as Matlab code (and as C++, C#, Python, JavaScript, and other languages).

Color coding for figures and variables:

Input (red): The initial input that is given
Output (green): The solution to be found
Function (blue): A function that is available from the downloadable code.

Example 1: A and B to delta

Given two positions A and B. Find the exact vector from A to B in meters north, east and down, and find the direction (azimuth/bearing) to B, relative to north. Use WGS-84 ellipsoid.

Problem

See vector symbols explained for more details about the mathematical notation.

Given two positions, A and B as latitudes, longitudes and depths (relative to Earth, E):

–

l a t_{E A}

l o n g_{E A}

and

z_{E A}

(depth = –height)

–

l a t_{E B}

l o n g_{E B}

and

z_{E B}

Find the exact vector between the two positions, given in meters north, east, and down, and find the direction (azimuth) to B, relative to north.

Details:

Assume WGS-84 ellipsoid. The given depths are from the ellipsoid surface.
Use position A to define north, east, and down directions. (Due to the curvature of Earth and different directions to the North Pole, the north, east, and down directions will change (relative to Earth) for different places. Position A must be outside the poles for the north and east directions to be defined.

Solution

Note: The solutions are also available as downloadable code in several different programming languages.

1. First, the given latitudes and longitudes are converted to n-vectors:

\begin{array}{l} n_{E A}^{E} = l a t_l o n g 2 n_E (l a t_{E A}, l o n g_{E A}) \\ n_{E B}^{E} = l a t_l o n g 2 n_E (l a t_{E B}, l o n g_{E B}) \end{array}

2. When the positions are given as n-vectors (and depths), it is easy to find the delta vector decomposed in E (using the first function mentioned at the end of n-vector explained):

p_{A B}^{E} = n_E A_E_a n d_n_E B_E 2 p_A B_E (n_{E A}^{E}, n_{E B}^{E}, z_{E A}, z_{E B})

No ellipsoid is specified when calling the function, thus WGS-84 (default) is used.

3. We now have the delta vector from A to B, but (as seen from the superscript) the three coordinates of the vector are along the Earth coordinate frame E, while we need the coordinates to be north, east and down. To get this, we define a North-East-Down coordinate frame called N, and then we need the rotation matrix (direction cosine matrix)

R_{E N}^{}

to go between E and N (as explained here). In our math library, we have a simple function that calculates

R_{E N}^{}

from n-vector, and we use this function (using n-vector at A-position):

R_{E N}^{} = n_E 2 R_E N (n_{E A}^{E})

4. Now the delta vector is easily decomposed in N. When deciding if we should premultiply

p_{A B}^{E}

with

R_{E N}^{}

R_{N E}^{}

, we use the "closest-rule" saying that the subscript that is closest to the vector should be the frame where the vector is decomposed (also explained here). Since the vector is decomposed in E, we must use

R_{N E}^{}

(

R_{N E}^{}

is the transpose of

R_{E N}^{}

p_{A B}^{N} = R_{N E}^{} p_{A B}^{E}

5. The three components of

p_{A B}^{N}

are the north, east and down displacements from A to B in meters. The azimuth is simply found from element 1 and 2 of the vector (the north and east components):

a z i m u t h = atan2 (p_{A B}^{N} (2), p_{A B}^{N} (1))

Example 2: B and delta to C

Given the position of vehicle B and a bearing and distance to an object C. Find the exact position of C. Use WGS-72 ellipsoid.

Problem

See vector symbols explained for more details about the mathematical notation.

A radar or sonar attached to a vehicle B (Body coordinate frame) measures the distance and direction to an object C. We assume that the distance and two angles measured by the sensor (typically bearing and elevation relative to B) are already converted (by converting from spherical to Cartesian coordinates) to the vector

p_{B C}^{B}

(i.e. the vector from B to C, decomposed in B). The position of B is given as

n_{E B}^{E}

and

z_{E B}^{}

, and the orientation (attitude) of B is given as

R_{N B}^{}

(this rotation matrix can be found from roll/pitch/yaw by using zyx2R.m).

Find the exact position of object C as n-vector and depth (

n_{E C}^{E}

and

z_{E C}^{}

), assuming Earth ellipsoid with semi-major axis a and flattening f. For WGS-72, use a = 6 378 135 m and f = 1/298.26.

Solution

Note: The solutions are also available as downloadable code in several different programming languages.

1. The delta vector is given in B. It should be decomposed in E before using it, and thus we need

R_{E B}^{}

. This matrix is found from the matrices

R_{E N}^{}

and

R_{N B}^{}

, and we need to find

R_{E N}^{}

, as in Example 1:

R_{E N}^{} = n_E 2 R_E N (n_{E B}^{E})

2. Now, we can find

R_{E B}^{}

by using that the closest frames cancel when multiplying two rotation matrices (i.e. N is cancelled here):

R_{E B}^{} = R_{E N}^{} R_{N B}^{}

3. The delta vector is now decomposed in E (details about decomposing are found here):

p_{B C}^{E} = R_{E B}^{} p_{B C}^{B}

4. It is now easy to find the position of C using the second function mentioned here (with custom ellipsoid overriding the default WGS-84).

[n_{E C}^{E}, z_{E C}] = n_E A_E_a n d_p_A B_E 2 n_E B_E (n_{E B}^{E}, p_{B C}^{E}, z_{E B}, a, f)

Details: For simplicity, the orientation of the vehicle in this example was given relative to N. If the navigation system of the vehicle was using the non-singular coordinate frame L instead, the solution would be the same, but replacing the function n_E2R_EN.m with n_E_and_wa2R_EL.m. See Table 2 in Gade (2010) and the help text in n_E_and_wa2R_EL.m for details.

Example 3: ECEF-vector to geodetic latitude

Given an ECEF-vector of a position. Find geodetic latitude, longitude and height (using WGS-84 ellipsoid).

Problem

See vector symbols explained for more details about the mathematical notation.

Position B is given as an “ECEF-vector”

p_{E B}^{E}

(i.e. a vector from E, the center of the Earth, to B, decomposed in E).

Find the geodetic latitude, longitude and height,

l a t_{E B}

l o n g_{E B}

and

h_{E B}

, assuming WGS-84 ellipsoid.

Solution

Note: The solutions are also available as downloadable code in several different programming languages.

1. We have a function that converts p-vector to n-vector, see Appendix B in Gade (2010):

[n_{E B}^{E}, z_{E B}] = p_E B_E 2 n_E B_E (p_{E B}^{E})

2. Find latitude, longitude and height from

n_{E B}^{E}

and

z_{E B}

[l a t_{E B}, l o n g_{E B}] = n_E 2 l a t_l o n g (n_{E B}^{E})

h_{E B}

= –

z_{E B}

Example 4: Geodetic latitude to ECEF-vector

Given geodetic latitude, longitude and height. Find the ECEF-vector (using WGS-84 ellipsoid).

Problem

See vector symbols explained for more details about the mathematical notation.

Geodetic latitude, longitude and height are given for position B as

l a t_{E B}

l o n g_{E B}

and

h_{E B}

, find the ECEF-vector for this position,

p_{E B}^{E}

Solution

Note: The solutions are also available as downloadable code in several different programming languages.

1. First, the given latitude and longitude are converted to n-vector:

n_{E B}^{E} = l a t_l o n g 2 n_E (l a t_{E B}, l o n g_{E B})

2. We want

p_{E B}^{E}

from

n_{E B}^{E}

, and have functions that convert between these two. Now we need n_EB_E2p_EB_E.m. See also Appendix B in Gade (2010).

p_{E B}^{E} = n_E B_E 2 p_E B_E (n_{E B}^{E}, - h_{E B})

In the remaining examples, we assume that a spherical Earth model is sufficiently accurate, and then the position functions used above are not needed. The examples demonstrate the simplicity of these calculations when using n-vector instead of latitude and longitude. The n-vector solutions are also non-singular for all Earth-positions and work across the Date Line (180th meridian).

Example 5: Surface distance

Given position A and B. Find the surface distance (i.e. great circle distance) and the Euclidean distance.

Problem

See vector symbols explained for more details about the mathematical notation.

Given two positions A and B as

n_{E A}^{E}

n_{E B}^{E}

Find the surface distance

s_{A B}

(i.e. great circle distance). The heights of A and B are not relevant (i.e. if they do not have zero height, we seek the distance between the points that are at the surface of the Earth, directly above/below A and B). The Euclidean distance (chord length)

d_{A B}

should also be found (for

d_{A B}

, nonzero heights can be used). Use Earth radius

r_{E a r t h}

Solution

Note: The solutions are also available as downloadable code in several different programming languages.

The surface distance can be found in three different ways; by utilizing the dot product of the vectors, the cross product, or a combination:

s_{A B}^{} = a c o s (n_{E A}^{E} \cdot n_{E B}^{E}) \cdot r_{E a r t h} = a s i n (|n_{E A}^{E} \times n_{E B}^{E}|) \cdot r_{E a r t h} = a t a n 2 (|n_{E A}^{E} \times n_{E B}^{E}|, n_{E A}^{E} \cdot n_{E B}^{E}) \cdot r_{E a r t h}

The atan2 expression is numerically well-conditioned for all distances, see also equation (16) in Gade (2010).

The Euclidean distance is found by

d_{A B}^{} = |n_{E A}^{E} - n_{E B}^{E}| \cdot r_{E a r t h}

If nonzero heights for A and B are wanted for the Euclidean distance, each n-vector should be multiplied with

r_{E a r t h}

h_{i}

, where

h_{i}

is the corresponding height of each n-vector. If ellipsoidal Earth model is wanted, use n_EA_E_and_n_EB_E2p_AB_E.m from Example 1, and take the norm of the output vector.

Example 6: Interpolated position

Given the position of B at time

t_{0}

and

t_{1}

. Find an interpolated position at time

t_{i}

Problem

See vector symbols explained for more details about the mathematical notation.

Given the position of B at time

t_{0}

and

t_{1}

n_{E B}^{E} (t_{0})

and

n_{E B}^{E} (t_{1})

Find an interpolated position at time

t_{i}

n_{E B}^{E} (t_{i})

. All positions are given as n-vectors.

Solution

Note: The solutions are also available as downloadable code in several different programming languages.

Standard interpolation can be used directly with n-vector:

n_{E B}^{E} (t_{i}) = u n i t (n_{E B}^{E} (t_{0}) + (t_{i} - t_{0}) \frac{n_{E B}^{E} (t_{1}) - n_{E B}^{E} (t_{0})}{t_{1} - t_{0}})

Example 7: Mean position/center

Given three positions A, B, and C. Find the mean position (center/midpoint).

Problem 7

See vector symbols explained for more details about the mathematical notation.

Three positions A, B, and C are given as n-vectors

n_{E A}^{E}

n_{E B}^{E}

, and

n_{E C}^{E}

. Find the mean position, M, given as

n_{E M}^{E}

. Note that the calculation is independent of the heights/depths of the positions.

Solution 7

Note: The solutions are also available as downloadable code in several different programming languages.

The mean position is simply given by the mean n-vector:

n_{E M}^{E} = u n i t (n_{E A}^{E} + n_{E B}^{E} + n_{E C}^{E}) = \frac{n_{E A}^{E} + n_{E B}^{E} + n_{E C}^{E}}{|n_{E A}^{E} + n_{E B}^{E} + n_{E C}^{E}|}

Details about the definition of the horizontal geographical mean position

Several definitions of a geographical mean/center/midpoint are possible. The solution given here can be described as the “center of gravity” of the given positions, and can be compared to the centroid of a geometrical shape. Informally, the found mean position can be described by the following: Assume a sphere that is a model of the Earth, and let the sphere roll freely on a horizontal plane (in a uniform gravitational field). Place weights at the given positions of the sphere, and the weights will pull one side of the sphere down due to the gravitation (assuming the positions are not antipodal). When the sphere is in equilibrium, the “center of gravity”-position is the tangent point between the sphere and the plane.

Another possible way to define a midpoint is the position where the sum of surface distances (great circle distances) from the original positions is at a minimum. This midpoint is undefined when only two positions are given. If three positions are given in one dimension as 0, 0, and 3, this midpoint is at 0, while the arithmetic mean is 1. Iterations are probably needed to calculate this midpoint, and code for this is not included at this web site.

Example 8: A and azimuth/distance to B

Given position A and an azimuth/bearing and a (great circle) distance. Find the destination point B.

Problem 8

See vector symbols explained for more details about the mathematical notation.

Position A is given as n-vector

n_{E A}^{E}

. We also have an initial direction of travel given as an azimuth (bearing) relative to north (clockwise), and finally the distance to travel along a great circle is given as

s_{A B}

. Use Earth radius

r_{E a r t h}

Find the destination point B, given as

n_{E B}^{E}

In geodesy, this is known as "The first geodetic problem" or "The direct geodetic problem" for a sphere, and we see that this is similar to Example 2, but now the delta is given as an azimuth and a great circle distance. "The second/inverse geodetic problem" for a sphere is already solved in Examples 1 and 5.

Solution 8

Note: The solutions are also available as downloadable code in several different programming languages.

The azimuth (relative to north) is a singular quantity (undefined at the Poles), but from this angle we can find a (non-singular) quantity that is more convenient when working with vector algebra: a vector

d^{E}

that points in the initial direction. We find this from azimuth by first finding the north and east vectors at the start point, with unit lengths (see also equations (9) and (10) in Gade (2010), and the figure below):

\begin{array}{l} k_{e a s t}^{E} = u n i t ([\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}] \times n_{E A}^{E}) \\ k_{n o r t h}^{E} = n_{E A}^{E} \times k_{e a s t}^{E} \end{array}

Here we have assumed that our coordinate frame E has its z-axis along the rotational axis of the Earth, pointing towards the North Pole. Hence, this axis is given by

[\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}]

nvector globe example 8 — The north and east vectors are simple to find from the n-vector. The gray Earth-rotation-axis vector is also drawn at position A for convenience (vectors can be drawn at any position since they only have direction and magnitude).

The two vectors

k_{n o r t h}^{E}

and

k_{e a s t}^{E}

are horizontal, orthogonal, and span the tangent plane at the initial position. A unit vector

d^{E}

in the direction of the azimuth is now given by:

d^{E} = k_{n o r t h}^{E} \cdot \cos (a z i m u t h) + k_{e a s t}^{E} \cdot \sin (a z i m u t h)

With the initial direction given as

d^{E}

instead of azimuth, it is now quite simple to find

n_{E B}^{E}

. We know that

d^{E}

and

n_{E A}^{E}

are orthogonal, and they will span the plane where

n_{E B}^{E}

will lie. Thus, we can use sin and cos in the same manner as above, with the angle travelled given by

s_{A B}

r_{E a r t h}

n_{E B}^{E} = n_{E A}^{E} \cdot \cos (\frac{s_{A B}}{r_{E a r t h}}) + d^{E} \cdot \sin (\frac{s_{A B}}{r_{E a r t h}})

Example 9: Intersection of two paths / triangulation

Given path A going through

A_{1}

and

A_{2}

, and path B going through

B_{1}

and

B_{2}

. Find the intersection of the two paths.

Variant: Each of the two paths is specified by start point and azimuth instead of two points, which corresponds to triangulation.

Credits: Some readers of the paper Gade (2010) have contacted the author and reported new examples (other than those included in the paper) they have solved by means of n-vector. The following example has been reported twice, and is thus included here. Thanks to Even Børhaug at Kongsberg Maritime and Chris Veness at Movable Type Ltd, for this example. Børhaug reported that the solution he found that was based on latitude and longitude used multiple pages of code, while Veness says that such a solution is "horribly complex". Using n-vector, the intersection is very simple to find.

Problem 9

See vector symbols explained for more details about the mathematical notation.

Define a path from two given positions (at the surface of a spherical Earth), as the great circle that goes through the two points (assuming that the two positions are not antipodal).

Path A is given by

n_{E A,1}^{E}

and

n_{E A,2}^{E}

, while path B is given by

n_{E B,1}^{E}

and

n_{E B,2}^{E}

Find the position C where the two paths intersect, as

n_{E C}^{E}

Solution 9

Note: The solutions are also available as downloadable code in several different programming languages.

A convenient way to represent a great circle is by its normal vector (i.e. the normal vector to the plane containing the great circle). This normal vector is simply found by taking the cross product of the two n-vectors defining the great circle (path). Having the normal vectors to both paths, the intersection is now simply found by taking the cross product of the two normal vectors:

n_{E C}^{E} = \pm u n i t ((n_{E A,1}^{E} \times n_{E A,2}^{E}) \times (n_{E B,1}^{E} \times n_{E B,2}^{E}))

Note that there will be two places where the great circles intersect, and thus two solutions are found. Selecting the solution that is closest to e.g.

n_{E A,1}^{E}

can be achieved by selecting the solution that has a positive dot product with

n_{E A,1}^{E}

(or the mean position from Example 7 could be used instead of

n_{E A,1}^{E}

). In program code, this can be implemented by selecting any of the two solutions above (+ or –), called C', and calculate:

n_{E C}^{E} = s i g n (n_{E C'}^{E} \cdot n_{E A,1}^{E}) n_{E C'}^{E}

Variant / triangulation: If each of the two paths is specified by start point and azimuth instead of two points, the normal vectors to the great circles can be found by the following:

1. Find the direction vector for each path in the same manner as

d^{E}

in Example 8.

2. Take the cross product of the direction vector and the n-vector for the initial position.

Example 10: Cross track distance (cross track error)

Given path A going through

A_{1}

and

A_{2}

, and a point B. Find the cross track distance/cross track error between B and the path.

Variant: The path is specified by start point and azimuth instead of two points.

(Along track distance (and the corresponding point on the path) can also be found quite easily)

Credits: Thanks to Chris Veness at Movable Type Ltd who first reported this example.

Problem 10

See vector symbols explained for more details about the mathematical notation.

Path A is given by the two n-vectors

n_{E A,1}^{E}

and

n_{E A,2}^{E}

(as in the previous example), and a position B is given by

n_{E B}^{E}

Find the cross track distance

s_{x t}

between the path A (i.e. the great circle through

n_{E A,1}^{E}

and

n_{E A,2}^{E}

) and the position B (i.e. the shortest distance at the surface, between the great circle and B). Also, find the Euclidean distance

d_{x t}

between B and the plane defined by the great circle. Use Earth radius

r_{E a r t h}

Solution 10

Note: The solutions are also available as downloadable code in several different programming languages.

First, find the normal

c^{E}

to the great circle, with direction given by the right hand rule and the direction of travel:

c^{E} = u n i t (n_{E A,1}^{E} \times n_{E A,2}^{E})

We now see that the cross track distance to B is the great circle distance from

c^{E}

minus 90° (π/2 radians):

s_{x t}^{} = (a c o s (c^{E} \cdot n_{E B}^{E}) - \frac{π}{2}) \cdot r_{E a r t h} = - a s i n (c^{E} \cdot n_{E B}^{E}) \cdot r_{E a r t h}

Finding the Euclidean distance is even simpler, since it is the projection of

n_{E B}^{E}

onto

c^{E}

, thus simply the dot product:

d_{x t}^{} = - c^{E} \cdot n_{E B}^{E} \cdot r_{E a r t h}

For both

s_{x t}

and

d_{x t}

, positive answers means that B is to the right of the track.

Variant: If the path is specified by start point and azimuth instead of two points, the normal vector to the great circle can be found by the following:

1. Find the direction vector for the path in the same manner as

d^{E}

in Example 8.

2. Take the cross product of the direction vector and the n-vector for the initial position.

Along track distance to the closest point on path A can be found by using the four vectors in this example, and then calculating the intersection point as in Example 9 (when the n-vector for the intersection point is found, the along track distance from

A_{1}

is found as in Example 5). Thanks to Enrico Spinielli at EUROCONTROL for reporting this.

About the solutions

The four first examples are simply solved by using the functions provided for download. The functions are exact, and can use ellipsoidal or spherical Earth model. They are without iterations (i.e. on closed form), and WGS-84 reference ellipsoid is used as default. Custom ellipsoids/spheres can be specified. The functions are based on n-vector instead of latitude and longitude.

The last six examples are solved for spherical Earth, and then no functions are needed. Using n-vector instead of latitude and longitude, these examples are quite straightforward to solve by simple vector algebra.

The solutions to all 10 examples are non-singular, and they work equally well for all global positions (i.e. they work at the Poles, close to the Poles, and across the Date Line). The solutions also have full numerical accuracy for both short and long distances (they are numerically well-conditioned).

Useful pages

The n-vector page

Solving ten common problems

Example 1: A and B to delta

Problem

Solution

Example 2: B and delta to C

Problem

Solution

Example 3: ECEF-vector to geodetic latitude

Problem

Solution

Example 4: Geodetic latitude to ECEF-vector

Problem

Solution

Example 5: Surface distance

Problem

Solution

Example 6: Interpolated position

Problem

Solution

Example 7: Mean position/center

Problem 7

Solution 7

Details about the definition of the horizontal geographical mean position

Example 8: A and azimuth/distance to B

Solution 8

Example 9: Intersection of two paths / triangulation

Problem 9

Solution 9

Example 10: Cross track distance (cross track error)

Problem 10

Solution 10

About the solutions

Useful pages

n-vector downloads

Vector symbols explained

n-vector explained