Movable Type Home Page

Movable Type Scripts


Vector-based geodesy

Vector-based methods provide an alternative approach to latitude/longitude geodesy calculations from the more common spherical trigonometry methods.

Many calculations, especially more complex ones, are far simpler using vector-based calculations. Also, special cases such as the poles and the date line create fewer problems.

These notes set out the maths and code for such calculations. This page is more condensed than the main latitude/longitude page, which includes far more explanations and context.

Points

great circle vector

A point on the earth’s surface can be represented by an ‘n-vector’1, normal to the surface of the earth.

On a spherical model earth, this will be a vector originating from the centre of the earth. On an ellipsoidal model earth, it will intersect the equatorial plane at an angle equal to the geodetic latitude.

The formulae for converting between latitude/longitude points and n-vectors apply equally to spherical and ellipsoidal model earths.

From latitude/longitude to n-vector

The n-vector for a point φ,λ on the earth’s surface, where latitude = φ and longitude = λ, is defined as

  cosφ·cosλ  
v{x,y,z} =   cosφ·sinλ  
  sinφ  

From n-vector to latitude/longitude

For a given n-vector, the latitude φ and longitude λ of the point it represents on the surface is defined as

  • φ = atan2(vz, √vx² + vy²)
  • λ = atan2(vy, vx)

Note that these formulations take x pointing to 0°N,0°E, y pointing to 0°N,90°E, and z pointing to 90°N (making a right-handed vector coordinate system). See toNvector() and toLatLon() in the source code below.

I generally use Greek letters (eg φ/λ) for angles in radians, English (eg lat/lon) for angles in degrees (having found mixing degrees & radians is often the easiest way to confusing bugs!) – if you encounter any problems, ensure your <head> includes <meta charset="utf-8">.

Spherical model earth

The remainder of this page relates to a spherical model earth. I will be completing some functions relating to an ellipsoidal model earth in the near future.

great circle vector

Great circles

A great circle is a plane through the earth’s centre, and can be represented by a vector normal to the plane of the great circle. The direction of the vector gives the direction of travel around the great circle.

A great circle can be defined by a pair of points a, b on the surface. In this case, the great circle is the (normalised) cross product of the n-vectors representing the points, ĉ = a × b.

A great circle can also be defined by a start point and bearing; a vector defining a great circle defined by a start point φ,λ and a bearing θ is given by1

  sinλ·cosθ − sinφ·cosλ·sinθ  
c{x,y,z} =   −cosλ·cosθ − sinφ·sinλ·sinθ  
  cosφ·sinθ  

(the bearing θ being measured geodetically clockwise from north, not mathematically anti-clockwise from the x axis).

See greatCircle() in the source code below.

Distance between two points

On a spherical earth model, the great circle distance between two points (the ‘inverse geodetic problem’) is the angle between them (in radians) multiplied by the radius of the earth.

For two points a and b, the angle between them is the the arccos of the dot product, α = acos(a·b), or the arcsin of the magnitude of the cross product, α = asin(|a×b|). Better conditioned for all angles, however, is the arctan version1

where R is the radius of the earth (normally taken as 6371 km).

See distanceTo() in the source code below.

initial bearing

Initial bearing

To get the initial bearing (sometimes called ‘forward azimuth’) between two points a and b, if we take N as a vector representing the North Pole, then we can take two great circles, one passing through a and b, the other passing through a and N. The compass bearing from the first point to the second is the (signed) angle between the vectors perpendicular to these two planes (that is, the vectors representing those two great circles).1 To get angles greater than ±π/2 (±90°), the vector a can be used as a reference to determine the direction of c1×c2 hence the sign of its length.

For final bearing , simply take the initial bearing from the end point to the start point (‘reverse azimuth’) and reverse it (using θ = (θ+180) % 360).

See bearingTo() in the source code below.

Midpoint between two points

The midpoint between two points is represented by the halfway vector between them:

