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 toVector()
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 toVector()
), 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 var eX = P2.minus(P1).unit(); // unit vector in x direction P1>P2 var i = eX.dot(P3.minus(P1)); // signed magnitude of x component of P1>P3 var eY = P3.minus(P1).minus(eX.times(i)).unit(); // unit vector in y direction var d = P2.minus(P1).length(); // distance P1>P2 var j = eY.dot(P3.minus(P1)); // signed magnitude of y component of P1>P3 var x = (r1*r1  r2*r2 + d*d) / (2*d); // x component of P1 > intersection var y = (r1*r1  r3*r3 + i*i + j*j) / (2*j)  x*i/j; // y component of P1 > intersection var eZ = eX.cross(eY); // unit vector perpendicular to plane var z = Math.sqrt(r1*r1  x*x  y*y); // z will be NaN for no intersections var Pi = P1.plus(eX.times(x)).plus(eY.times(y)); // note don't use z component; assume points at same height
... 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.
See below for the JavaScript source code, also available on GitHub. 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.
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.
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 would like to show your appreciation and support continued development of these scripts, I would most gratefully accept donations.
If you have any queries or find any problems, contact me at ku.oc.epytelbavom@oegstpircs.
© 20112016 Chris Veness