Vectorbased 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 vectorbased 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.
A point on the earth’s surface can be represented by an ‘nvector’^{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 nvectors apply equally to spherical and ellipsoidal model earths.
From latitude/longitude to nvectorThe nvector for a point φ,λ on the earth’s surface, where latitude = φ and longitude = λ, is defined as

From nvector to latitude/longitudeFor a given nvector, the latitude φ and longitude λ of the point it represents on the surface is defined as

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 righthanded 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="utf8">
.
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.
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 nvectors 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 by^{1}
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 anticlockwise from the x axis).
See greatCircle()
in the source code below.
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 version^{1}
where R is the radius of the earth (normally taken as 6371 km).
See distanceTo()
in the source code below.
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 c_{1}×c_{2} 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.
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 m̂ = a + f·(b−a)).
See midpointTo()
, intermediatePointOnChordTo()
in the source code below.
The destination point (the ‘direct geodetic problem’) can be found by adding a direction vector to the start point nvector.
Taking a as the nvector representing the start point, δ the angular distance travelled, and θ the bearing (from north), the destination point nvector 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.
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 startpoint & 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 nvector defining the intersection point. The only complex part is determining which of the candidate (antipodal) intersection points is correct.
either:  n_{i} = c_{1} × c_{2} 
or:  n_{i} = c_{2} × c_{1} 
Using the bearing c × n, if c × n ⋅ i is positive, the bearing is pointing to that intersection point; otherwise, use the other point.
See intersection()
in the source code below.
The crosstrack distance is the distance of a point from a greatcircle 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 startpoint and bearing, or from a pair of points.
See crossTrackDistanceTo()
in the source code below.
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:
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) nvectors p_{1}, p_{2}, p_{3},
to the three points (using toNvector()
), and (angular) distances r_{1},
r_{2}, r_{3} from the point of unknown location to those points, then
z will have 0, 1, or 2 solutions depending whether the spheres defined by r_{1}, r_{2}, r_{3} 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
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 R²), E is the spherical excess, n is the number of sides of the polygon, θᵢ is the nth interior angle, and R is the earth’s radius.
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 nonobvious 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 semiglobal cases, but will be fine for normal cases. See GitHub source code for a complete solution.
... to verify the code is correct!
Enter the coordinates into the text boxes to try out the calculations. A variety of formats are accepted, principally:
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 nvectors (and normalizing the result).
Since the geographic mean is so simple, I did wonder whether a bestfit 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.
More information is available in Kenneth Gade’s paper A Nonsingular Horizontal Position Representation, The Journal of Navigation, volume 63, issue 03, The Royal Institute of Navigation, 2010; with further explanations on NavLab’s nvector page, and MATLAB code on the Forsvarets forskningsinstitutt nvector 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 Cstyle syntax, JavaScript reads remarkably close to pseudocode: 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 asis 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.
I offer these scripts for free use and adaptation to balance my debt to the opensource infoverse. You are welcome to reuse 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.epytelbavom@oegstpircs.
© 20112022 Chris Veness