This can be easily extended to an intermediate point at a fraction f along a chord between the points by multiplying the vectors by the fraction:

(or equivalently = a + f·(ba)).

See midpointTo(), intermediatePointOnChordTo() in the source code below.

initial bearing

Destination point given start point, distance & bearing

The destination point (the ‘direct geodetic problem’) can be found by adding a direction vector to the start point n-vector.

Taking a as the n-vector representing the start point, δ the angular distance travelled, and θ the bearing (from north), the destination point n-vector b can be found from

The direction vector is also the cross product of the great circle and the start point, c×a.

See destinationPoint() in the source code below.

Intersection of two paths

This is horribly complex (at least to me) using spherical trigonometry, but quite straightforward using vectors.

Whether the paths are defined by pairs of points, or by start-point & bearing, the great circles for each of the paths can be obtained as explained above.

The cross product of the vectors defining the great circles will then be the n-vector defining the intersection point. The only complex part is determining which of the candidate (antipodal) intersection points is correct.

either: ni = c1 × c2
or: ni = c2 × c1

Using the bearing c × n, if c × ni is positive, the bearing is pointing to that intersection point; otherwise, use the other point.

See intersection() in the source code below.

Cross-track distance

The cross-track distance is the distance of a point from a great-circle path – the deviation from the path one should be following.

Since a great circle is represented by a vector perpendicular to the plane, if α is the angle between the vector representing the point and vector representing the great circle, then the distance from the point to the great circle is d = R · (90° − α).

The great circle vector can be obtained from a start-point and bearing, or from a pair of points.

See crossTrackDistanceTo() in the source code below.

triangulation

Triangulation

The location of a point can be determined from bearings from two known points.

The method is related to destination point given start point distance & bearing, except that two direction vectors are derived, from those a pair of great circle vectors, and then the intersection point is the cross product of the two great circles:

Trilateration

The location of a point can be determined if the distances to three other points are known; the derivation is given in Wikipedia.

Given three (geocentric) n-vectors p1, p2, p3, to the three points (using toNvector()), and (angular) distances r1, r2, r3 from the point of unknown location to those points, then

z will have 0, 1, or 2 solutions depending whether the spheres defined by r1, r2, r3 intersect. If we assume the points are on the plane of the sphere and ignore z, we will get a single solution even if the distances are slightly incorrect.

JavaScript:

    // the following uses x,y coordinate system with origin at P1, x axis P1->P2
    const eX = P2.minus(P1).unit();                        // unit vector in x direction P1->P2
    const i = eX.dot(P3.minus(P1));                        // signed magnitude of x component of P1->P3
    const eY = P3.minus(P1).minus(eX.times(i)).unit();     // unit vector in y direction
    const d = P2.minus(P1).length();                       // distance P1->P2
    const j = eY.dot(P3.minus(P1));                        // signed magnitude of y component of P1->P3
    const x = (r1*r1 - r2*r2 + d*d) / (2*d);               // x component of P1 -> intersection
    const y = (r1*r1 - r3*r3 + i*i + j*j) / (2*j) - x*i/j; // y component of P1 -> intersection
    const eZ = eX.cross(eY);                               // unit vector perpendicular to plane
    const z = Math.sqrt(r1*r1 - x*x - y*y);                // z will be NaN for no intersections

    const Pi = P1.plus(eX.times(x)).plus(eY.times(y)); // note don't use z component; assume points at same height
    

Area of a spherical polygon

The area of a spherical polygon can be found from the ‘spherical excess’ using Girard’s theorem:

where A is the area (in units of ), E is the spherical excess, n is the number of sides of the polygon, θᵢ is the n-th interior angle, and R is the earth’s radius.

Centre of a spherical polygon

As a software developer, I love it when maths that is rather beyond my comprehension results in a a solution which is really simple to code.

Apparently, it ‘follows from a non-obvious application of Stokes’ theorem* that the centre of a spherical polygon is

