# Geodesic calculations¶

Contents

## Introduction¶

Consider a 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\).

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 from GeographicLib. 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}\). Refer to the

for full documentation. A brief summary of the routines is given in geodesic(3).

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
meters
^{2}. - \(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 :math`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\).

## Background¶

The algorithms implemented by this package are given in [Karney2013] 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 then 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 [GeodesicBib], [GeodesicWiki], and [Karney2011].