diff options
| author | Even Rouault <even.rouault@mines-paris.org> | 2015-02-09 09:01:31 +0000 |
|---|---|---|
| committer | Even Rouault <even.rouault@mines-paris.org> | 2015-02-09 09:01:31 +0000 |
| commit | 862c4682e67ab06b60951b57bf6634789651677a (patch) | |
| tree | 9f44b080df92cc50fef9a0a8cb6855b0594dc966 /src/geodesic.h | |
| parent | 41f59237fa4d2064098fddb2e65729ddf2eacc35 (diff) | |
| download | PROJ-862c4682e67ab06b60951b57bf6634789651677a.tar.gz PROJ-862c4682e67ab06b60951b57bf6634789651677a.zip | |
Fix NaN handling by geod_inverse and incorrect area computation by geod_polygon_addedge if the longitude extent of the added edge exceeds 180 degrees (patch by Charles Karney, #251, #253)
git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2600 4e78687f-474d-0410-85f9-8d5e500ac6b2
Diffstat (limited to 'src/geodesic.h')
| -rw-r--r-- | src/geodesic.h | 118 |
1 files changed, 77 insertions, 41 deletions
diff --git a/src/geodesic.h b/src/geodesic.h index 69985009..6b2afc5a 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -4,10 +4,10 @@ * * This an implementation in C of the geodesic algorithms described in * - C. F. F. Karney, - * <a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> + * <a href="https://dx.doi.org/10.1007/s00190-012-0578-z"> * Algorithms for geodesics</a>, * J. Geodesy <b>87</b>, 43--55 (2013); - * DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> + * DOI: <a href="https://dx.doi.org/10.1007/s00190-012-0578-z"> * 10.1007/s00190-012-0578-z</a>; * addenda: <a href="http://geographiclib.sf.net/geod-addenda.html"> * geod-addenda.html</a>. @@ -75,18 +75,18 @@ * (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: - * - \e lat1 = −\e lat2 (with neither at a pole). If \e azi1 = \e - * azi2, the geodesic is unique. Otherwise there are two geodesics + * - \e lat1 = −\e lat2 (with neither point at a pole). If \e azi1 = + * \e azi2, the geodesic is unique. Otherwise there are two geodesics * and the second one is obtained by setting [\e azi1, \e azi2] = [\e * azi2, \e azi1], [\e M12, \e M21] = [\e M21, \e M12], \e S12 = * −\e S12. (This occurs when the longitude difference is near * ±180° for oblate ellipsoids.) - * - \e lon2 = \e lon1 ± 180° (with neither at a pole). If - * \e azi1 = 0° or ±180°, the geodesic is unique. + * - \e lon2 = \e lon1 ± 180° (with neither point at a pole). + * If \e azi1 = 0° or ±180°, the geodesic is unique. * Otherwise there are two geodesics and the second one is obtained by - * setting [\e azi1, \e azi2] = [−\e azi1, −\e azi2], \e - * S12 = −\e S12. (This occurs when the \e lat2 is near - * −\e lat1 for prolate ellipsoids.) + * setting [\e azi1, \e azi2] = [−\e azi1, −\e azi2], \e S12 + * = −\e S12. (This occurs when \e lat2 is near −\e lat1 for + * prolate ellipsoids.) * - Points 1 and 2 at opposite poles. There are infinitely many * geodesics which can be generated by setting [\e azi1, \e azi2] = * [\e azi1, \e azi2] + [\e d, −\e d], for arbitrary \e d. (For @@ -103,17 +103,17 @@ * to the member functions. Most of the internal comments have been retained. * However, in the process of transcription some documentation has been lost * and the documentation for the C++ classes, GeographicLib::Geodesic, - * GeographicLib::GeodesicLine, and GeographicLib::PolygonArea, should be + * GeographicLib::GeodesicLine, and GeographicLib::PolygonAreaT, should be * consulted. The C++ code remains the "reference implementation". Think * twice about restructuring the internals of the C code since this may make * porting fixes from the C++ code more difficult. * - * Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and licensed + * Copyright (c) Charles Karney (2012-2014) <charles@karney.com> and licensed * under the MIT/X11 License. For more information, see * http://geographiclib.sourceforge.net/ * * This library was distributed with - * <a href="../index.html">GeographicLib</a> 1.32. + * <a href="../index.html">GeographicLib</a> 1.40. **********************************************************************/ #if !defined(GEODESIC_H) @@ -128,7 +128,7 @@ * The minor version of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_MINOR 32 +#define GEODESIC_VERSION_MINOR 40 /** * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) @@ -381,11 +381,13 @@ extern "C" { * @param[in] lat1 latitude of point 1 (degrees). * @param[in] lon1 longitude of point 1 (degrees). * @param[in] azi1 azimuth at point 1 (degrees). - * @param[in] arcmode flag determining the meaning of the \e - * s12_a12. - * @param[in] s12_a12 if \e arcmode is 0, this is the distance between - * point 1 and point 2 (meters); otherwise it is the arc length between - * point 1 and point 2 (degrees); it can be negative. + * @param[in] flags bitor'ed combination of geod_flags(); \e flags & + * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags & + * GEOD_LONG_NOWRAP prevents the value of \e lon2 being wrapped into + * the range [−180°, 180°). + * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the distance + * between point 1 and point 2 (meters); otherwise it is the arc length + * between point 1 and point 2 (degrees); it can be negative. * @param[out] plat2 pointer to the latitude of point 2 (degrees). * @param[out] plon2 pointer to the longitude of point 2 (degrees). * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). @@ -402,14 +404,21 @@ extern "C" { * * \e g must have been initialized with a call to geod_init(). \e lat1 * should be in the range [−90°, 90°]; \e lon1 and \e azi1 - * should be in the range [−540°, 540°). The function value \e - * a12 equals \e s12_a12 is \e arcmode is non-zero. Any of the "return" - * arguments \e plat2, etc., may be replaced by 0, if you do not need some - * quantities computed. + * should be in the range [−540°, 540°). The function + * value \e a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the + * "return" arguments \e plat2, etc., may be replaced by 0, if you do not + * need some quantities computed. + * + * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 − + * \e lon1 indicates how many times the geodesic wrapped around the + * ellipsoid. Because \e lon2 might be outside the normal allowed range + * for longitudes, [−540°, 540°), be sure to normalize it, + * e.g., with fmod(\e lon2, 360.0) before using it in subsequent + * calculations **********************************************************************/ double geod_gendirect(const struct geod_geodesic* g, double lat1, double lon1, double azi1, - int arcmode, double s12_a12, + unsigned flags, double s12_a12, double* plat2, double* plon2, double* pazi2, double* ps12, double* pm12, double* pM12, double* pM21, double* pS12); @@ -453,12 +462,16 @@ extern "C" { * * @param[in] l a pointer to the geod_geodesicline object specifying the * geodesic line. - * @param[in] arcmode flag determining the meaning of the second parameter; - * if arcmode is 0, then \e l must have been initialized with \e caps |= + * @param[in] flags bitor'ed combination of geod_flags(); \e flags & + * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags & + * GEOD_LONG_NOWRAP prevents the value of \e lon2 being wrapped into + * the range [−180°, 180°); if \e flags & GEOD_ARCMODE is + * 0, then \e l must have been initialized with \e caps |= * GEOD_DISTANCE_IN. - * @param[in] s12_a12 if \e arcmode is 0, this is the distance between - * point 1 and point 2 (meters); otherwise it is the arc length between - * point 1 and point 2 (degrees); it can be negative. + * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the + * distance between point 1 and point 2 (meters); otherwise it is the + * arc length between point 1 and point 2 (degrees); it can be + * negative. * @param[out] plat2 pointer to the latitude of point 2 (degrees). * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires * that \e l was initialized with \e caps |= GEOD_LONGITUDE. @@ -480,11 +493,18 @@ extern "C" { * @return \e a12 arc length of between point 1 and point 2 (degrees). * * \e l must have been initialized with a call to geod_lineinit() with \e - * caps |= GEOD_DISTANCE_IN. The values of \e lon2 and \e azi2 returned are - * in the range [−180°, 180°). Any of the "return" arguments - * \e plat2, etc., may be replaced by 0, if you do not need some quantities - * computed. Requesting a value which \e l is not capable of computing is - * not an error; the corresponding argument will not be altered. + * caps |= GEOD_DISTANCE_IN. The value \e azi2 returned is in the range + * [−180°, 180°). Any of the "return" arguments \e plat2, + * etc., may be replaced by 0, if you do not need some quantities + * computed. Requesting a value which \e l is not capable of computing + * is not an error; the corresponding argument will not be altered. + * + * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 − + * \e lon1 indicates how many times the geodesic wrapped around the + * ellipsoid. Because \e lon2 might be outside the normal allowed range + * for longitudes, [−540°, 540°), be sure to normalize it, + * e.g., with fmod(\e lon2, 360.0) before using it in subsequent + * calculations * * Example, compute way points between JFK and Singapore Changi Airport * using geod_genposition(). In this example, the points are evenly space in @@ -494,7 +514,7 @@ extern "C" { @code struct geod_geodesic g; struct geod_geodesicline l; - double a12, azi1, lat[101],lon[101]; + double a12, azi1, lat[101], lon[101]; int i; geod_init(&g, 6378137, 1/298.257223563); a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99, @@ -508,7 +528,7 @@ extern "C" { @endcode **********************************************************************/ double geod_genposition(const struct geod_geodesicline* l, - int arcmode, double s12_a12, + unsigned flags, double s12_a12, double* plat2, double* plon2, double* pazi2, double* ps12, double* pm12, double* pM12, double* pM21, @@ -526,6 +546,10 @@ extern "C" { * polylinep is non-zero, then the vertices and edges define a polyline and * only the perimeter is returned by geod_polygon_compute(). * + * The area and perimeter are accumulated at two times the standard floating + * point precision to guard against the loss of accuracy with many-sided + * polygons. At any point you can ask for the perimeter and area so far. + * * An example of the use of this function is given in the documentation for * geod_polygon_compute(). **********************************************************************/ @@ -592,10 +616,12 @@ extern "C" { * polyline (meters). * @return the number of points. * - * Only simple polygons (which are not self-intersecting) are allowed. - * There's no need to "close" the polygon by repeating the first vertex. Set - * \e pA or \e pP to zero, if you do not want the corresponding quantity - * returned. + * The area and perimeter are accumulated at two times the standard floating + * point precision to guard against the loss of accuracy with many-sided + * polygons. Only simple polygons (which are not self-intersecting) are + * allowed. There's no need to "close" the polygon by repeating the first + * vertex. Set \e pA or \e pP to zero, if you do not want the corresponding + * quantity returned. * * Example, compute the perimeter and area of the geodesic triangle with * vertices (0°N,0°E), (0°N,90°E), (90°N,0°E). @@ -704,7 +730,7 @@ extern "C" { * area returned is signed with counter-clockwise traversal being treated as * positive. * - * Example, compute the area of Antarctic: + * Example, compute the area of Antarctica: @code double lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7, @@ -723,7 +749,7 @@ extern "C" { double* pA, double* pP); /** - * mask values for the the \e caps argument to geod_lineinit(). + * mask values for the \e caps argument to geod_lineinit(). **********************************************************************/ enum geod_mask { GEOD_NONE = 0U, /**< Calculate nothing */ @@ -738,6 +764,16 @@ extern "C" { GEOD_ALL = 0x7F80U| 0x1FU /**< Calculate everything */ }; + /** + * flag values for the \e flags argument to geod_gendirect() and + * geod_genposition() + **********************************************************************/ + enum geod_flags { + GEOD_NOFLAGS = 0U, /**< No flags */ + GEOD_ARCMODE = 1U<<0, /**< Position given in terms of arc distance */ + GEOD_LONG_NOWRAP = 1U<<15 /**< Don't wrap longitude */ + }; + #if defined(__cplusplus) } #endif |