(actually, the ‘centre’ being the projection of the moment of the 3d centroid onto the unit sphere, but that amounts to the centre of the polygon).

JavaScript – given an array of LatLon vertices defining the spherical polygon:

        let centreV = new NvectorSpherical(0, 0, 0);
        for (let vertex=0; vertex<polygon.length; vertex++) {
            const a = polygon[vertex].toNvector();                      // current vertex
            const b = polygon[(vertex+1) % polygon.length].toNvector(); // next vertex
            const v = a.cross(b).unit().times(a.angleTo(b)/2);          // a×b / |a×b| ⋅ θab/2
            centreV = centreV.plus(v);
        }
    

(this could be done in a functional approach using Array.reduce(), but in JavaScript it seems clearer procedurally).

As a ‘polygon' on a sphere always has two interpretations (depending in this case on the direction of travel around the polygon), the simplicity of this formulation is compromised by having to select the appropriate polygon, to avoid getting a ‘centre’ on the opposite side of the globe. This could be rigorously achieved by taking the polygon with the smaller area; however, this involves a lot of work – a simpler solution is to reverse the centre vector if it is pointing in the opposite direction from the 1st vertex (this might give the ‘wrong’ solution for pathelogical semi-global cases, but will be fine for normal cases. See GitHub source code for a complete solution.

Live examples

... to verify the code is correct!

Enter the co-ordinates into the text boxes to try out the calculations. A variety of formats are accepted, principally:

Pair of points
Point 1: ,
Point 2: ,
Distance: km
Initial bearing:
Final bearing:
Midpoint:
Destination point along great-circle given distance and bearing from start point
Start point: ,
Bearing:
Distance: km
Destination point:
Final bearing:
Intersection of two great-circle paths
Point 1: , Brng 1:
Point 2: , Brng 2:
Intersection point:
Cross-track distance
Start point: ,
Bearing:
Current point: ,
Cross-track distance: km

Further functions

Geographic mean

The geographic mean of a number of latitude/longitude points (which would be complex at best using spherical trigonometry) can be obtained by summing the n-vectors (and normalizing the result).

Best-fit path through points

Since the geographic mean is so simple, I did wonder whether a best-fit RMS regression line through a set of points would be approachable; however, it is highly complex in 3 dimensions, even using vectors. The solution is apparently an eigenvector problem using singular value decomposition (SVD). 1, 2 I’ll leave that to the experts.

Further reading

More information is available in Kenneth Gade’s paper A Non-singular Horizontal Position Representation, The Journal of Navigation, volume 63, issue 03, The Royal Institute of Navigation, 2010; with further explanations on NavLab’s n-vector page, and MATLAB code on the Forsvarets forskningsinstitutt n-vector GitHub repository.


See below for the JavaScript source code, also available on GitHub. Full documentation is available, as well as a test suite.

Note I use Greek letters in variables representing maths symbols conventionally presented as Greek letters: I value the great benefit in legibility over the minor inconvenience in typing.

With its untyped C-style syntax, JavaScript reads remarkably close to pseudo-code: exposing the algorithms with a minimum of syntactic distractions. These functions should be simple to translate into other languages if required, though can also be used as-is in browsers and Node.js.

For convenience & clarity, I have extended the base JavaScript Number object with toRadians() and toDegrees() methods: I don’t see great likelihood of conflicts, as these are ubiquitous operations.

February 2019: I have refactored the library to use ES modules, as well as extending it in scope; see the GitHub README and CHANGELOG for details.

Other languages: I cannot support translations into other languages, but if you have ported the code to another language, I am happy to provide links here.

OSI MIT License I offer these scripts for free use and adaptation to balance my debt to the open-source info-verse. You are welcome to re-use these scripts [under an MIT licence, without any warranty express or implied] provided solely that you retain my copyright notice and a link to this page.

If you have any queries or find any problems, contact me at ku.oc.epyt-elbavom@oeg-stpircs.

© 2011-2022 Chris Veness