Geodesic calculations

Introduction

Consider an ellipsoid of revolution with equatorial radius \(a\), polar semi-axis \(b\), and flattening \(f=(a-b)/a\). Points on the surface of the ellipsoid are characterized by their latitude \(\phi\) and longitude \(\lambda\). (Note that latitude here means the geographical latitude, the angle between the normal to the ellipsoid and the equatorial plane).

The shortest path between two points on the ellipsoid at \((\phi_1,\lambda_1)\) and \((\phi_2,\lambda_2)\) is called the geodesic. Its length is \(s_{12}\) and the geodesic from point 1 to point 2 has forward azimuths \(\alpha_1\) and \(\alpha_2\) at the two end points. In this figure, we have \(\lambda_{12}=\lambda_2-\lambda_1\).

Figure from wikipedia

A geodesic can be extended indefinitely by requiring that any sufficiently small segment is a shortest path; geodesics are also the straightest curves on the surface.

Solution of geodesic problems

Traditionally two geodesic problems are considered:

  • the direct problem — given \(\phi_1\), \(\lambda_1\), \(\alpha_1\), \(s_{12}\), determine \(\phi_2\), \(\lambda_2\), \(\alpha_2\).

  • the inverse problem — given \(\phi_1\), \(\lambda_1\), \(\phi_2\), \(\lambda_2\), determine \(s_{12}\), \(\alpha_1\), \(\alpha_2\).

PROJ incorporates C library for Geodesics. This library provides routines to solve the direct and inverse geodesic problems. Full double precision accuracy is maintained provided that \(\lvert f\rvert<\frac1{50}\).

The interface to the geodesic routines differ in two respects from the rest of PROJ:

  • angles (latitudes, longitudes, and azimuths) are in degrees (instead of in radians);

  • the shape of ellipsoid is specified by the flattening \(f\); this can be negative to denote a prolate ellipsoid; setting \(f=0\) corresponds to a sphere, in which case the geodesic becomes a great circle.

PROJ also includes a command line tool, geod(1), for performing simple geodesic calculations.

Additional properties

The routines also calculate several other quantities of interest

  • \(S_{12}\) is the area between the geodesic from point 1 to point 2 and the equator; i.e., it is the area, measured counter-clockwise, of the quadrilateral with corners \((\phi_1,\lambda_1)\), \((0,\lambda_1)\), \((0,\lambda_2)\), and \((\phi_2,\lambda_2)\). It is given in meters2.

  • \(m_{12}\), the reduced length of the geodesic is defined such that if the initial azimuth is perturbed by \(d\alpha_1\) (radians) then the second point is displaced by \(m_{12}\,d\alpha_1\) in the direction perpendicular to the geodesic. \(m_{12}\) is given in meters. On a curved surface the reduced length obeys a symmetry relation, \(m_{12}+m_{21}=0\). On a flat surface, we have \(m_{12}=s_{12}\).

  • \(M_{12}\) and \(M_{21}\) are geodesic scales. If two geodesics are parallel at point 1 and separated by a small distance \(dt\), then they are separated by a distance \(M_{12}\,dt\) at point 2. \(M_{21}\) is defined similarly (with the geodesics being parallel to one another at point 2). \(M_{12}\) and \(M_{21}\) are dimensionless quantities. On a flat surface, we have \(M_{12}=M_{21}=1\).

  • \(\sigma_{12}\) is the arc length on the auxiliary sphere. This is a construct for converting the problem to one in spherical trigonometry. The spherical arc length from one equator crossing to the next is always \(180^\circ\).

If points 1, 2, and 3 lie on a single geodesic, then the following addition rules hold:

  • \(s_{13}=s_{12}+s_{23}\),

  • \(\sigma_{13}=\sigma_{12}+\sigma_{23}\),

  • \(S_{13}=S_{12}+S_{23}\),

  • \(m_{13}=m_{12}M_{23}+m_{23}M_{21}\),

  • \(M_{13}=M_{12}M_{23}-(1-M_{12}M_{21})m_{23}/m_{12}\),

  • \(M_{31}=M_{32}M_{21}-(1-M_{23}M_{32})m_{12}/m_{23}\).

Multiple shortest geodesics

The shortest distance found by solving the inverse problem is (obviously) uniquely defined. However, in a few special cases there are multiple azimuths which yield the same shortest distance. Here is a catalog of those cases:

  • \(\phi_1=-\phi_2\) (with neither point at a pole). If \(\alpha_1=\alpha_2\), the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting \([\alpha_1,\alpha_2]\leftarrow[\alpha_2,\alpha_1]\), \([M_{12},M_{21}]\leftarrow[M_{21},M_{12}]\), \(S_{12}\leftarrow-S_{12}\). (This occurs when the longitude difference is near \(\pm180^\circ\) for oblate ellipsoids.)

  • \(\lambda_2=\lambda_1\pm180^\circ\) (with neither point at a pole). If \(\alpha_1=0^\circ\) or \(\pm180^\circ\), the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting \([\alpha_1,\alpha_2]\leftarrow[-\alpha_1,-\alpha_2]\), \(S_{12}\leftarrow-S_{12}\). (This occurs when \(\phi_2\) is near \(-\phi_1\) for prolate ellipsoids.)

  • Points 1 and 2 at opposite poles. There are infinitely many geodesics which can be generated by setting \([\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,-\delta]\), for arbitrary \(\delta\). (For spheres, this prescription applies when points 1 and 2 are antipodal.)

  • \(s_{12}=0\) (coincident points). There are infinitely many geodesics which can be generated by setting \([\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,\delta]\), for arbitrary \(\delta\).

Area of a polygon

The area of a geodesic polygon can be determined by summing \(-S_{12}\) for successive edges of the polygon (\(S_{12}\) is negated so that clockwise traversal of a polygon gives a positive area). However, if the polygon encircles a pole, the sum must be adjusted by \(\pm A/2\), where \(A\) is the area of the full ellipsoid, with the sign chosen to place the result in \((-A/2, A/2]\).

Background

The algorithms implemented by this package are given in [Karney2013] (addenda) and are based on [Bessel1825] and [Helmert1880]; the algorithm for areas is based on [Danielsen1989]. These improve on the work of [Vincenty1975] in the following respects:

  • The results are accurate to round-off for terrestrial ellipsoids (the error in the distance is less than 15 nanometers, compared to 0.1 mm for Vincenty).

  • The solution of the inverse problem is always found. (Vincenty's method fails to converge for nearly antipodal points.)

  • The routines calculate differential and integral properties of a geodesic. This allows, for example, the area of a geodesic polygon to be computed.

Additional background material is provided in GeographicLib's geodesic bibliography, Wikipedia's article "Geodesics on an ellipsoid", and [Karney2011] (errata).