# 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

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.

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

## Background¶

The algorithms implemented by this package are given in (addenda) and are based on and ; the algorithm for areas is based on . These improve on the work of 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 (errata).