Geodesic algorithms

Ellipsoid

An ellipsoid is defined by its primary parameters, i.e. its flattening f and the length of its major axis a.

The WGS-84 ellipsoid is an oblate ellipsoid, the shape formed by rotating an ellipse around its minor axis, see Figure 1.

Primary parameters

Secondary parameters

Secondary parameters are derived from the primary parameters:

Geodesic paths

An oblate ellipsoid is a sphere flattened at its poles.

Meridional geodesic segments on an oblate ellipsoid (e.g. paths that run along lines of longitude) directly correspond to great circle arcs on a sphere.

However, non-meridional geodesic paths lie further from the equator than their equivalent great circles.

The most extreme examples of geodesic segments diverging from great circle arcs is between nearly antipodal points, see Figure 2.

Symbols

A geodesic can be defined by its azimuth \(\alpha_0\) at the point where it crosses the equator, see E in Figure 3.

Auxiliary sphere

A common method to solve geodesic problems problems uses an auxiliary sphere, since the path of a great circle on an auxiliary sphere corresponds to the path of a geodesic on an ellipsoid, see Figure 4.

C. F. F. Karney published accurate geodesic algorithms based on the use of an auxiliary sphere in Algorithms for Geodesics [2].

Geodesic values are mapped to their equivalent values on an auxiliary sphere:

• Geodetic Latitude \(\phi \longleftrightarrow \beta\) Parametric Latitude
• Azimuth \(\alpha \longleftrightarrow \alpha\) Azimuth
• Geodesic Distance \(s \longleftrightarrow \sigma\) Spherical Distance
• Geographic Longitude \(\lambda \longleftrightarrow \omega\) Spherical Longitude

The following table shows parameters on an ellipsoid and their equivalent parameters on an auxiliary sphere:

Table 4 Auxiliary sphere equivalent parameters.

Latitude

Conventional geodetic latitude \(\phi\) is measured between the equatorial plane and a normal to surface of the ellipsoid.
Parametric latitude \(\beta\) is measured from the centre of the Earth to a point projected onto an auxiliary sphere, see Figure 5.
Their conversion is given by:

Eq 1

Azimuth

Azimuth is preserved between geodesic and spherical paths.

Clairaut (1735) discovered the following invariant for a geodesic:

Eq 2

\(\alpha_0\) is the equatorial azimuth, i.e. the azimuth where a geodesic crosses the equator, see E in Figures 3 & 4.

Note: Eq 2 can also be used to calculate the maximum/minimum latitude of a geodesic path, \(\beta_0\).

Equatorial azimuth defines a geodesic path. It is used to calculate the expansion parameter \(\epsilon\) used in series expansions to calculate geodesic distance and longitude:

Eq 3, [2] equation (9)

Eq 4, [2] equation (16)

Distance and longitude

Conversion between spherical and geodesic distances and longitudes is performed relative to the point where the geodesic path crosses the equator, E in Figure 6.

Eq 5, \(f\) is [2] equation (15)

