Movable Type Scripts

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 by¹

(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 version¹

• d = R · atan2(|a×b|, a·b)

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

See distanceTo() in the source code below.

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).¹ To get angles greater than ±π/2 (±90°), the vector a can be used as a reference to determine the direction of c₁×c₂ hence the sign of its length.

• N = {0,0,1}
• c₁ = a×b
• c₂ = a×N
• sinθ = |c₁×c₂| · sgn(c₁×c₂ · a)
• cosθ = c₁·c₂
• θ = atan2(sinθ, cosθ)

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:

• m̂ = a + b

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:

• m̂ = a·(1−f) + b·f

(or equivalently m̂ = a + f·(b−a)).

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

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

• N = {0,0,1} – vector representing north pole
• d̂_e = N×a – east vector at a
• d_n = a×d_e – north vector at a
• d = d_n·cosθ + d_e·sinθ – direction vector in dir’n of θ
• b = a·cosδ + d·sinδ

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.

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.

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

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:

• N = {0,0,1} – vector representing north pole
•  
• d̂_ae = N×a – east vector at a
• d_an = a×d_ae – north vector at a
• d_a = d_an·cosθ_a + d_ae·sinθ_a – direction vector from a
• c_a = a×d_a – great circle from a
•  
• d̂_be = N×b – east vector at b
• d_bn = b×d_be – north vector at b
• d_b = d_bn·cosθ_b + d_be·sinθ_b – direction vector from b
• c_b = b×d_b – great circle from b
•  
• c = c_a×c_b – n-vector of destination point

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 p₁, p₂, p₃, to the three points (using toNvector()), and (angular) distances r₁, r₂, r₃ from the point of unknown location to those points, then

• ê_x = p₂ − p₁
• i = ê_x ⋅ (p₃ − p₁)
• ê_y = p₃ − p₁ − ê_x⋅i
• d = | p₂ − p₁ |
• j = ê_y ⋅ (p₃ − p₁)
• x = (r₁² − r₂² + d²) / 2⋅d
• y = (r₁² − r₃² + i² + j²) / 2⋅j − x⋅i/j
• ê_z = ê_x × ê_y
• z = ±√r₁² − x² − y²
• p_i = p₁ + ê_x ⋅ x + ê_y⋅y + ê_z⋅z

z will have 0, 1, or 2 solutions depending whether the spheres defined by r₁, r₂, r₃ 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:

• E = Σθᵢ − (n−2)·π
• A = E·R²

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 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

• ∑ a×b / |a×b| ⋅ θ_ab/2, for each pair of consecutive vertices a,b

(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.

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.