From a141977eb6a476b4dfc7e35aaa26531f01f2fff9 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Mon, 27 Nov 2017 11:36:31 -0500 Subject: Update documentation on geodesics replacing an ancient page on geodesic calculations with something more helpful to the beginning user. --- docs/source/apps/geod.rst | 10 +- docs/source/faq.rst | 2 +- docs/source/geodesic.rst | 424 +++++++++++++++++++++++----------------------- 3 files changed, 216 insertions(+), 220 deletions(-) (limited to 'docs/source') diff --git a/docs/source/apps/geod.rst b/docs/source/apps/geod.rst index 4620a6cc..49fa4925 100644 --- a/docs/source/apps/geod.rst +++ b/docs/source/apps/geod.rst @@ -163,11 +163,9 @@ which gives: Further reading *************** -#. `GeographicLib `_ - -#. `C. F. F. Karney, Algorithms for Geodesics, J. Geodesy 87, 43-55 (2013) `_. - `Addendum `_ - -#. `The online geodesic bibliography `_ +#. `GeographicLib `_. +#. C. F. F. Karney, `Algorithms for Geodesics `_, J. Geodesy **87**\ (1), 43–55 (2013); + `addenda `_. +#. `A geodesic bibliography `_. diff --git a/docs/source/faq.rst b/docs/source/faq.rst index 0240df1d..14ceda85 100644 --- a/docs/source/faq.rst +++ b/docs/source/faq.rst @@ -237,7 +237,7 @@ How do I calculate distances/directions on the surface of the earth? -------------------------------------------------------------------------------- These are called geodesic calculations. There is a page about it here: -[wiki:GeodesicCalculations] +:ref:`geodesic`. What is the HEALPix projection and how can I use it? -------------------------------------------------------------------------------- diff --git a/docs/source/geodesic.rst b/docs/source/geodesic.rst index fac3649c..29327ae8 100644 --- a/docs/source/geodesic.rst +++ b/docs/source/geodesic.rst @@ -1,219 +1,217 @@ .. _geodesic: -================================================================================ -Geodesic Calculations -================================================================================ +Geodesic calculations +===================== .. contents:: Contents - :depth: 3 + :depth: 2 :backlinks: none - -Geodesic Calculations --------------------------------------------------------------------------------- - -Geodesic calculations are calculations along lines (great circle) on the -surface of the earth. They can answer questions like: - - * What is the distance between these two points? - * If I travel X meters from point A at bearing phi, where will I be. They are - done in native lat-long coordinates, rather than in projected coordinates. - -Relevant mailing list threads -................................................................................ - - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3361 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3375 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3435 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3588 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/3925 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/4047 - * http://thread.gmane.org/gmane.comp.gis.proj-4.devel/4083 - -Terminology --------------------------------------------------------------------------------- - -The shortest distance on the surface of a solid is generally termed a geodesic, -be it an ellipsoid of revolution, aposphere, etc. On a sphere, the geodesic is -termed a Great Circle. - -HOWEVER, when computing the distance between two points using a projected -coordinate system, that is a conformal projection such as Transverse Mercator, -Oblique Mercator, Normal Mercator, Stereographic, or Lambert Conformal Conic - -that then is a GRID distance which can be converted to an equivalent GEODETIC -distance using the function for "Scale Factor at a Point." The conversion is -then termed "Grid Distance to Geodetic Distance," even though it will not be as -exactly correct as a true ellipsoidal geodesic. Closer to the truth with a TM -than with a Lambert or other conformal projection, but still not exactly "on." - - -So, it can be termed "geodetic distance" or a "geodesic distance," depending -on just how you got there ... - - -The Math --------------------------------------------------------------------------------- - -Spherical Approximation -................................................................................ - -The simplest way to compute geodesics is using a sphere as an approximation for -the earth. This from Mikael Rittri on the Proj mailing list: - - If 1 percent accuracy is enough, I think you can use spherical formulas - with a fixed Earth radius. You can find good formulas in the Aviation - Formulary of Ed Williams, http://williams.best.vwh.net/avform.htm. - - For the fixed Earth radius, I would choose the average of the: - - :math:`c` = radius of curvature at the poles, - - :math:`\frac{b^2}{a}` = radius of curvature in a meridian plane at the equator, - - since these are the extreme values for the local radius of curvature of the - earth ellipsoid. - - If your coordinates are given in WGS84, then - - :math:`c` = 6 399 593.626 m, - - :math:`\frac{b^2}{a}` = 6 335 439.327 m, - - (see http://home.online.no/~sigurdhu/WGS84_Eng.html) so their average is 6,367,516.477 m. - The maximal error for distance calculation should then be less than 0.51 percent. - - When computing the azimuth between two points by the spherical formulas, I - think the maximal error on WGS84 will be 0.2 degrees, at least if the - points are not too far away (less than 1000 km apart, say). The error - should be maximal near the equator, for azimuths near northeast etc. - -I am not sure about the spherical errors for the forward geodetic problem: -point positioning given initial point, distance and azimuth. - -Ellipsoidal Approximation -................................................................................ - -For more accuracy, the earth can be approximated with an ellipsoid, -complicating the math somewhat. See the wikipedia page, `Geodesics on an -ellipsoid `__, for -more information. - -Thaddeus Vincenty's method, April 1975 -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -For a very good procedure to calculate inter point distances see: - -http://www.ngs.noaa.gov/PC_PROD/Inv_Fwd/ (Fortran code, DOS executables, and an online app) - -and algorithm details published in: `Vincenty, T. (1975) `__ - -Javascript code -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Chris Veness has coded Vincenty's formulas as !JavaScript. - -distance: http://www.movable-type.co.uk/scripts/latlong-vincenty.html - -direct: http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html - -C code -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -From Gerald Evenden: a library of the converted NGS Vincenty geodesic procedure -and an application program, 'geodesic'. In the case of a spherical earth -[Snyder1987]_'s preferred equations are used. - -* http://article.gmane.org/gmane.comp.gis.proj-4.devel/3588/ - -The link in this message is broken. The correct URL is -http://home.comcast.net/~gevenden56/proj/ - -Earlier Mr. Evenden had posted to the PROJ.4 mailing list this code for -determination of true distance and respective forward and back azimuths between -two points on the ellipsoid. Good for any pair of points that are not -antipodal. -Later he posted that this was not in fact the translation of NGS FORTRAN code, -but something else. But, for what it's worth, here is the posted code (source -unknown): - -* http://article.gmane.org/gmane.comp.gis.proj-4.devel/3478 - - -PROJ.4 - geod program -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - -The PROJ.4 [wiki:man_geod geod] program can be used for great circle distances -on an ellipsoid. As of proj version 4.9.0, this uses a translation of -GeographicLib::Geodesic (see below) into C. The underlying geodesic -calculation API is exposed as part of the PROJ.4 library (via the geodesic.h -header). Prior to version 4.9.0, the algorithm documented here was used: -` -Paul D. Thomas, 1970 -Spheroidal Geodesics, Reference Systems, and Local Geometry" -U.S. Naval Oceanographic Office, p. 162 -Engineering Library 526.3 T36s - -http://handle.dtic.mil/100.2/AD0703541 - -GeographicLib::Geodesic -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -Charles Karney has written a C++ class to do geodesic calculations and a -utility GeodSolve to call it. See - -* http://geographiclib.sourceforge.net/geod.html - -An online version of GeodSolve is available at - -* http://geographiclib.sourceforge.net/cgi-bin/GeodSolve - -This is an attempt to do geodesic calculations "right", i.e., - -* accurate to round-off (i.e., about 15 nm); -* inverse solution always succeeds (even for near anti-podal points); -* reasonably fast (comparable in speed to Vincenty); -* differential properties of geodesics are computed (these give the scales of - geodesic projections); -* the area between a geodesic and the equator is computed (allowing the - area of geodesic polygons to be found); -* included also is an implementation in terms of elliptic integrals which - can deal with ellipsoids with 0.01 < b/a < 100. - -A JavaScript implementation is included, see - -* `geo-calc `__, - a text interface to geodesic calculations; -* `geod-google `__, - a tool for drawing geodesics on Google Maps. - -Implementations in `Python `__, -`Matlab `__, -`C `__, -`Fortran `__ , and -`Java `__ are also available. - -The algorithms are described in - * C. F. F. Karney, `Algorithms for gedesics `__, - J. Geodesy '''87'''(1), 43-55 (2013), - DOI: `10.1007/s00190-012-0578-z `__; `geo-addenda.html `__. - -Triaxial Ellipsoid -................................................................................ - -A triaxial ellipsoid is a marginally better approximation to the shape of the earth -than an ellipsoid of revolution. -The problem of geodesics on a triaxial ellipsoid was solved by Jacobi in 1838. -For a discussion of this problem see -* http://geographiclib.sourceforge.net/html/triaxial.html -* the wikipedia entry: `Geodesics on a triaxial ellipsoid `__ - -The History --------------------------------------------------------------------------------- - -The bibliography of papers on the geodesic problem for an ellipsoid is -available at - -* http://geographiclib.sourceforge.net/geodesic-papers/biblio.html - -this includes links to online copies of the papers. +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 + +
+ 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 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 +`_ from `GeographicLib +`_. 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 + `_ + +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) + `_, + 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 + `_, + (Teubner, Leipzig, 1880), Chaps. 5–7. +* T. Vincenty, + `Direct and inverse solutions of geodesics on the ellipsoid with + application of nested equations + `_, + Survey Review **23**\ (176), 88–93 (1975). +* J. Danielsen, + `The area under the geodesic + `_, + Survey Review **30**\ (232), 61–66 (1989). +* C. F. F. Karney, + `Algorithms for geodesics + `_, + J. Geodesy **87**\ (1) 43–55 (2013); + `addenda `_. +* C. F. F. Karney, + `Geodesics on an ellipsoid of revolution + `_, + Feb. 2011; + `errata + `_. +* `A geodesic bibliography + `_. +* The wikipedia page, + `Geodesics on an ellipsoid + `_. -- cgit v1.2.3 From a4308d9c529fd913cd8d50813fc913bf80552a38 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Sun, 3 Dec 2017 17:06:02 -0500 Subject: Fixes to geodesic documention: * replace SVG figure by PNG version, * rewrite maths using :math:, * put references in the main reference section, * reformat references for consistency, * put references in alphabetical order, * use USGS URL for Snyder (1987). --- docs/source/geodesic.rst | 187 +++++++++++++++++---------------------------- docs/source/references.rst | 134 +++++++++++++++++++++++++------- 2 files changed, 179 insertions(+), 142 deletions(-) (limited to 'docs/source') diff --git a/docs/source/geodesic.rst b/docs/source/geodesic.rst index 29327ae8..d54212ca 100644 --- a/docs/source/geodesic.rst +++ b/docs/source/geodesic.rst @@ -10,53 +10,52 @@ Geodesic calculations 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 +Consider a ellipsoid of revolution with equatorial radius :math:`a`, polar +semi-axis :math:`b`, and flattening :math:`f=(a−b)/a`. Points on +the surface of the ellipsoid are characterized by their latitude :math:`\phi` +and longitude :math:`\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 -(φ\ :sub:`1`, λ\ :sub:`1`) and (φ\ :sub:`2`, λ\ :sub:`2`) +:math:`(\phi_1,\lambda_1)` and :math:`(\phi_2,\lambda_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`. +:math:`s_{12}` and the geodesic from point 1 to point 2 has forward +azimuths :math:`\alpha_1` and :math:`\alpha_2` at the two end +points. In this figure, we have :math:`\lambda_{12}=\lambda_2-\lambda_1`. .. raw:: html
- Figure from wikipedia + width="320">
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 +Solution of geodesic problems ----------------------------- 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 direct problem — given :math:`\phi_1`, + :math:`\lambda_1`, :math:`\alpha_1`, :math:`s_{12}`, + determine :math:`\phi_2`, :math:`\lambda_2`, :math:`\alpha_2`. -* the inverse problem — given φ\ :sub:`1`, - λ\ :sub:`1`, φ\ :sub:`2`, λ\ :sub:`2`, - determine *s*\ :sub:`12`, α\ :sub:`1`, and - α\ :sub:`2`. +* the inverse problem — given :math:`\phi_1`, + :math:`\lambda_1`, :math:`\phi_2`, :math:`\lambda_2`, + determine :math:`s_{12}`, :math:`\alpha_1`, + :math:`\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 −0.02 < *f* < 0.02. Refer +precision accuracy is maintained provided that +:math:`\lvert f\rvert<\frac1{50}`. Refer to the `application programming interface @@ -70,9 +69,9 @@ 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. +* the shape of ellipsoid is specified by the flattening :math:`f`; this can + be negative to denote a prolate ellipsoid; setting :math:`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. @@ -82,47 +81,41 @@ Additional properties The routines also calculate several other quantities of interest -* *S*\ :sub:`12` is the area between the geodesic from point 1 to +* :math:`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 - (φ\ :sub:`1`,λ\ :sub:`1`), (0,λ\ :sub:`1`), - (0,λ\ :sub:`2`), and - (φ\ :sub:`2`,λ\ :sub:`2`). It is given in + :math:`(\phi_1,\lambda_1)`, :math:`(0,\lambda_1)`, + :math:`(0,\lambda_2)`, and + :math:`(\phi_2,\lambda_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 +* :math:`m_{12}`, the reduced length of the geodesic is defined such + that if the initial azimuth is perturbed by :math:`d\alpha_1` + (radians) then the second point is displaced by :math:`m_{12}\,d\alpha_1` + in the direction perpendicular to the + geodesic. :math:`m_{12}` is given in meters. On a curved surface + the reduced length obeys a symmetry relation, :math:`m_{12}+m_{21}=0`. + On a flat surface, we have :math:`m_{12}=s_{12}`. +* :math:`M_{12}` and :math:`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*\ :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. + :\math`dt`, then they are separated by a distance :math:`M_{12}\,dt` at + point 2. :math:`M_{21}` is defined similarly (with the geodesics + being parallel to one another at point 2). :math:`M_{12}` and + :math:`M_{21}` are dimensionless quantities. On a flat surface, + we have :math:`M_{12}=M_{21}=1`. +* :math:`\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°. + the next is always :math:`180^\circ`. 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` +* :math:`s_{13}=s_{12}+s_{23}`, +* :math:`\sigma_{13}=\sigma_{12}+\sigma_{23}`, +* :math:`S_{13}=S_{12}+S_{23}`, +* :math:`m_{13}=m_{12}M_{23}+m_{23}M_{21}`, +* :math:`M_{13}=M_{12}M_{23}-(1-M_{12}M_{21})m_{23}/m_{12}`, +* :math:`M_{31}=M_{32}M_{21}-(1-M_{23}M_{32})m_{12}/m_{23}`. Multiple shortest geodesics --------------------------- @@ -132,42 +125,41 @@ The shortest distance found by solving the inverse problem is 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 +* :math:`\phi_1=-\phi_2` (with neither point at + a pole). If :math:`\alpha_1=\alpha_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 + :math:`[\alpha_1,\alpha_2]\leftarrow[\alpha_2,\alpha_1]`, + :math:`[M_{12},M_{21}]\leftarrow[M_{21},M_{12}]`, + :math:`S_{12}\leftarrow-S_{12}`. + (This occurs when the longitude difference is near :math:`\pm180^\circ` + for oblate ellipsoids.) +* :math:`\lambda_2=\lambda_1\pm180^\circ` (with + neither point at a pole). If :math:`\alpha_1=0^\circ` or + :math:`\pm180^\circ`, 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 + :math:`[\alpha_1,\alpha_2]\leftarrow[-\alpha_1,-\alpha_2]`, + :math:`S_{12}\leftarrow-S_{12}`. (This occurs when + :math:`\phi_2` is near :math:`-\phi_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 δ. + :math:`[\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,-\delta]`, + for arbitrary :math:`\delta`. (For spheres, this prescription applies when points 1 and 2 are antipodal.) -* *s*\ :sub:`12` = 0 (coincident points). There are infinitely many +* :math:`s_{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 δ. + :math:`[\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,\delta]`, + for arbitrary :math:`\delta`. 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 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 @@ -178,40 +170,5 @@ Vincenty (1975) in the following respects: 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) - `_, - 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 - `_, - (Teubner, Leipzig, 1880), Chaps. 5–7. -* T. Vincenty, - `Direct and inverse solutions of geodesics on the ellipsoid with - application of nested equations - `_, - Survey Review **23**\ (176), 88–93 (1975). -* J. Danielsen, - `The area under the geodesic - `_, - Survey Review **30**\ (232), 61–66 (1989). -* C. F. F. Karney, - `Algorithms for geodesics - `_, - J. Geodesy **87**\ (1) 43–55 (2013); - `addenda `_. -* C. F. F. Karney, - `Geodesics on an ellipsoid of revolution - `_, - Feb. 2011; - `errata - `_. -* `A geodesic bibliography - `_. -* The wikipedia page, - `Geodesics on an ellipsoid - `_. +Additional background material is provided in [GeodesicBib]_, +[GeodesicWiki]_, and [Karney2011]_. diff --git a/docs/source/references.rst b/docs/source/references.rst index 31100ded..e80e8caf 100644 --- a/docs/source/references.rst +++ b/docs/source/references.rst @@ -5,30 +5,110 @@ References ================================================================================ -.. [Evenden1995] Evenden, G. I., 1995, `Cartograpic Projection Procedures for the UNIX Environment - A User's Manual `_ - -.. [Evenden2005] Evenden, G. I., 2005, `libproj4: A Comprehensive Library of Cartographic Projection Functions (Preliminary Draft) `_ - -.. [Steers1970] Steers, J.A., 1970, An introduction to the study of map projections (15th ed.): London, Univ. London Press, p. 229 - -.. [Snyder1987] Snyder. John P. 1987. `Map Projections - A Working Manual `_. US. Geological Survey professional paper; 1395. - -.. [Snyder1993] Snyder, 1993, Flattening the Earth, Chicago and London, The university of Chicago press - -.. [EberHewitt1979] Eber, L.E., and R.P. Hewitt. 1979. `Conversion algorithms for the CALCOFI station grid `__. California Cooperative Oceanic Fisheries Investigations Reports 20:135-137. - -.. [WeberMoore2013] Weber, E.D., and T.J. Moore. 2013. `Corrected Conversion Algorithms For The Calcofi Station Grid And Their Implementation In Several Computer Languages `__. California Cooperative Oceanic Fisheries Investigations Reports 54. - -.. [CalabrettaGreisen2002] M. Calabretta and E. Greisen, 2002, "Representations of celestial coordinates in FITS". Astronomy & Astrophysics 395, 3, 1077–1122. - -.. [ChanONeil1975] F. Chan and E.M.O'Neill, 1975, "Feasibility Study of a Quadrilateralized Spherical Cube Earth Data Base". Tech. Rep. EPRF 2-75 (CSC), Environmental Prediction Research Facility. - -.. [ONeilLaubscher1976] E.M. O'Neill and R.E. Laubscher, 1976, "Extended Studies of a Quadrilateralized Spherical Cube Earth Data Base". Tech. Rep. NEPRF 3-76 (CSC), Naval Environmental Prediction Research Facility. - -.. [LambersKolb2012] M. Lambers and A. Kolb, 2012, "Ellipsoidal Cube Maps for Accurate Rendering of Planetary-Scale Terrain Data", Proc. Pacfic Graphics (Short Papers). - -.. [Zajac1978] A. Zajac, 1978, "Atlas of distribution of vascular plants in Poland (ATPOL)". Taxon 27(5/6), 481–484. - -.. [Komsta2016] L. Komsta, 2016, `ATPOL geobotanical grid revisited – a proposal of coordinate conversion algorithms `__. Annales UMCS Sectio E Agricultura 71(1), 31-37. - -.. [Verey2017] M. Verey, 2017, `Theoretical analysis and practical consequences of adopting an ATPOL grid model as a conical projection, defining the conversion of plane coordinates to the WGS - 84 ellipsoid `__. Fragmenta Floristica et Geobotanica Polonica (preprint submitted). +.. [Bessel1825] F. W. Bessel, 1825, + `The calculation of longitude and latitude from geodesic measurements + `_, + Astron. Nachr. **331**\ (8), 852–861 (2010), + translated by C. F. F. Karney and R. E. Deakin. + +.. [CalabrettaGreisen2002] M. Calabretta and E. Greisen, 2002, + "Representations of celestial coordinates in FITS". + Astronomy & Astrophysics 395, 3, 1077–1122. + +.. [ChanONeil1975] F. Chan and E. M. O'Neill, 1975, + "Feasibility Study of a Quadrilateralized Spherical Cube Earth Data Base", + Tech. Rep. EPRF 2-75 (CSC), Environmental Prediction Research Facility. + +.. [Danielsen1989] J. Danielsen, 1989, + `The area under the geodesic + `_, + Survey Review **30**\ (232), 61–66. + +.. [EberHewitt1979] L. E. Eber and R.P. Hewitt, 1979, + `Conversion algorithms for the CALCOFI station grid + `_, + California Cooperative Oceanic Fisheries Investigations Reports 20:135-137. + +.. [Evenden1995] G. I. Evenden, 1995, + `Cartograpic Projection Procedures for the UNIX Environment - + A User's Manual + `_. + +.. [Evenden2005] G. I. Evenden, 2005, + `libproj4: A Comprehensive Library of Cartographic Projection Functions + (Preliminary Draft) + `_. + +.. [GeodesicBib] `A geodesic bibliography + `_. + +.. [GeodesicWiki] The wikipedia page, + `Geodesics on an ellipsoid + `_. + +.. [Helmert1880] F. R. Helmert, 1880, + `Mathematical and Physical Theories of Higher Geodesy, Vol 1 + `_, + (Teubner, Leipzig), Chaps. 5–7. + +.. [Karney2011] C. F. F. Karney, 2011, + `Geodesics on an ellipsoid of revolution + `_; + `errata + `_. + +.. [Karney2013] C. F. F. Karney, 2013, + `Algorithms for geodesics + `_, + J. Geodesy **87**\ (1) 43–55; + `addenda `_. + +.. [Komsta2016] L. Komsta, 2016, + `ATPOL geobotanical grid revisited - a proposal of coordinate conversion + algorithms + `_, + Annales UMCS Sectio E Agricultura 71(1), 31–37. + +.. [LambersKolb2012] M. Lambers and A. Kolb, 2012, + "Ellipsoidal Cube Maps for Accurate Rendering of Planetary-Scale + Terrain Data", Proc. Pacfic Graphics (Short Papers). + +.. [ONeilLaubscher1976] E. M. O'Neill and R. E. Laubscher, 1976, + "Extended Studies of a Quadrilateralized Spherical Cube Earth Data Base", + Tech. Rep. NEPRF 3-76 (CSC), + Naval Environmental Prediction Research Facility. + +.. [Snyder1987] J. P. Snyder, 1987, + `Map Projections - A Working Manual + `_. + U.S. Geological Survey professional paper 1395. + +.. [Snyder1993] J. P. Snyder, 1993, + Flattening the Earth, Chicago and London, The University of Chicago press. + +.. [Steers1970] J. A. Steers, 1970, + An introduction to the study of map projections (15th ed.), + London Univ. Press, p. 229. + +.. [Verey2017] M. Verey, 2017, + `Theoretical analysis and practical consequences of adopting an ATPOL + grid model as a conical projection, defining the conversion of plane + coordinates to the WGS-84 ellipsoid + `_, + Fragmenta Floristica et Geobotanica Polonica (preprint submitted). + +.. [Vincenty1975] T. Vincenty, 1975, + `Direct and inverse solutions of geodesics on the ellipsoid with + application of nested equations + `_, + Survey Review **23**\ (176), 88–93. + +.. [WeberMoore2013] E. D. Weber and T.J. Moore, 2013, + `Corrected Conversion Algorithms For The Calcofi Station Grid And Their + Implementation In Several Computer Languages + `_, + California Cooperative Oceanic Fisheries Investigations Reports 54. + +.. [Zajac1978] A. Zajac, 1978, + "Atlas of distribution of vascular plants in Poland (ATPOL)", + Taxon 27(5/6), 481–484. -- cgit v1.2.3