+

+ +

+

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

+

+
+ +

+

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 application programming interface +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².
\(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\).

+

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

+

+ + +

+