Eq 6, \(f'\) is [2] equation (20)

Eq 7, \(g\) is [2] equation (23)

The coefficients used in equations \(f, f'\) and \(g\) are defined in [2].

Inverse problem

The inverse problem of geodesy is to find the azimuth and distance between two points given their initial geodetic latitudes: \(\phi_1\) and \(\phi_2\) and longitude difference \(\lambda_{12}\), see Figure 7.

The distance and longitude equations derive their coefficients from the equatorial azimuth, \(\alpha_0\) using Eq 3 and 4.

Initial azimuth

The fundamental issue of the inverse problem is that azimuth is not known; only the geodesic longitude difference, which depends on the start azimuth.

Karney’s solution to the inverse problem uses Newton’s method to converge to an accurate value for the start azimuth from an estimated initial azimuth, see [2] section 4 and [4] Appendix A: Addenda to Karney (2013).

The key to Karney’s solution is to use a relatively accurate estimate for the initial azimuth, see [2] section 5.

Geodesic intersections

C. F. F. Karney’s 2023 paper Geodesic Intersections [3] is an elaboration of Baselga and Martinez Llario’s 2017 paper: Intersection and point-to-line solutions for geodesics on the ellipsoid [4].

Karney’s paper provides complete solutions to the intersection of geodesic paths. Not just finding the closest intersection, but all possible intersections between two geodesic paths. It elaborates the spherical trigonometry of Baselga and Martinez Llario’s paper and fixes a few issues, i.e.: finding the closest intersection and handling coincident geodesics.

Like Baselga and Martinez Llario’s paper, Karney’s paper calculates distances to intersections along both geodesic paths using spherical trigonometry. Iterating until the geodesic distance between the intersection points calculated along each geodesic path is within a given threshold.

Geodesic segment intersections

We are interested in the closest intersection of two geodesic segments (equivalent to line segments and great circle arcs); we are not concerned about any other geodesic path intersections.

By representing geodesic segments as great circle arcs on a unit sphere, our algorithm uses spherical vector geometry instead of spherical trigonometry to calculate distances along geodesic segments to their closest intersection point.

A unit sphere is obtained by stretching the ellipsoid along its axis of rotation (i.e. the z axis). Great circles on a sphere correspond to great ellipses on an ellipsoid NOT geodesics, so the algorithm is iterative.

First we calculate a tentative intersection point (see [3] Appendix A) from the points and poles at the arc mid-points:

• an intersection point is calculated from the poles;
• the arcs centroid is calculated from their mid-points;
• the closest intersection point to the centroid is chosen, see Figure 9.

Then we:

1. calculate great circle distances along each arc to the tentative intersection point;
2. calculate new points and poles along each arc at the distances;
3. calculate the distance between the new points;
4. if the new points are NOT close enough to each other, a new tentative intersection point is calculated from the poles and we return to step 1.
5. Otherwise, the algorithm stops because it has found the intersection.

The intersection point can be calculated from the great circle distance along either arc see, Figure 10 .

Note: the test for whether the points are close enough can just compare the square of the Euclidean distance against the threshold. Calculating the square of the Euclidean distance is much quicker than calculating the geodesic distance.

Coincident geodesics

If geodesic segments lie on the same geodesic path, they are said to be coincident. Coincident geodesics do not intersect, they overlap, i.e. they effectively intersect everywhere

Karney [3] handles coincident geodesics by calculating a new geodesic segment between the geodesic start points, forming a triangle, see Figure 11.

If the triangle is degenerate, i.e. all three points are inline, then the geodesics are coincident.

We also construct a geodesic segment between the geodesic segment start points and determine whether the geodesic segments are coincident from their azimuth differences: \(\mu_x\) and \(\mu_y\).

Along track and cross track distance

Both Karney [3] Baselga and Martinez Llario [4] describe calculating along track and cross track distance as solving the “interception” problem, see [3] Appendix B.

Like with intersections, we calculate a tentative interception point from the along track distance of the point along the great circle arc on the unit sphere.

Then, we iterate. Calculating the new along track distance to the interception point from the last calculated interception point using the new great circle pole vector.

The interception point is found when the new along track distance is within the required threshold.

The geodesic along track distance is the sum of the calculated great circle along track distances on the unit sphere converted to geodesic distance.

The cross track distance is the geodesic distance between the point and interception point.

The sign of the cross track distance is given by the dot product of the point and the great circle pole vector at the interception point, see cross track distance.

References

1 World Geodetic System 1984 Manual

ICAO Doc 9674 2nd Edition (2002).

2 Karney 2013

C. F. F. Karney, 2013, Algorithms for geodesics, doi:10.1007/s00190-012-0578-z

3 Karney 2023

C. F. F. Karney, 2023, Geodesic intersections, doi:10.1061/JSUED2.SUENG-1483

4 Baselga and Martinez Llario 2017

Baselga, Sergio & Martinez-Llario, Jose, 2017, Intersection and point-to-line solutions for geodesics on the ellipsoid, doi:10.1007/s11200-017-1020-z