1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
|
.. _geodesic:
Geodesic calculations
=====================
.. contents:: Contents
:depth: 2
:backlinks: none
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 φ
and longitude λ. (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
(φ\ :sub:`1`, λ\ :sub:`1`) and (φ\ :sub:`2`, λ\ :sub:`2`)
is called the geodesic. Its length is
*s*\ :sub:`12` and the geodesic from point 1 to point 2 has forward
azimuths α\ :sub:`1` and α\ :sub:`2` at the two end
points. In this figure, we have λ\ :sub:`12` =
λ\ :sub:`2` − λ\ :sub:`1`.
.. raw:: html
<center>
<img src="https://upload.wikimedia.org/wikipedia/commons/c/cb/Geodesic_problem_on_an_ellipsoid.svg"
alt="Figure from wikipedia"
width="250">
</center>
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 programs
-----------------------------
Traditionally two geodesic problems are considered:
* the direct problem — given φ\ :sub:`1`,
λ\ :sub:`1`, α\ :sub:`1`, *s*\ :sub:`12`,
determine φ\ :sub:`2`, λ\ :sub:`2`, and
α\ :sub:`2`.
* the inverse problem — given φ\ :sub:`1`,
λ\ :sub:`1`, φ\ :sub:`2`, λ\ :sub:`2`,
determine *s*\ :sub:`12`, α\ :sub:`1`, and
α\ :sub:`2`.
PROJ incorporates `C library for Geodesics
<https://geographiclib.sourceforge.io/1.49/C/>`_ from `GeographicLib
<https://geographiclib.sourceforge.io>`_. This library provides
routines to solve the direct and inverse geodesic problems. Full double
precision accuracy is maintained provided that −0.02 < *f* < 0.02. Refer
to the
`application programming interface
<https://geographiclib.sourceforge.io/1.49/C/geodesic_8h.html>`_
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, :ref:`geod`\ (1), for performing
simple geodesic calculations.
Additional properties
---------------------
The routines also calculate several other quantities of interest
* *S*\ :sub:`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
(φ\ :sub:`1`,λ\ :sub:`1`), (0,λ\ :sub:`1`),
(0,λ\ :sub:`2`), and
(φ\ :sub:`2`,λ\ :sub:`2`). It is given in
meters\ :sup:`2`.
* *m*\ :sub:`12`, the reduced length of the geodesic is defined such
that if the initial azimuth is perturbed by *d*\ α\ :sub:`1`
(radians) then the second point is displaced by *m*\ :sub:`12`
*d*\ α\ :sub:`1` in the direction perpendicular to the
geodesic. *m*\ :sub:`12` is given in meters. On a curved surface
the reduced length obeys a symmetry relation, *m*\ :sub:`12` +
*m*\ :sub:`21` = 0. On a flat surface, we have *m*\ :sub:`12` =
*s*\ :sub:`12`.
* *M*\ :sub:`12` and *M*\ :sub:`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*\ :sub:`12` *dt* at
point 2. *M*\ :sub:`21` is defined similarly (with the geodesics
being parallel to one another at point 2). *M*\ :sub:`12` and
*M*\ :sub:`21` are dimensionless quantities. On a flat surface,
we have *M*\ :sub:`12` = *M*\ :sub:`21` = 1.
* σ\ :sub:`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°.
If points 1, 2, and 3 lie on a single geodesic, then the following
addition rules hold:
* *s*\ :sub:`13` = *s*\ :sub:`12` + *s*\ :sub:`23`
* σ\ :sub:`13` = σ\ :sub:`12` + σ\ :sub:`23`
* *S*\ :sub:`13` = *S*\ :sub:`12` + *S*\ :sub:`23`
* *m*\ :sub:`13` = *m*\ :sub:`12`\ *M*\ :sub:`23` +
*m*\ :sub:`23`\ *M*\ :sub:`21`
* *M*\ :sub:`13` = *M*\ :sub:`12`\ *M*\ :sub:`23` −
(1 − *M*\ :sub:`12`\ *M*\ :sub:`21`)
*m*\ :sub:`23`/*m*\ :sub:`12`
* *M*\ :sub:`31` = *M*\ :sub:`32`\ *M*\ :sub:`21` −
(1 − *M*\ :sub:`23`\ *M*\ :sub:`32`)
*m*\ :sub:`12`/*m*\ :sub:`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:
* φ\ :sub:`1` = −φ\ :sub:`2` (with neither point at
a pole). If α\ :sub:`1` = α\ :sub:`2`, the geodesic
is unique. Otherwise there are two geodesics and the second one is
obtained by setting
[α\ :sub:`1`,α\ :sub:`2`] ← [α\ :sub:`2`,α\ :sub:`1`],
[*M*\ :sub:`12`,\ *M*\ :sub:`21`] ← [*M*\ :sub:`21`,\ *M*\ :sub:`12`],
*S*\ :sub:`12` ← −\ *S*\ :sub:`12`.
(This occurs when the longitude difference is near ±180° for oblate
ellipsoids.)
* λ\ :sub:`2` = λ\ :sub:`1` ± 180° (with
neither point at a pole). If α\ :sub:`1` = 0° or
±180°, the geodesic is unique. Otherwise there are two
geodesics and the second one is obtained by setting
[α\ :sub:`1`,α\ :sub:`2`] ← [−α\ :sub:`1`,−α\ :sub:`2`],
*S*\ :sub:`12` ← −\ *S*\ :sub:`12`. (This occurs when
φ\ :sub:`2` is near −φ\ :sub:`1` for prolate
ellipsoids.)
* Points 1 and 2 at opposite poles. There are infinitely many
geodesics which can be generated by setting
[α\ :sub:`1`,α\ :sub:`2`] ←
[α\ :sub:`1`,α\ :sub:`2`] + [δ,−δ], for arbitrary δ.
(For spheres, this prescription applies when points 1 and 2 are
antipodal.)
* *s*\ :sub:`12` = 0 (coincident points). There are infinitely many
geodesics which can be generated by setting
[α\ :sub:`1`,α\ :sub:`2`]_←
[α\ :sub:`1`,α\ :sub:`2`]_+ [δ,δ], for
arbitrary δ.
Background
----------
The algorithms implemented by this package are given in Karney (2013)
and are based on Bessel (1825) and Helmert (1880); the algorithm for
areas is based on Danielsen (1989). These improve on the work of
Vincenty (1975) 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.
References
----------
* F. W. Bessel,
`The calculation of longitude and latitude from geodesic measurements (1825)
<https://arxiv.org/abs/0908.1824>`_,
Astron. Nachr. **331**\ (8), 852–861 (2010),
translated by C. F. F. Karney and R. E. Deakin.
* F. R. Helmert,
`Mathematical and Physical Theories of Higher Geodesy, Vol 1
<https://doi.org/10.5281/zenodo.32050>`_,
(Teubner, Leipzig, 1880), Chaps. 5–7.
* T. Vincenty,
`Direct and inverse solutions of geodesics on the ellipsoid with
application of nested equations
<http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf>`_,
Survey Review **23**\ (176), 88–93 (1975).
* J. Danielsen,
`The area under the geodesic
<https://doi.org/10.1179/003962689791474267>`_,
Survey Review **30**\ (232), 61–66 (1989).
* C. F. F. Karney,
`Algorithms for geodesics
<https://doi.org/10.1007/s00190-012-0578-z>`_,
J. Geodesy **87**\ (1) 43–55 (2013);
`addenda <https://geographiclib.sourceforge.io/geod-addenda.html>`_.
* C. F. F. Karney,
`Geodesics on an ellipsoid of revolution
<https://arxiv.org/abs/1102.1215v1>`_,
Feb. 2011;
`errata
<https://geographiclib.sourceforge.io/geod-addenda.html#geod-errata>`_.
* `A geodesic bibliography
<https://geographiclib.sourceforge.io/geodesic-papers/biblio.html>`_.
* The wikipedia page,
`Geodesics on an ellipsoid
<https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid>`_.
|
