From 1ee854e1cf70a48f7aee01ce1f56864d7f56f22f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matej=20Fran=C4=8De=C5=A1kin?= Date: Sun, 10 Nov 2019 16:58:08 +0100 Subject: [PATCH] MGRS - Added GeographicLib --- qgroundcontrol.pro | 22 +- src/Geo/Constants.hpp | 417 +++++++++++++++++++ src/Geo/MGRS.cpp | 463 ++++++++++++++++++++++ src/Geo/MGRS.hpp | 361 +++++++++++++++++ src/Geo/Math.cpp | 597 ++++++++++++++++++++++++++++ src/Geo/Math.hpp | 647 ++++++++++++++++++++++++++++++ src/Geo/PolarStereographic.cpp | 109 +++++ src/Geo/PolarStereographic.hpp | 160 ++++++++ src/{ => Geo}/QGCGeo.cc | 0 src/{ => Geo}/QGCGeo.h | 0 src/Geo/TransverseMercator.cpp | 596 ++++++++++++++++++++++++++++ src/Geo/TransverseMercator.hpp | 206 ++++++++++ src/{ => Geo}/UTM.cpp | 0 src/{ => Geo}/UTM.h | 0 src/Geo/UTMUPS.cpp | 297 ++++++++++++++ src/Geo/UTMUPS.hpp | 428 ++++++++++++++++++++ src/Geo/Utility.cpp | 61 +++ src/Geo/Utility.h | 704 +++++++++++++++++++++++++++++++++ src/Geo/Utility.hpp | 704 +++++++++++++++++++++++++++++++++ 19 files changed, 5768 insertions(+), 4 deletions(-) create mode 100644 src/Geo/Constants.hpp create mode 100644 src/Geo/MGRS.cpp create mode 100644 src/Geo/MGRS.hpp create mode 100644 src/Geo/Math.cpp create mode 100644 src/Geo/Math.hpp create mode 100644 src/Geo/PolarStereographic.cpp create mode 100644 src/Geo/PolarStereographic.hpp rename src/{ => Geo}/QGCGeo.cc (100%) rename src/{ => Geo}/QGCGeo.h (100%) create mode 100644 src/Geo/TransverseMercator.cpp create mode 100644 src/Geo/TransverseMercator.hpp rename src/{ => Geo}/UTM.cpp (100%) rename src/{ => Geo}/UTM.h (100%) create mode 100644 src/Geo/UTMUPS.cpp create mode 100644 src/Geo/UTMUPS.hpp create mode 100644 src/Geo/Utility.cpp create mode 100644 src/Geo/Utility.h create mode 100644 src/Geo/Utility.hpp diff --git a/qgroundcontrol.pro b/qgroundcontrol.pro index 536f6155b..5df3a1ab2 100644 --- a/qgroundcontrol.pro +++ b/qgroundcontrol.pro @@ -406,6 +406,7 @@ INCLUDEPATH += \ src/FlightMap \ src/FlightMap/Widgets \ src/FollowMe \ + src/Geo \ src/GPS \ src/Joystick \ src/PlanView \ @@ -627,12 +628,20 @@ HEADERS += \ src/MissionManager/VisualMissionItem.h \ src/PositionManager/PositionManager.h \ src/PositionManager/SimulatedPosition.h \ + src/Geo/QGCGeo.h \ + src/Geo/UTM.h \ + src/Geo/Constants.hpp \ + src/Geo/Math.hpp \ + src/Geo/Utility.hpp \ + src/Geo/UTMUPS.hpp \ + src/Geo/MGRS.hpp \ + src/Geo/TransverseMercator.hpp \ + src/Geo/PolarStereographic.hpp \ src/QGC.h \ src/QGCApplication.h \ src/QGCComboBox.h \ src/QGCConfig.h \ src/QGCFileDownload.h \ - src/QGCGeo.h \ src/QGCLoggingCategory.h \ src/QGCMapPalette.h \ src/QGCPalette.h \ @@ -689,7 +698,6 @@ HEADERS += \ src/uas/UAS.h \ src/uas/UASInterface.h \ src/uas/UASMessageHandler.h \ - src/UTM.h \ src/AnalyzeView/GeoTagController.h \ src/AnalyzeView/ExifParser.h \ src/uas/FileManager.h \ @@ -857,11 +865,18 @@ SOURCES += \ src/MissionManager/VisualMissionItem.cc \ src/PositionManager/PositionManager.cpp \ src/PositionManager/SimulatedPosition.cc \ + src/Geo/QGCGeo.cc \ + src/Geo/UTM.cpp \ + src/Geo/Math.cpp \ + src/Geo/Utility.cpp \ + src/Geo/UTMUPS.cpp \ + src/Geo/MGRS.cpp \ + src/Geo/TransverseMercator.cpp \ + src/Geo/PolarStereographic.cpp \ src/QGC.cc \ src/QGCApplication.cc \ src/QGCComboBox.cc \ src/QGCFileDownload.cc \ - src/QGCGeo.cc \ src/QGCLoggingCategory.cc \ src/QGCMapPalette.cc \ src/QGCPalette.cc \ @@ -917,7 +932,6 @@ SOURCES += \ src/main.cc \ src/uas/UAS.cc \ src/uas/UASMessageHandler.cc \ - src/UTM.cpp \ src/AnalyzeView/GeoTagController.cc \ src/AnalyzeView/ExifParser.cc \ src/uas/FileManager.cc \ diff --git a/src/Geo/Constants.hpp b/src/Geo/Constants.hpp new file mode 100644 index 000000000..adb2d46e8 --- /dev/null +++ b/src/Geo/Constants.hpp @@ -0,0 +1,417 @@ +/** + * \file Constants.hpp + * \brief Header for GeographicLib::Constants class + * + * Copyright (c) Charles Karney (2008-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#if !defined(GEOGRAPHICLIB_CONSTANTS_HPP) +#define GEOGRAPHICLIB_CONSTANTS_HPP 1 + +// This will be overwritten by ./configure + +#define GEOGRAPHICLIB_VERSION_STRING "1.50" +#define GEOGRAPHICLIB_VERSION_MAJOR 1 +#define GEOGRAPHICLIB_VERSION_MINOR 50 +#define GEOGRAPHICLIB_VERSION_PATCH 0 + +// Undefine HAVE_LONG_DOUBLE if this type is unknown to the compiler +#define GEOGRAPHICLIB_HAVE_LONG_DOUBLE 1 + +typedef double real; + +// Define WORDS_BIGENDIAN to be 1 if your machine is big endian +/* #undef WORDS_BIGENDIAN */ + +#include "Math.hpp" + +/** + * @relates GeographicLib::Constants + * Pack the version components into a single integer. Users should not rely on + * this particular packing of the components of the version number; see the + * documentation for GEOGRAPHICLIB_VERSION, below. + **********************************************************************/ +#define GEOGRAPHICLIB_VERSION_NUM(a,b,c) ((((a) * 10000 + (b)) * 100) + (c)) + +/** + * @relates GeographicLib::Constants + * The version of GeographicLib as a single integer, packed as MMmmmmpp where + * MM is the major version, mmmm is the minor version, and pp is the patch + * level. Users should not rely on this particular packing of the components + * of the version number. Instead they should use a test such as \code + #if GEOGRAPHICLIB_VERSION >= GEOGRAPHICLIB_VERSION_NUM(1,37,0) + ... + #endif + * \endcode + **********************************************************************/ +#define GEOGRAPHICLIB_VERSION \ + GEOGRAPHICLIB_VERSION_NUM(GEOGRAPHICLIB_VERSION_MAJOR, \ + GEOGRAPHICLIB_VERSION_MINOR, \ + GEOGRAPHICLIB_VERSION_PATCH) + +/** + * @relates GeographicLib::Constants + * Is the C++11 static_assert available? + **********************************************************************/ +#if !defined(GEOGRAPHICLIB_HAS_STATIC_ASSERT) +# if __cplusplus >= 201103 || defined(__GXX_EXPERIMENTAL_CXX0X__) +# define GEOGRAPHICLIB_HAS_STATIC_ASSERT 1 +# elif defined(_MSC_VER) && _MSC_VER >= 1600 +// For reference, here is a table of Visual Studio and _MSC_VER +// correspondences: +// +// _MSC_VER Visual Studio +// 1100 vc5 +// 1200 vc6 +// 1300 vc7 +// 1310 vc7.1 (2003) +// 1400 vc8 (2005) +// 1500 vc9 (2008) +// 1600 vc10 (2010) +// 1700 vc11 (2012) +// 1800 vc12 (2013) First version of VS to include enough C++11 support +// 1900 vc14 (2015) +// 191[0-9] vc15 (2017) +// 192[0-9] vc16 (2019) +# define GEOGRAPHICLIB_HAS_STATIC_ASSERT 1 +# else +# define GEOGRAPHICLIB_HAS_STATIC_ASSERT 0 +# endif +#endif + +/** + * @relates GeographicLib::Constants + * A compile-time assert. Use C++11 static_assert, if available. + **********************************************************************/ +#if !defined(GEOGRAPHICLIB_STATIC_ASSERT) +# if GEOGRAPHICLIB_HAS_STATIC_ASSERT +# define GEOGRAPHICLIB_STATIC_ASSERT static_assert +# else +# define GEOGRAPHICLIB_STATIC_ASSERT(cond,reason) \ + { enum{ GEOGRAPHICLIB_STATIC_ASSERT_ENUM = 1/int(cond) }; } +# endif +#endif + +#if defined(_MSC_VER) && defined(GEOGRAPHICLIB_SHARED_LIB) && \ + GEOGRAPHICLIB_SHARED_LIB +# if GEOGRAPHICLIB_SHARED_LIB > 1 +# error GEOGRAPHICLIB_SHARED_LIB must be 0 or 1 +# elif defined(GeographicLib_SHARED_EXPORTS) +# define GEOGRAPHICLIB_EXPORT __declspec(dllexport) +# else +# define GEOGRAPHICLIB_EXPORT __declspec(dllimport) +# endif +#else +# define GEOGRAPHICLIB_EXPORT +#endif + +// Use GEOGRAPHICLIB_DEPRECATED to mark functions, types or variables as +// deprecated. Code inspired by Apache Subversion's svn_types.h file (via +// MPFR). +#if defined(__GNUC__) +# if __GNUC__ > 4 +# define GEOGRAPHICLIB_DEPRECATED(msg) __attribute__((deprecated(msg))) +# else +# define GEOGRAPHICLIB_DEPRECATED(msg) __attribute__((deprecated)) +# endif +#elif defined(_MSC_VER) && _MSC_VER >= 1300 +# define GEOGRAPHICLIB_DEPRECATED(msg) __declspec(deprecated(msg)) +#else +# define GEOGRAPHICLIB_DEPRECATED(msg) +#endif + +#include +#include + +/** + * \brief Namespace for %GeographicLib + * + * All of %GeographicLib is defined within the GeographicLib namespace. In + * addition all the header files are included via %GeographicLib/Class.hpp. + * This minimizes the likelihood of conflicts with other packages. + **********************************************************************/ +namespace GeographicLib { + + /** + * \brief %Constants needed by %GeographicLib + * + * Define constants specifying the WGS84 ellipsoid, the UTM and UPS + * projections, and various unit conversions. + * + * Example of use: + * \include example-Constants.cpp + **********************************************************************/ + class GEOGRAPHICLIB_EXPORT Constants { + private: + Constants(); // Disable constructor + + public: + /** + * A synonym for Math::degree(). + **********************************************************************/ + static real degree() { return Math::degree(); } + /** + * @return the number of radians in an arcminute. + **********************************************************************/ + static real arcminute() + { return Math::degree() / 60; } + /** + * @return the number of radians in an arcsecond. + **********************************************************************/ + static real arcsecond() + { return Math::degree() / 3600; } + + /** \name Ellipsoid parameters + **********************************************************************/ + ///@{ + /** + * @tparam T the type of the returned value. + * @return the equatorial radius of WGS84 ellipsoid (6378137 m). + **********************************************************************/ + template static T WGS84_a() + { return 6378137 * meter(); } + /** + * A synonym for WGS84_a(). + **********************************************************************/ + static real WGS84_a() { return WGS84_a(); } + /** + * @tparam T the type of the returned value. + * @return the flattening of WGS84 ellipsoid (1/298.257223563). + **********************************************************************/ + template static T WGS84_f() { + // Evaluating this as 1000000000 / T(298257223563LL) reduces the + // round-off error by about 10%. However, expressing the flattening as + // 1/298.257223563 is well ingrained. + return 1 / ( T(298257223563LL) / 1000000000 ); + } + /** + * A synonym for WGS84_f(). + **********************************************************************/ + static real WGS84_f() { return WGS84_f(); } + /** + * @tparam T the type of the returned value. + * @return the gravitational constant of the WGS84 ellipsoid, \e GM, in + * m3 s−2. + **********************************************************************/ + template static T WGS84_GM() + { return T(3986004) * 100000000 + 41800000; } + /** + * A synonym for WGS84_GM(). + **********************************************************************/ + static real WGS84_GM() { return WGS84_GM(); } + /** + * @tparam T the type of the returned value. + * @return the angular velocity of the WGS84 ellipsoid, ω, in rad + * s−1. + **********************************************************************/ + template static T WGS84_omega() + { return 7292115 / (T(1000000) * 100000); } + /** + * A synonym for WGS84_omega(). + **********************************************************************/ + static real WGS84_omega() { return WGS84_omega(); } + /** + * @tparam T the type of the returned value. + * @return the equatorial radius of GRS80 ellipsoid, \e a, in m. + **********************************************************************/ + template static T GRS80_a() + { return 6378137 * meter(); } + /** + * A synonym for GRS80_a(). + **********************************************************************/ + static real GRS80_a() { return GRS80_a(); } + /** + * @tparam T the type of the returned value. + * @return the gravitational constant of the GRS80 ellipsoid, \e GM, in + * m3 s−2. + **********************************************************************/ + template static T GRS80_GM() + { return T(3986005) * 100000000; } + /** + * A synonym for GRS80_GM(). + **********************************************************************/ + static real GRS80_GM() { return GRS80_GM(); } + /** + * @tparam T the type of the returned value. + * @return the angular velocity of the GRS80 ellipsoid, ω, in rad + * s−1. + * + * This is about 2 π 366.25 / (365.25 × 24 × 3600) rad + * s−1. 365.25 is the number of days in a Julian year and + * 365.35/366.25 converts from solar days to sidereal days. Using the + * number of days in a Gregorian year (365.2425) results in a worse + * approximation (because the Gregorian year includes the precession of the + * earth's axis). + **********************************************************************/ + template static T GRS80_omega() + { return 7292115 / (T(1000000) * 100000); } + /** + * A synonym for GRS80_omega(). + **********************************************************************/ + static real GRS80_omega() { return GRS80_omega(); } + /** + * @tparam T the type of the returned value. + * @return the dynamical form factor of the GRS80 ellipsoid, + * J2. + **********************************************************************/ + template static T GRS80_J2() + { return T(108263) / 100000000; } + /** + * A synonym for GRS80_J2(). + **********************************************************************/ + static real GRS80_J2() { return GRS80_J2(); } + /** + * @tparam T the type of the returned value. + * @return the central scale factor for UTM (0.9996). + **********************************************************************/ + template static T UTM_k0() + {return T(9996) / 10000; } + /** + * A synonym for UTM_k0(). + **********************************************************************/ + static real UTM_k0() { return UTM_k0(); } + /** + * @tparam T the type of the returned value. + * @return the central scale factor for UPS (0.994). + **********************************************************************/ + template static T UPS_k0() + { return T(994) / 1000; } + /** + * A synonym for UPS_k0(). + **********************************************************************/ + static real UPS_k0() { return UPS_k0(); } + ///@} + + /** \name SI units + **********************************************************************/ + ///@{ + /** + * @tparam T the type of the returned value. + * @return the number of meters in a meter. + * + * This is unity, but this lets the internal system of units be changed if + * necessary. + **********************************************************************/ + template static T meter() { return T(1); } + /** + * A synonym for meter(). + **********************************************************************/ + static real meter() { return meter(); } + /** + * @return the number of meters in a kilometer. + **********************************************************************/ + static real kilometer() + { return 1000 * meter(); } + /** + * @return the number of meters in a nautical mile (approximately 1 arc + * minute) + **********************************************************************/ + static real nauticalmile() + { return 1852 * meter(); } + + /** + * @tparam T the type of the returned value. + * @return the number of square meters in a square meter. + * + * This is unity, but this lets the internal system of units be changed if + * necessary. + **********************************************************************/ + template static T square_meter() + { return meter() * meter(); } + /** + * A synonym for square_meter(). + **********************************************************************/ + static real square_meter() + { return square_meter(); } + /** + * @return the number of square meters in a hectare. + **********************************************************************/ + static real hectare() + { return 10000 * square_meter(); } + /** + * @return the number of square meters in a square kilometer. + **********************************************************************/ + static real square_kilometer() + { return kilometer() * kilometer(); } + /** + * @return the number of square meters in a square nautical mile. + **********************************************************************/ + static real square_nauticalmile() + { return nauticalmile() * nauticalmile(); } + ///@} + + /** \name Anachronistic British units + **********************************************************************/ + ///@{ + /** + * @return the number of meters in an international foot. + **********************************************************************/ + static real foot() + { return real(254 * 12) / 10000 * meter(); } + /** + * @return the number of meters in a yard. + **********************************************************************/ + static real yard() { return 3 * foot(); } + /** + * @return the number of meters in a fathom. + **********************************************************************/ + static real fathom() { return 2 * yard(); } + /** + * @return the number of meters in a chain. + **********************************************************************/ + static real chain() { return 22 * yard(); } + /** + * @return the number of meters in a furlong. + **********************************************************************/ + static real furlong() { return 10 * chain(); } + /** + * @return the number of meters in a statute mile. + **********************************************************************/ + static real mile() { return 8 * furlong(); } + /** + * @return the number of square meters in an acre. + **********************************************************************/ + static real acre() { return chain() * furlong(); } + /** + * @return the number of square meters in a square statute mile. + **********************************************************************/ + static real square_mile() { return mile() * mile(); } + ///@} + + /** \name Anachronistic US units + **********************************************************************/ + ///@{ + /** + * @return the number of meters in a US survey foot. + **********************************************************************/ + static real surveyfoot() + { return real(1200) / 3937 * meter(); } + ///@} + }; + + /** + * \brief Exception handling for %GeographicLib + * + * A class to handle exceptions. It's derived from std::runtime_error so it + * can be caught by the usual catch clauses. + * + * Example of use: + * \include example-GeographicErr.cpp + **********************************************************************/ + class GeographicErr : public std::runtime_error { + public: + + /** + * Constructor + * + * @param[in] msg a string message, which is accessible in the catch + * clause via what(). + **********************************************************************/ + GeographicErr(const std::string& msg) : std::runtime_error(msg) {} + }; + +} // namespace GeographicLib + +#endif // GEOGRAPHICLIB_CONSTANTS_HPP diff --git a/src/Geo/MGRS.cpp b/src/Geo/MGRS.cpp new file mode 100644 index 000000000..da65f6639 --- /dev/null +++ b/src/Geo/MGRS.cpp @@ -0,0 +1,463 @@ +/** + * \file MGRS.cpp + * \brief Implementation for GeographicLib::MGRS class + * + * Copyright (c) Charles Karney (2008-2017) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#include "MGRS.hpp" +#include "Utility.hpp" + +namespace GeographicLib { + + using namespace std; + + const char* const MGRS::hemispheres_ = "SN"; + const char* const MGRS::utmcols_[] = { "ABCDEFGH", "JKLMNPQR", "STUVWXYZ" }; + const char* const MGRS::utmrow_ = "ABCDEFGHJKLMNPQRSTUV"; + const char* const MGRS::upscols_[] = + { "JKLPQRSTUXYZ", "ABCFGHJKLPQR", "RSTUXYZ", "ABCFGHJ" }; + const char* const MGRS::upsrows_[] = + { "ABCDEFGHJKLMNPQRSTUVWXYZ", "ABCDEFGHJKLMNP" }; + const char* const MGRS::latband_ = "CDEFGHJKLMNPQRSTUVWX"; + const char* const MGRS::upsband_ = "ABYZ"; + const char* const MGRS::digits_ = "0123456789"; + + const int MGRS::mineasting_[] = + { minupsSind_, minupsNind_, minutmcol_, minutmcol_ }; + const int MGRS::maxeasting_[] = + { maxupsSind_, maxupsNind_, maxutmcol_, maxutmcol_ }; + const int MGRS::minnorthing_[] = + { minupsSind_, minupsNind_, + minutmSrow_, minutmSrow_ - (maxutmSrow_ - minutmNrow_) }; + const int MGRS::maxnorthing_[] = + { maxupsSind_, maxupsNind_, + maxutmNrow_ + (maxutmSrow_ - minutmNrow_), maxutmNrow_ }; + + void MGRS::Forward(int zone, bool northp, real x, real y, real lat, + int prec, std::string& mgrs) { + // The smallest angle s.t., 90 - angeps() < 90 (approx 50e-12 arcsec) + // 7 = ceil(log_2(90)) + static const real angeps = ldexp(real(1), -(Math::digits() - 7)); + if (zone == UTMUPS::INVALID || + Math::isnan(x) || Math::isnan(y) || Math::isnan(lat)) { + mgrs = "INVALID"; + return; + } + bool utmp = zone != 0; + CheckCoords(utmp, northp, x, y); + if (!(zone >= UTMUPS::MINZONE && zone <= UTMUPS::MAXZONE)) + throw GeographicErr("Zone " + Utility::str(zone) + " not in [0,60]"); + if (!(prec >= -1 && prec <= maxprec_)) + throw GeographicErr("MGRS precision " + Utility::str(prec) + + " not in [-1, " + + Utility::str(int(maxprec_)) + "]"); + // Fixed char array for accumulating string. Allow space for zone, 3 block + // letters, easting + northing. Don't need to allow for terminating null. + char mgrs1[2 + 3 + 2 * maxprec_]; + int + zone1 = zone - 1, + z = utmp ? 2 : 0, + mlen = z + 3 + 2 * prec; + if (utmp) { + mgrs1[0] = digits_[ zone / base_ ]; + mgrs1[1] = digits_[ zone % base_ ]; + // This isn't necessary...! Keep y non-neg + // if (!northp) y -= maxutmSrow_ * tile_; + } + // The C++ standard mandates 64 bits for long long. But + // check, to make sure. + GEOGRAPHICLIB_STATIC_ASSERT(numeric_limits::digits >= 44, + "long long not wide enough to store 10e12"); + long long + ix = (long long)(floor(x * mult_)), + iy = (long long)(floor(y * mult_)), + m = (long long)(mult_) * (long long)(tile_); + int xh = int(ix / m), yh = int(iy / m); + if (utmp) { + int + // Correct fuzziness in latitude near equator + iband = abs(lat) > angeps ? LatitudeBand(lat) : (northp ? 0 : -1), + icol = xh - minutmcol_, + irow = UTMRow(iband, icol, yh % utmrowperiod_); + if (irow != yh - (northp ? minutmNrow_ : maxutmSrow_)) + throw GeographicErr("Latitude " + Utility::str(lat) + + " is inconsistent with UTM coordinates"); + mgrs1[z++] = latband_[10 + iband]; + mgrs1[z++] = utmcols_[zone1 % 3][icol]; + mgrs1[z++] = utmrow_[(yh + (zone1 & 1 ? utmevenrowshift_ : 0)) + % utmrowperiod_]; + } else { + bool eastp = xh >= upseasting_; + int iband = (northp ? 2 : 0) + (eastp ? 1 : 0); + mgrs1[z++] = upsband_[iband]; + mgrs1[z++] = upscols_[iband][xh - (eastp ? upseasting_ : + (northp ? minupsNind_ : + minupsSind_))]; + mgrs1[z++] = upsrows_[northp][yh - (northp ? minupsNind_ : minupsSind_)]; + } + if (prec > 0) { + ix -= m * xh; iy -= m * yh; + long long d = (long long)(pow(real(base_), maxprec_ - prec)); + ix /= d; iy /= d; + for (int c = prec; c--;) { + mgrs1[z + c ] = digits_[ix % base_]; ix /= base_; + mgrs1[z + c + prec] = digits_[iy % base_]; iy /= base_; + } + } + mgrs.resize(mlen); + copy(mgrs1, mgrs1 + mlen, mgrs.begin()); + } + + void MGRS::Forward(int zone, bool northp, real x, real y, + int prec, std::string& mgrs) { + real lat, lon; + if (zone > 0) { + // Does a rough estimate for latitude determine the latitude band? + real ys = northp ? y : y - utmNshift_; + // A cheap calculation of the latitude which results in an "allowed" + // latitude band would be + // lat = ApproxLatitudeBand(ys) * 8 + 4; + // + // Here we do a more careful job using the band letter corresponding to + // the actual latitude. + ys /= tile_; + if (abs(ys) < 1) + lat = real(0.9) * ys; // accurate enough estimate near equator + else { + real + // The poleward bound is a fit from above of lat(x,y) + // for x = 500km and y = [0km, 950km] + latp = real(0.901) * ys + (ys > 0 ? 1 : -1) * real(0.135), + // The equatorward bound is a fit from below of lat(x,y) + // for x = 900km and y = [0km, 950km] + late = real(0.902) * ys * (1 - real(1.85e-6) * ys * ys); + if (LatitudeBand(latp) == LatitudeBand(late)) + lat = latp; + else + // bounds straddle a band boundary so need to compute lat accurately + UTMUPS::Reverse(zone, northp, x, y, lat, lon); + } + } else + // Latitude isn't needed for UPS specs or for INVALID + lat = 0; + Forward(zone, northp, x, y, lat, prec, mgrs); + } + + void MGRS::Reverse(const std::string& mgrs, + int& zone, bool& northp, real& x, real& y, + int& prec, bool centerp) { + int + p = 0, + len = int(mgrs.length()); + if (len >= 3 && + toupper(mgrs[0]) == 'I' && + toupper(mgrs[1]) == 'N' && + toupper(mgrs[2]) == 'V') { + zone = UTMUPS::INVALID; + northp = false; + x = y = Math::NaN(); + prec = -2; + return; + } + int zone1 = 0; + while (p < len) { + int i = Utility::lookup(digits_, mgrs[p]); + if (i < 0) + break; + zone1 = 10 * zone1 + i; + ++p; + } + if (p > 0 && !(zone1 >= UTMUPS::MINUTMZONE && zone1 <= UTMUPS::MAXUTMZONE)) + throw GeographicErr("Zone " + Utility::str(zone1) + " not in [1,60]"); + if (p > 2) + throw GeographicErr("More than 2 digits at start of MGRS " + + mgrs.substr(0, p)); + if (len - p < 1) + throw GeographicErr("MGRS string too short " + mgrs); + bool utmp = zone1 != UTMUPS::UPS; + int zonem1 = zone1 - 1; + const char* band = utmp ? latband_ : upsband_; + int iband = Utility::lookup(band, mgrs[p++]); + if (iband < 0) + throw GeographicErr("Band letter " + Utility::str(mgrs[p-1]) + " not in " + + (utmp ? "UTM" : "UPS") + " set " + band); + bool northp1 = iband >= (utmp ? 10 : 2); + if (p == len) { // Grid zone only (ignore centerp) + // Approx length of a degree of meridian arc in units of tile. + real deg = real(utmNshift_) / (90 * tile_); + zone = zone1; + northp = northp1; + if (utmp) { + // Pick central meridian except for 31V + x = ((zone == 31 && iband == 17) ? 4 : 5) * tile_; + // Pick center of 8deg latitude bands + y = floor(8 * (iband - real(9.5)) * deg + real(0.5)) * tile_ + + (northp ? 0 : utmNshift_); + } else { + // Pick point at lat 86N or 86S + x = ((iband & 1 ? 1 : -1) * floor(4 * deg + real(0.5)) + + upseasting_) * tile_; + // Pick point at lon 90E or 90W. + y = upseasting_ * tile_; + } + prec = -1; + return; + } else if (len - p < 2) + throw GeographicErr("Missing row letter in " + mgrs); + const char* col = utmp ? utmcols_[zonem1 % 3] : upscols_[iband]; + const char* row = utmp ? utmrow_ : upsrows_[northp1]; + int icol = Utility::lookup(col, mgrs[p++]); + if (icol < 0) + throw GeographicErr("Column letter " + Utility::str(mgrs[p-1]) + + " not in " + + (utmp ? "zone " + mgrs.substr(0, p-2) : + "UPS band " + Utility::str(mgrs[p-2])) + + " set " + col ); + int irow = Utility::lookup(row, mgrs[p++]); + if (irow < 0) + throw GeographicErr("Row letter " + Utility::str(mgrs[p-1]) + " not in " + + (utmp ? "UTM" : + "UPS " + Utility::str(hemispheres_[northp1])) + + " set " + row); + if (utmp) { + if (zonem1 & 1) + irow = (irow + utmrowperiod_ - utmevenrowshift_) % utmrowperiod_; + iband -= 10; + irow = UTMRow(iband, icol, irow); + if (irow == maxutmSrow_) + throw GeographicErr("Block " + mgrs.substr(p-2, 2) + + " not in zone/band " + mgrs.substr(0, p-2)); + + irow = northp1 ? irow : irow + 100; + icol = icol + minutmcol_; + } else { + bool eastp = iband & 1; + icol += eastp ? upseasting_ : (northp1 ? minupsNind_ : minupsSind_); + irow += northp1 ? minupsNind_ : minupsSind_; + } + int prec1 = (len - p)/2; + real + unit = 1, + x1 = icol, + y1 = irow; + for (int i = 0; i < prec1; ++i) { + unit *= base_; + int + ix = Utility::lookup(digits_, mgrs[p + i]), + iy = Utility::lookup(digits_, mgrs[p + i + prec1]); + if (ix < 0 || iy < 0) + throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p)); + x1 = base_ * x1 + ix; + y1 = base_ * y1 + iy; + } + if ((len - p) % 2) { + if (Utility::lookup(digits_, mgrs[len - 1]) < 0) + throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p)); + else + throw GeographicErr("Not an even number of digits in " + + mgrs.substr(p)); + } + if (prec1 > maxprec_) + throw GeographicErr("More than " + Utility::str(2*maxprec_) + + " digits in " + mgrs.substr(p)); + if (centerp) { + unit *= 2; x1 = 2 * x1 + 1; y1 = 2 * y1 + 1; + } + zone = zone1; + northp = northp1; + x = (tile_ * x1) / unit; + y = (tile_ * y1) / unit; + prec = prec1; + } + + void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) { + // Limits are all multiples of 100km and are all closed on the lower end + // and open on the upper end -- and this is reflected in the error + // messages. However if a coordinate lies on the excluded upper end (e.g., + // after rounding), it is shifted down by eps. This also folds UTM + // northings to the correct N/S hemisphere. + + // The smallest length s.t., 1.0e7 - eps() < 1.0e7 (approx 1.9 nm) + // 25 = ceil(log_2(2e7)) -- use half circumference here because + // northing 195e5 is a legal in the "southern" hemisphere. + static const real eps = ldexp(real(1), -(Math::digits() - 25)); + int + ix = int(floor(x / tile_)), + iy = int(floor(y / tile_)), + ind = (utmp ? 2 : 0) + (northp ? 1 : 0); + if (! (ix >= mineasting_[ind] && ix < maxeasting_[ind]) ) { + if (ix == maxeasting_[ind] && x == maxeasting_[ind] * tile_) + x -= eps; + else + throw GeographicErr("Easting " + Utility::str(int(floor(x/1000))) + + "km not in MGRS/" + + (utmp ? "UTM" : "UPS") + " range for " + + (northp ? "N" : "S" ) + " hemisphere [" + + Utility::str(mineasting_[ind]*tile_/1000) + + "km, " + + Utility::str(maxeasting_[ind]*tile_/1000) + + "km)"); + } + if (! (iy >= minnorthing_[ind] && iy < maxnorthing_[ind]) ) { + if (iy == maxnorthing_[ind] && y == maxnorthing_[ind] * tile_) + y -= eps; + else + throw GeographicErr("Northing " + Utility::str(int(floor(y/1000))) + + "km not in MGRS/" + + (utmp ? "UTM" : "UPS") + " range for " + + (northp ? "N" : "S" ) + " hemisphere [" + + Utility::str(minnorthing_[ind]*tile_/1000) + + "km, " + + Utility::str(maxnorthing_[ind]*tile_/1000) + + "km)"); + } + + // Correct the UTM northing and hemisphere if necessary + if (utmp) { + if (northp && iy < minutmNrow_) { + northp = false; + y += utmNshift_; + } else if (!northp && iy >= maxutmSrow_) { + if (y == maxutmSrow_ * tile_) + // If on equator retain S hemisphere + y -= eps; + else { + northp = true; + y -= utmNshift_; + } + } + } + } + + int MGRS::UTMRow(int iband, int icol, int irow) { + // Input is iband = band index in [-10, 10) (as returned by LatitudeBand), + // icol = column index in [0,8) with origin of easting = 100km, and irow = + // periodic row index in [0,20) with origin = equator. Output is true row + // index in [-90, 95). Returns maxutmSrow_ = 100, if irow and iband are + // incompatible. + + // Estimate center row number for latitude band + // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles + real c = 100 * (8 * iband + 4)/real(90); + bool northp = iband >= 0; + // These are safe bounds on the rows + // iband minrow maxrow + // -10 -90 -81 + // -9 -80 -72 + // -8 -71 -63 + // -7 -63 -54 + // -6 -54 -45 + // -5 -45 -36 + // -4 -36 -27 + // -3 -27 -18 + // -2 -18 -9 + // -1 -9 -1 + // 0 0 8 + // 1 8 17 + // 2 17 26 + // 3 26 35 + // 4 35 44 + // 5 44 53 + // 6 53 62 + // 7 62 70 + // 8 71 79 + // 9 80 94 + int + minrow = iband > -10 ? + int(floor(c - real(4.3) - real(0.1) * northp)) : -90, + maxrow = iband < 9 ? + int(floor(c + real(4.4) - real(0.1) * northp)) : 94, + baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2; + // Offset irow by the multiple of utmrowperiod_ which brings it as close as + // possible to the center of the latitude band, (minrow + maxrow) / 2. + // (Add maxutmSrow_ = 5 * utmrowperiod_ to ensure operand is positive.) + irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow; + if (!( irow >= minrow && irow <= maxrow )) { + // Outside the safe bounds, so need to check... + // Northing = 71e5 and 80e5 intersect band boundaries + // y = 71e5 in scol = 2 (x = [3e5,4e5] and x = [6e5,7e5]) + // y = 80e5 in scol = 1 (x = [2e5,3e5] and x = [7e5,8e5]) + // This holds for all the ellipsoids given in NGA.SIG.0012_2.0.0_UTMUPS. + // The following deals with these special cases. + int + // Fold [-10,-1] -> [9,0] + sband = iband >= 0 ? iband : -iband - 1, + // Fold [-90,-1] -> [89,0] + srow = irow >= 0 ? irow : -irow - 1, + // Fold [4,7] -> [3,0] + scol = icol < 4 ? icol : -icol + 7; + // For example, the safe rows for band 8 are 71 - 79. However row 70 is + // allowed if scol = [2,3] and row 80 is allowed if scol = [0,1]. + if ( ! ( (srow == 70 && sband == 8 && scol >= 2) || + (srow == 71 && sband == 7 && scol <= 2) || + (srow == 79 && sband == 9 && scol >= 1) || + (srow == 80 && sband == 8 && scol <= 1) ) ) + irow = maxutmSrow_; + } + return irow; + } + + void MGRS::Check() { + real lat, lon, x, y, t = tile_; int zone; bool northp; + UTMUPS::Reverse(31, true , 1*t, 0*t, lat, lon); + if (!( lon < 0 )) + throw GeographicErr("MGRS::Check: equator coverage failure"); + UTMUPS::Reverse(31, true , 1*t, 95*t, lat, lon); + if (!( lat > 84 )) + throw GeographicErr("MGRS::Check: UTM doesn't reach latitude = 84"); + UTMUPS::Reverse(31, false, 1*t, 10*t, lat, lon); + if (!( lat < -80 )) + throw GeographicErr("MGRS::Check: UTM doesn't reach latitude = -80"); + UTMUPS::Forward(56, 3, zone, northp, x, y, 32); + if (!( x > 1*t )) + throw GeographicErr("MGRS::Check: Norway exception creates a gap"); + UTMUPS::Forward(72, 21, zone, northp, x, y, 35); + if (!( x > 1*t )) + throw GeographicErr("MGRS::Check: Svalbard exception creates a gap"); + UTMUPS::Reverse(0, true , 20*t, 13*t, lat, lon); + if (!( lat < 84 )) + throw + GeographicErr("MGRS::Check: North UPS doesn't reach latitude = 84"); + UTMUPS::Reverse(0, false, 20*t, 8*t, lat, lon); + if (!( lat > -80 )) + throw + GeographicErr("MGRS::Check: South UPS doesn't reach latitude = -80"); + // Entries are [band, x, y] either side of the band boundaries. Units for + // x, y are t = 100km. + const short tab[] = { + 0, 5, 0, 0, 9, 0, // south edge of band 0 + 0, 5, 8, 0, 9, 8, // north edge of band 0 + 1, 5, 9, 1, 9, 9, // south edge of band 1 + 1, 5, 17, 1, 9, 17, // north edge of band 1 + 2, 5, 18, 2, 9, 18, // etc. + 2, 5, 26, 2, 9, 26, + 3, 5, 27, 3, 9, 27, + 3, 5, 35, 3, 9, 35, + 4, 5, 36, 4, 9, 36, + 4, 5, 44, 4, 9, 44, + 5, 5, 45, 5, 9, 45, + 5, 5, 53, 5, 9, 53, + 6, 5, 54, 6, 9, 54, + 6, 5, 62, 6, 9, 62, + 7, 5, 63, 7, 9, 63, + 7, 5, 70, 7, 7, 70, 7, 7, 71, 7, 9, 71, // y = 71t crosses boundary + 8, 5, 71, 8, 6, 71, 8, 6, 72, 8, 9, 72, // between bands 7 and 8. + 8, 5, 79, 8, 8, 79, 8, 8, 80, 8, 9, 80, // y = 80t crosses boundary + 9, 5, 80, 9, 7, 80, 9, 7, 81, 9, 9, 81, // between bands 8 and 9. + 9, 5, 95, 9, 9, 95, // north edge of band 9 + }; + const int bandchecks = sizeof(tab) / (3 * sizeof(short)); + for (int i = 0; i < bandchecks; ++i) { + UTMUPS::Reverse(38, true, tab[3*i+1]*t, tab[3*i+2]*t, lat, lon); + if (!( LatitudeBand(lat) == tab[3*i+0] )) + throw GeographicErr("MGRS::Check: Band error, b = " + + Utility::str(tab[3*i+0]) + ", x = " + + Utility::str(tab[3*i+1]) + "00km, y = " + + Utility::str(tab[3*i+2]) + "00km"); + } + } + +} // namespace GeographicLib diff --git a/src/Geo/MGRS.hpp b/src/Geo/MGRS.hpp new file mode 100644 index 000000000..697ab72e0 --- /dev/null +++ b/src/Geo/MGRS.hpp @@ -0,0 +1,361 @@ +/** + * \file MGRS.hpp + * \brief Header for GeographicLib::MGRS class + * + * Copyright (c) Charles Karney (2008-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#if !defined(GEOGRAPHICLIB_MGRS_HPP) +#define GEOGRAPHICLIB_MGRS_HPP 1 + +#include "Constants.hpp" +#include "UTMUPS.hpp" + +#if defined(_MSC_VER) +// Squelch warnings about dll vs string +# pragma warning (push) +# pragma warning (disable: 4251) +#endif + +namespace GeographicLib { + + /** + * \brief Convert between UTM/UPS and %MGRS + * + * MGRS is defined in Chapter 3 of + * - J. W. Hager, L. L. Fry, S. S. Jacks, D. R. Hill, + * + * Datums, Ellipsoids, Grids, and Grid Reference Systems, + * Defense Mapping Agency, Technical Manual TM8358.1 (1990). + * . + * This document has been updated by the two NGA documents + * - + * Universal Grids and Grid Reference Systems, + * NGA.STND.0037_2.0.0_GRIDS (2014). + * - + * The Universal Grids and the Transverse Mercator and Polar Stereographic + * Map Projections, NGA.SIG.0012_2.0.0_UTMUPS (2014). + * + * This implementation has the following properties: + * - The conversions are closed, i.e., output from Forward is legal input for + * Reverse and vice versa. Conversion in both directions preserve the + * UTM/UPS selection and the UTM zone. + * - Forward followed by Reverse and vice versa is approximately the + * identity. (This is affected in predictable ways by errors in + * determining the latitude band and by loss of precision in the MGRS + * coordinates.) + * - The trailing digits produced by Forward are consistent as the precision + * is varied. Specifically, the digits are obtained by operating on the + * easting with ⌊106 x⌋ and extracting the + * required digits from the resulting number (and similarly for the + * northing). + * - All MGRS coordinates truncate to legal 100 km blocks. All MGRS + * coordinates with a legal 100 km block prefix are legal (even though the + * latitude band letter may now belong to a neighboring band). + * - The range of UTM/UPS coordinates allowed for conversion to MGRS + * coordinates is the maximum consistent with staying within the letter + * ranges of the MGRS scheme. + * - All the transformations are implemented as static methods in the MGRS + * class. + * + * The NGA software package + * geotrans + * also provides conversions to and from MGRS. Version 3.0 (and earlier) + * suffers from some drawbacks: + * - Inconsistent rules are used to determine the whether a particular MGRS + * coordinate is legal. A more systematic approach is taken here. + * - The underlying projections are not very accurately implemented. + * + * Example of use: + * \include example-MGRS.cpp + **********************************************************************/ + class GEOGRAPHICLIB_EXPORT MGRS { + private: + typedef Math::real real; + static const char* const hemispheres_; + static const char* const utmcols_[3]; + static const char* const utmrow_; + static const char* const upscols_[4]; + static const char* const upsrows_[2]; + static const char* const latband_; + static const char* const upsband_; + static const char* const digits_; + + static const int mineasting_[4]; + static const int maxeasting_[4]; + static const int minnorthing_[4]; + static const int maxnorthing_[4]; + enum { + base_ = 10, + // Top-level tiles are 10^5 m = 100 km on a side + tilelevel_ = 5, + // Period of UTM row letters + utmrowperiod_ = 20, + // Row letters are shifted by 5 for even zones + utmevenrowshift_ = 5, + // Maximum precision is um + maxprec_ = 5 + 6, + // For generating digits at maxprec + mult_ = 1000000, + }; + static void CheckCoords(bool utmp, bool& northp, real& x, real& y); + static int UTMRow(int iband, int icol, int irow); + + friend class UTMUPS; // UTMUPS::StandardZone calls LatitudeBand + // Return latitude band number [-10, 10) for the given latitude (degrees). + // The bands are reckoned in include their southern edges. + static int LatitudeBand(real lat) { + using std::floor; + int ilat = int(floor(lat)); + return (std::max)(-10, (std::min)(9, (ilat + 80)/8 - 10)); + } + // Return approximate latitude band number [-10, 10) for the given northing + // (meters). With this rule, each 100km tile would have a unique band + // letter corresponding to the latitude at the center of the tile. This + // function isn't currently used. + static int ApproxLatitudeBand(real y) { + // northing at tile center in units of tile = 100km + using std::floor; using std::abs; + real ya = floor( (std::min)(real(88), abs(y/tile_)) ) + + real(0.5); + // convert to lat (mult by 90/100) and then to band (divide by 8) + // the +1 fine tunes the boundary between bands 3 and 4 + int b = int(floor( ((ya * 9 + 1) / 10) / 8 )); + // For the northern hemisphere we have + // band rows num + // N 0 0:8 9 + // P 1 9:17 9 + // Q 2 18:26 9 + // R 3 27:34 8 + // S 4 35:43 9 + // T 5 44:52 9 + // U 6 53:61 9 + // V 7 62:70 9 + // W 8 71:79 9 + // X 9 80:94 15 + return y >= 0 ? b : -(b + 1); + } + // UTMUPS access these enums + enum { + tile_ = 100000, // Size MGRS blocks + minutmcol_ = 1, + maxutmcol_ = 9, + minutmSrow_ = 10, + maxutmSrow_ = 100, // Also used for UTM S false northing + minutmNrow_ = 0, // Also used for UTM N false northing + maxutmNrow_ = 95, + minupsSind_ = 8, // These 4 ind's apply to easting and northing + maxupsSind_ = 32, + minupsNind_ = 13, + maxupsNind_ = 27, + upseasting_ = 20, // Also used for UPS false northing + utmeasting_ = 5, // UTM false easting + // Difference between S hemisphere northing and N hemisphere northing + utmNshift_ = (maxutmSrow_ - minutmNrow_) * tile_ + }; + MGRS(); // Disable constructor + + public: + + /** + * Convert UTM or UPS coordinate to an MGRS coordinate. + * + * @param[in] zone UTM zone (zero means UPS). + * @param[in] northp hemisphere (true means north, false means south). + * @param[in] x easting of point (meters). + * @param[in] y northing of point (meters). + * @param[in] prec precision relative to 100 km. + * @param[out] mgrs MGRS string. + * @exception GeographicErr if \e zone, \e x, or \e y is outside its + * allowed range. + * @exception GeographicErr if the memory for the MGRS string can't be + * allocated. + * + * \e prec specifies the precision of the MGRS string as follows: + * - \e prec = −1 (min), only the grid zone is returned + * - \e prec = 0, 100 km + * - \e prec = 1, 10 km + * - \e prec = 2, 1 km + * - \e prec = 3, 100 m + * - \e prec = 4, 10 m + * - \e prec = 5, 1 m + * - \e prec = 6, 0.1 m + * - … + * - \e prec = 11 (max), 1 μm + * + * UTM eastings are allowed to be in the range [100 km, 900 km], northings + * are allowed to be in in [0 km, 9500 km] for the northern hemisphere and + * in [1000 km, 10000 km] for the southern hemisphere. (However UTM + * northings can be continued across the equator. So the actual limits on + * the northings are [−9000 km, 9500 km] for the "northern" + * hemisphere and [1000 km, 19500 km] for the "southern" hemisphere.) + * + * UPS eastings/northings are allowed to be in the range [1300 km, 2700 km] + * in the northern hemisphere and in [800 km, 3200 km] in the southern + * hemisphere. + * + * The ranges are 100 km more restrictive than for the conversion between + * geographic coordinates and UTM and UPS given by UTMUPS. These + * restrictions are dictated by the allowed letters in MGRS coordinates. + * The choice of 9500 km for the maximum northing for northern hemisphere + * and of 1000 km as the minimum northing for southern hemisphere provide + * at least 0.5 degree extension into standard UPS zones. The upper ends + * of the ranges for the UPS coordinates is dictated by requiring symmetry + * about the meridians 0E and 90E. + * + * All allowed UTM and UPS coordinates may now be converted to legal MGRS + * coordinates with the proviso that eastings and northings on the upper + * boundaries are silently reduced by about 4 nm (4 nanometers) to place + * them \e within the allowed range. (This includes reducing a southern + * hemisphere northing of 10000 km by 4 nm so that it is placed in latitude + * band M.) The UTM or UPS coordinates are truncated to requested + * precision to determine the MGRS coordinate. Thus in UTM zone 38n, the + * square area with easting in [444 km, 445 km) and northing in [3688 km, + * 3689 km) maps to MGRS coordinate 38SMB4488 (at \e prec = 2, 1 km), + * Khulani Sq., Baghdad. + * + * The UTM/UPS selection and the UTM zone is preserved in the conversion to + * MGRS coordinate. Thus for \e zone > 0, the MGRS coordinate begins with + * the zone number followed by one of [C--M] for the southern + * hemisphere and [N--X] for the northern hemisphere. For \e zone = + * 0, the MGRS coordinates begins with one of [AB] for the southern + * hemisphere and [XY] for the northern hemisphere. + * + * The conversion to the MGRS is exact for prec in [0, 5] except that a + * neighboring latitude band letter may be given if the point is within 5nm + * of a band boundary. For prec in [6, 11], the conversion is accurate to + * roundoff. + * + * If \e prec = −1, then the "grid zone designation", e.g., 18T, is + * returned. This consists of the UTM zone number (absent for UPS) and the + * first letter of the MGRS string which labels the latitude band for UTM + * and the hemisphere for UPS. + * + * If \e x or \e y is NaN or if \e zone is UTMUPS::INVALID, the returned + * MGRS string is "INVALID". + * + * Return the result via a reference argument to avoid the overhead of + * allocating a potentially large number of small strings. If an error is + * thrown, then \e mgrs is unchanged. + **********************************************************************/ + static void Forward(int zone, bool northp, real x, real y, + int prec, std::string& mgrs); + + /** + * Convert UTM or UPS coordinate to an MGRS coordinate when the latitude is + * known. + * + * @param[in] zone UTM zone (zero means UPS). + * @param[in] northp hemisphere (true means north, false means south). + * @param[in] x easting of point (meters). + * @param[in] y northing of point (meters). + * @param[in] lat latitude (degrees). + * @param[in] prec precision relative to 100 km. + * @param[out] mgrs MGRS string. + * @exception GeographicErr if \e zone, \e x, or \e y is outside its + * allowed range. + * @exception GeographicErr if \e lat is inconsistent with the given UTM + * coordinates. + * @exception std::bad_alloc if the memory for \e mgrs can't be allocated. + * + * The latitude is ignored for \e zone = 0 (UPS); otherwise the latitude is + * used to determine the latitude band and this is checked for consistency + * using the same tests as Reverse. + **********************************************************************/ + static void Forward(int zone, bool northp, real x, real y, real lat, + int prec, std::string& mgrs); + + /** + * Convert a MGRS coordinate to UTM or UPS coordinates. + * + * @param[in] mgrs MGRS string. + * @param[out] zone UTM zone (zero means UPS). + * @param[out] northp hemisphere (true means north, false means south). + * @param[out] x easting of point (meters). + * @param[out] y northing of point (meters). + * @param[out] prec precision relative to 100 km. + * @param[in] centerp if true (default), return center of the MGRS square, + * else return SW (lower left) corner. + * @exception GeographicErr if \e mgrs is illegal. + * + * All conversions from MGRS to UTM/UPS are permitted provided the MGRS + * coordinate is a possible result of a conversion in the other direction. + * (The leading 0 may be dropped from an input MGRS coordinate for UTM + * zones 1--9.) In addition, MGRS coordinates with a neighboring + * latitude band letter are permitted provided that some portion of the + * 100 km block is within the given latitude band. Thus + * - 38VLS and 38WLS are allowed (latitude 64N intersects the square + * 38[VW]LS); but 38VMS is not permitted (all of 38WMS is north of 64N) + * - 38MPE and 38NPF are permitted (they straddle the equator); but 38NPE + * and 38MPF are not permitted (the equator does not intersect either + * block). + * - Similarly ZAB and YZB are permitted (they straddle the prime + * meridian); but YAB and ZZB are not (the prime meridian does not + * intersect either block). + * + * The UTM/UPS selection and the UTM zone is preserved in the conversion + * from MGRS coordinate. The conversion is exact for prec in [0, 5]. With + * \e centerp = true, the conversion from MGRS to geographic and back is + * stable. This is not assured if \e centerp = false. + * + * If a "grid zone designation" (for example, 18T or A) is given, then some + * suitable (but essentially arbitrary) point within that grid zone is + * returned. The main utility of the conversion is to allow \e zone and \e + * northp to be determined. In this case, the \e centerp parameter is + * ignored and \e prec is set to −1. + * + * If the first 3 characters of \e mgrs are "INV", then \e x and \e y are + * set to NaN, \e zone is set to UTMUPS::INVALID, and \e prec is set to + * −2. + * + * If an exception is thrown, then the arguments are unchanged. + **********************************************************************/ + static void Reverse(const std::string& mgrs, + int& zone, bool& northp, real& x, real& y, + int& prec, bool centerp = true); + + /** \name Inspector functions + **********************************************************************/ + ///@{ + /** + * @return \e a the equatorial radius of the WGS84 ellipsoid (meters). + * + * (The WGS84 value is returned because the UTM and UPS projections are + * based on this ellipsoid.) + **********************************************************************/ + static Math::real EquatorialRadius() { return UTMUPS::EquatorialRadius(); } + + /** + * @return \e f the flattening of the WGS84 ellipsoid. + * + * (The WGS84 value is returned because the UTM and UPS projections are + * based on this ellipsoid.) + **********************************************************************/ + static Math::real Flattening() { return UTMUPS::Flattening(); } + + /** + * \deprecated An old name for EquatorialRadius(). + **********************************************************************/ + // GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()") + static Math::real MajorRadius() { return EquatorialRadius(); } + ///@} + + /** + * Perform some checks on the UTMUPS coordinates on this ellipsoid. Throw + * an error if any of the assumptions made in the MGRS class is not true. + * This check needs to be carried out if the ellipsoid parameters (or the + * UTM/UPS scales) are ever changed. + **********************************************************************/ + static void Check(); + + }; + +} // namespace GeographicLib + +#if defined(_MSC_VER) +# pragma warning (pop) +#endif + +#endif // GEOGRAPHICLIB_MGRS_HPP diff --git a/src/Geo/Math.cpp b/src/Geo/Math.cpp new file mode 100644 index 000000000..4b8226cec --- /dev/null +++ b/src/Geo/Math.cpp @@ -0,0 +1,597 @@ +/** + * \file Math.cpp + * \brief Implementation for GeographicLib::Math class + * + * Copyright (c) Charles Karney (2015-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#include "Math.hpp" + +#if defined(_MSC_VER) +// Squelch warnings about constant conditional expressions +# pragma warning (disable: 4127) +#endif + +namespace GeographicLib { + + using namespace std; + + void Math::dummy() { + GEOGRAPHICLIB_STATIC_ASSERT(GEOGRAPHICLIB_PRECISION >= 1 && + GEOGRAPHICLIB_PRECISION <= 5, + "Bad value of precision"); + } + + int Math::digits() { +#if GEOGRAPHICLIB_PRECISION != 5 + return std::numeric_limits::digits; +#else + return std::numeric_limits::digits(); +#endif + } + + int Math::set_digits(int ndigits) { +#if GEOGRAPHICLIB_PRECISION != 5 + (void)ndigits; +#else + mpfr::mpreal::set_default_prec(ndigits >= 2 ? ndigits : 2); +#endif + return digits(); + } + + int Math::digits10() { +#if GEOGRAPHICLIB_PRECISION != 5 + return std::numeric_limits::digits10; +#else + return std::numeric_limits::digits10(); +#endif + } + + int Math::extra_digits() { + return + digits10() > std::numeric_limits::digits10 ? + digits10() - std::numeric_limits::digits10 : 0; + } + + template T Math::hypot(T x, T y) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::hypot; return hypot(x, y); +#else + x = abs(x); y = abs(y); + if (x < y) std::swap(x, y); // Now x >= y >= 0 + y /= (x != 0 ? x : 1); + return x * sqrt(1 + y * y); + // For an alternative (square-root free) method see + // C. Moler and D. Morrision (1983) https://doi.org/10.1147/rd.276.0577 + // and A. A. Dubrulle (1983) https://doi.org/10.1147/rd.276.0582 +#endif + } + + template T Math::expm1(T x) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::expm1; return expm1(x); +#else + GEOGRAPHICLIB_VOLATILE T + y = exp(x), + z = y - 1; + // The reasoning here is similar to that for log1p. The expression + // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y - + // 1)/log(y) is a slowly varying quantity near y = 1 and is accurately + // computed. + return abs(x) > 1 ? z : (z == 0 ? x : x * z / log(y)); +#endif + } + + template T Math::log1p(T x) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::log1p; return log1p(x); +#else + GEOGRAPHICLIB_VOLATILE T + y = 1 + x, + z = y - 1; + // Here's the explanation for this magic: y = 1 + z, exactly, and z + // approx x, thus log(y)/z (which is nearly constant near z = 0) returns + // a good approximation to the true log(1 + x)/x. The multiplication x * + // (log(y)/z) introduces little additional error. + return z == 0 ? x : x * log(y) / z; +#endif + } + + template T Math::asinh(T x) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::asinh; return asinh(x); +#else + T y = abs(x); // Enforce odd parity + y = log1p(y * (1 + y/(hypot(T(1), y) + 1))); + return x > 0 ? y : (x < 0 ? -y : x); // asinh(-0.0) = -0.0 +#endif + } + + template T Math::atanh(T x) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::atanh; return atanh(x); +#else + T y = abs(x); // Enforce odd parity + y = log1p(2 * y/(1 - y))/2; + return x > 0 ? y : (x < 0 ? -y : x); // atanh(-0.0) = -0.0 +#endif + } + + template T Math::copysign(T x, T y) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::copysign; return copysign(x, y); +#else + // NaN counts as positive + return abs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1); +#endif + } + + template T Math::cbrt(T x) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::cbrt; return cbrt(x); +#else + T y = pow(abs(x), 1/T(3)); // Return the real cube root + return x > 0 ? y : (x < 0 ? -y : x); // cbrt(-0.0) = -0.0 +#endif + } + + template T Math::remainder(T x, T y) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::remainder; return remainder(x, y); +#else + y = abs(y); // The result doesn't depend on the sign of y + T z = fmod(x, y); + if (z == 0) + // This shouldn't be necessary. However, before version 14 (2015), + // Visual Studio had problems dealing with -0.0. Specifically + // VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0 + // python 2.7 on Windows 32-bit machines has the same problem. + z = copysign(z, x); + else if (2 * abs(z) == y) + z -= fmod(x, 2 * y) - z; // Implement ties to even + else if (2 * abs(z) > y) + z += (z < 0 ? y : -y); // Fold remaining cases to (-y/2, y/2) + return z; +#endif + } + + template T Math::remquo(T x, T y, int* n) { + // boost::math::remquo doesn't handle nans correctly +#if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 + using std::remquo; return remquo(x, y, n); +#else + T z = remainder(x, y); + if (n) { + T + a = remainder(x, 2 * y), + b = remainder(x, 4 * y), + c = remainder(x, 8 * y); + *n = (a > z ? 1 : (a < z ? -1 : 0)); + *n += (b > a ? 2 : (b < a ? -2 : 0)); + *n += (c > b ? 4 : (c < b ? -4 : 0)); + if (y < 0) *n *= -1; + if (y != 0) { + if (x/y > 0 && *n <= 0) + *n += 8; + else if (x/y < 0 && *n >= 0) + *n -= 8; + } + } + return z; +#endif + } + + template T Math::round(T x) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::round; return round(x); +#else + // The handling of corner cases is copied from boost; see + // https://github.com/boostorg/math/pull/8 + // with improvements to return -0 when appropriate. + if (0 < x && x < T(0.5)) + return +T(0); + else if (0 > x && x > -T(0.5)) + return -T(0); + else if (x > 0) { + T t = ceil(x); + return t - x > T(0.5) ? t - 1 : t; + } else if (x < 0) { + T t = floor(x); + return x - t > T(0.5) ? t + 1 : t; + } else // +/-0 and NaN + return x; // Retain sign of 0 +#endif + } + + template long Math::lround(T x) { +#if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION != 5 + using std::lround; return lround(x); +#else + // Default value for overflow + NaN + (x == LONG_MIN) + long r = std::numeric_limits::min(); + x = round(x); + if (abs(x) < -T(r)) // Assume T(LONG_MIN) is exact + r = long(x); + return r; +#endif + } + + template T Math::fma(T x, T y, T z) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::fma; return fma(x, y, z); +#else + return x * y + z; +#endif + } + + template T Math::sum(T u, T v, T& t) { + GEOGRAPHICLIB_VOLATILE T s = u + v; + GEOGRAPHICLIB_VOLATILE T up = s - v; + GEOGRAPHICLIB_VOLATILE T vpp = s - up; + up -= u; + vpp -= v; + t = -(up + vpp); + // u + v = s + t + // = round(u + v) + t + return s; + } + + template T Math::AngRound(T x) { + static const T z = 1/T(16); + if (x == 0) return 0; + GEOGRAPHICLIB_VOLATILE T y = abs(x); + // The compiler mustn't "simplify" z - (z - y) to y + y = y < z ? z - (z - y) : y; + return x < 0 ? -y : y; + } + + template void Math::sincosd(T x, T& sinx, T& cosx) { + // In order to minimize round-off errors, this function exactly reduces + // the argument to the range [-45, 45] before converting it to radians. + T r; int q; + // N.B. the implementation of remquo in glibc pre 2.22 were buggy. See + // https://sourceware.org/bugzilla/show_bug.cgi?id=17569 + // This was fixed in version 2.22 on 2015-08-05 + r = remquo(x, T(90), &q); // now abs(r) <= 45 + r *= degree(); + // g++ -O turns these two function calls into a call to sincos + T s = sin(r), c = cos(r); +#if defined(_MSC_VER) && _MSC_VER < 1900 + // Before version 14 (2015), Visual Studio had problems dealing + // with -0.0. Specifically + // VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0 + // VC 12 and 64-bit compile: sin(-0.0) -> +0.0 + // AngNormalize has a similar fix. + // python 2.7 on Windows 32-bit machines has the same problem. + if (x == 0) s = x; +#endif + switch (unsigned(q) & 3U) { + case 0U: sinx = s; cosx = c; break; + case 1U: sinx = c; cosx = -s; break; + case 2U: sinx = -s; cosx = -c; break; + default: sinx = -c; cosx = s; break; // case 3U + } + // Set sign of 0 results. -0 only produced for sin(-0) + if (x != 0) { sinx += T(0); cosx += T(0); } + } + + template T Math::sind(T x) { + // See sincosd + T r; int q; + r = remquo(x, T(90), &q); // now abs(r) <= 45 + r *= degree(); + unsigned p = unsigned(q); + r = p & 1U ? cos(r) : sin(r); + if (p & 2U) r = -r; + if (x != 0) r += T(0); + return r; + } + + template T Math::cosd(T x) { + // See sincosd + T r; int q; + r = remquo(x, T(90), &q); // now abs(r) <= 45 + r *= degree(); + unsigned p = unsigned(q + 1); + r = p & 1U ? cos(r) : sin(r); + if (p & 2U) r = -r; + return T(0) + r; + } + + template T Math::tand(T x) { + static const T overflow = 1 / sq(std::numeric_limits::epsilon()); + T s, c; + sincosd(x, s, c); + return c != 0 ? s / c : (s < 0 ? -overflow : overflow); + } + + template T Math::atan2d(T y, T x) { + // In order to minimize round-off errors, this function rearranges the + // arguments so that result of atan2 is in the range [-pi/4, pi/4] before + // converting it to degrees and mapping the result to the correct + // quadrant. + int q = 0; + if (abs(y) > abs(x)) { std::swap(x, y); q = 2; } + if (x < 0) { x = -x; ++q; } + // here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4] + T ang = atan2(y, x) / degree(); + switch (q) { + // Note that atan2d(-0.0, 1.0) will return -0. However, we expect that + // atan2d will not be called with y = -0. If need be, include + // + // case 0: ang = 0 + ang; break; + // + // and handle mpfr as in AngRound. + case 1: ang = (y >= 0 ? 180 : -180) - ang; break; + case 2: ang = 90 - ang; break; + case 3: ang = -90 + ang; break; + } + return ang; + } + + template T Math::atand(T x) + { return atan2d(x, T(1)); } + + template T Math::eatanhe(T x, T es) { + return es > T(0) ? es * atanh(es * x) : -es * atan(es * x); + } + + template T Math::taupf(T tau, T es) { + T tau1 = hypot(T(1), tau), + sig = sinh( eatanhe(tau / tau1, es ) ); + return hypot(T(1), sig) * tau - sig * tau1; + } + + template T Math::tauf(T taup, T es) { + const int numit = 5; + const T tol = sqrt(numeric_limits::epsilon()) / T(10); + T e2m = T(1) - sq(es), + // To lowest order in e^2, taup = (1 - e^2) * tau = _e2m * tau; so use + // tau = taup/_e2m as a starting guess. (This starting guess is the + // geocentric latitude which, to first order in the flattening, is equal + // to the conformal latitude.) Only 1 iteration is needed for |lat| < + // 3.35 deg, otherwise 2 iterations are needed. If, instead, tau = taup + // is used the mean number of iterations increases to 1.99 (2 iterations + // are needed except near tau = 0). + tau = taup/e2m, + stol = tol * max(T(1), abs(taup)); + // min iterations = 1, max iterations = 2; mean = 1.94 + for (int i = 0; i < numit || GEOGRAPHICLIB_PANIC; ++i) { + T taupa = taupf(tau, es), + dtau = (taup - taupa) * (1 + e2m * sq(tau)) / + ( e2m * hypot(T(1), tau) * hypot(T(1), taupa) ); + tau += dtau; + if (!(abs(dtau) >= stol)) + break; + } + return tau; + } + + template bool Math::isfinite(T x) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::isfinite; return isfinite(x); +#else +#if defined(_MSC_VER) + return abs(x) <= (std::numeric_limits::max)(); +#else + // There's a problem using MPFR C++ 3.6.3 and g++ -std=c++14 (reported on + // 2015-05-04) with the parens around std::numeric_limits::max. Of + // course, these parens are only needed to deal with Windows stupidly + // defining max as a macro. So don't insert the parens on non-Windows + // platforms. + return abs(x) <= std::numeric_limits::max(); +#endif +#endif + } + + template T Math::NaN() { +#if defined(_MSC_VER) + return std::numeric_limits::has_quiet_NaN ? + std::numeric_limits::quiet_NaN() : + (std::numeric_limits::max)(); +#else + return std::numeric_limits::has_quiet_NaN ? + std::numeric_limits::quiet_NaN() : + std::numeric_limits::max(); +#endif + } + + template bool Math::isnan(T x) { +#if GEOGRAPHICLIB_CXX11_MATH + using std::isnan; return isnan(x); +#else + return x != x; +#endif + } + + template T Math::infinity() { +#if defined(_MSC_VER) + return std::numeric_limits::has_infinity ? + std::numeric_limits::infinity() : + (std::numeric_limits::max)(); +#else + return std::numeric_limits::has_infinity ? + std::numeric_limits::infinity() : + std::numeric_limits::max(); +#endif + } + + /// \cond SKIP + // Instantiate + template Math::real GEOGRAPHICLIB_EXPORT + Math::hypot(Math::real, Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::expm1(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::log1p(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::asinh(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::atanh(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::cbrt(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::remainder(Math::real, Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::remquo(Math::real, Math::real, int*); + template Math::real GEOGRAPHICLIB_EXPORT + Math::round(Math::real); + template long GEOGRAPHICLIB_EXPORT + Math::lround(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::copysign(Math::real, Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::fma(Math::real, Math::real, Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::sum(Math::real, Math::real, Math::real&); + template Math::real GEOGRAPHICLIB_EXPORT + Math::AngRound(Math::real); + template void GEOGRAPHICLIB_EXPORT + Math::sincosd(Math::real, Math::real&, Math::real&); + template Math::real GEOGRAPHICLIB_EXPORT + Math::sind(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::cosd(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::tand(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::atan2d(Math::real, Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::atand(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::eatanhe(Math::real, Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::taupf(Math::real, Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::tauf(Math::real, Math::real); + template bool GEOGRAPHICLIB_EXPORT + Math::isfinite(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::NaN(); + template bool GEOGRAPHICLIB_EXPORT + Math::isnan(Math::real); + template Math::real GEOGRAPHICLIB_EXPORT + Math::infinity(); +#if GEOGRAPHICLIB_PRECISION != 2 + // Always have double versions available + template double GEOGRAPHICLIB_EXPORT + Math::hypot(double, double); + template double GEOGRAPHICLIB_EXPORT + Math::expm1(double); + template double GEOGRAPHICLIB_EXPORT + Math::log1p(double); + template double GEOGRAPHICLIB_EXPORT + Math::asinh(double); + template double GEOGRAPHICLIB_EXPORT + Math::atanh(double); + template double GEOGRAPHICLIB_EXPORT + Math::cbrt(double); + template double GEOGRAPHICLIB_EXPORT + Math::remainder(double, double); + template double GEOGRAPHICLIB_EXPORT + Math::remquo(double, double, int*); + template double GEOGRAPHICLIB_EXPORT + Math::round(double); + template long GEOGRAPHICLIB_EXPORT + Math::lround(double); + template double GEOGRAPHICLIB_EXPORT + Math::copysign(double, double); + template double GEOGRAPHICLIB_EXPORT + Math::fma(double, double, double); + template double GEOGRAPHICLIB_EXPORT + Math::sum(double, double, double&); + template double GEOGRAPHICLIB_EXPORT + Math::AngRound(double); + template void GEOGRAPHICLIB_EXPORT + Math::sincosd(double, double&, double&); + template double GEOGRAPHICLIB_EXPORT + Math::sind(double); + template double GEOGRAPHICLIB_EXPORT + Math::cosd(double); + template double GEOGRAPHICLIB_EXPORT + Math::tand(double); + template double GEOGRAPHICLIB_EXPORT + Math::atan2d(double, double); + template double GEOGRAPHICLIB_EXPORT + Math::atand(double); + template double GEOGRAPHICLIB_EXPORT + Math::eatanhe(double, double); + template double GEOGRAPHICLIB_EXPORT + Math::taupf(double, double); + template double GEOGRAPHICLIB_EXPORT + Math::tauf(double, double); + template bool GEOGRAPHICLIB_EXPORT + Math::isfinite(double); + template double GEOGRAPHICLIB_EXPORT + Math::NaN(); + template bool GEOGRAPHICLIB_EXPORT + Math::isnan(double); + template double GEOGRAPHICLIB_EXPORT + Math::infinity(); +#endif +#if GEOGRAPHICLIB_HAVE_LONG_DOUBLE && GEOGRAPHICLIB_PRECISION != 3 + // And always have long double versions available (as long as long double is + // a really different from double). + template long double GEOGRAPHICLIB_EXPORT + Math::hypot(long double, long double); + template long double GEOGRAPHICLIB_EXPORT + Math::expm1(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::log1p(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::asinh(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::atanh(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::cbrt(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::remainder(long double, long double); + template long double GEOGRAPHICLIB_EXPORT + Math::remquo(long double, long double, int*); + template long double GEOGRAPHICLIB_EXPORT + Math::round(long double); + template long GEOGRAPHICLIB_EXPORT + Math::lround(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::copysign(long double, long double); + template long double GEOGRAPHICLIB_EXPORT + Math::fma(long double, long double, long double); + template long double GEOGRAPHICLIB_EXPORT + Math::sum(long double, long double, long double&); + template long double GEOGRAPHICLIB_EXPORT + Math::AngRound(long double); + template void GEOGRAPHICLIB_EXPORT + Math::sincosd(long double, long double&, long double&); + template long double GEOGRAPHICLIB_EXPORT + Math::sind(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::cosd(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::tand(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::atan2d(long double, long double); + template long double GEOGRAPHICLIB_EXPORT + Math::atand(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::eatanhe(long double, long double); + template long double GEOGRAPHICLIB_EXPORT + Math::taupf(long double, long double); + template long double GEOGRAPHICLIB_EXPORT + Math::tauf(long double, long double); + template bool GEOGRAPHICLIB_EXPORT + Math::isfinite(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::NaN(); + template bool GEOGRAPHICLIB_EXPORT + Math::isnan(long double); + template long double GEOGRAPHICLIB_EXPORT + Math::infinity(); +#endif + // Also we need int versions for Utility::nummatch + template int GEOGRAPHICLIB_EXPORT Math::NaN(); + template int GEOGRAPHICLIB_EXPORT Math::infinity(); + /// \endcond + +} // namespace GeographicLib diff --git a/src/Geo/Math.hpp b/src/Geo/Math.hpp new file mode 100644 index 000000000..4dbd1457e --- /dev/null +++ b/src/Geo/Math.hpp @@ -0,0 +1,647 @@ +/** + * \file Math.hpp + * \brief Header for GeographicLib::Math class + * + * Copyright (c) Charles Karney (2008-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +// Constants.hpp includes Math.hpp. Place this include outside Math.hpp's +// include guard to enforce this ordering. +#include "Constants.hpp" + +#if !defined(GEOGRAPHICLIB_MATH_HPP) +#define GEOGRAPHICLIB_MATH_HPP 1 + +/** + * Are C++11 math functions available? + **********************************************************************/ +#if !defined(GEOGRAPHICLIB_CXX11_MATH) +// Recent versions of g++ -std=c++11 (4.7 and later?) set __cplusplus to 201103 +// and support the new C++11 mathematical functions, std::atanh, etc. However +// the Android toolchain, which uses g++ -std=c++11 (4.8 as of 2014-03-11, +// according to Pullan Lu), does not support std::atanh. Android toolchains +// might define __ANDROID__ or ANDROID; so need to check both. With OSX the +// version is GNUC version 4.2 and __cplusplus is set to 201103, so remove the +// version check on GNUC. +# if defined(__GNUC__) && __cplusplus >= 201103 && \ + !(defined(__ANDROID__) || defined(ANDROID) || defined(__CYGWIN__)) +# define GEOGRAPHICLIB_CXX11_MATH 1 +// Visual C++ 12 supports these functions +# elif defined(_MSC_VER) && _MSC_VER >= 1800 +# define GEOGRAPHICLIB_CXX11_MATH 1 +# else +# define GEOGRAPHICLIB_CXX11_MATH 0 +# endif +#endif + +#if !defined(GEOGRAPHICLIB_WORDS_BIGENDIAN) +# define GEOGRAPHICLIB_WORDS_BIGENDIAN 0 +#endif + +#if !defined(GEOGRAPHICLIB_HAVE_LONG_DOUBLE) +# define GEOGRAPHICLIB_HAVE_LONG_DOUBLE 0 +#endif + +#if !defined(GEOGRAPHICLIB_PRECISION) +/** + * The precision of floating point numbers used in %GeographicLib. 1 means + * float (single precision); 2 (the default) means double; 3 means long double; + * 4 is reserved for quadruple precision. Nearly all the testing has been + * carried out with doubles and that's the recommended configuration. In order + * for long double to be used, GEOGRAPHICLIB_HAVE_LONG_DOUBLE needs to be + * defined. Note that with Microsoft Visual Studio, long double is the same as + * double. + **********************************************************************/ +# define GEOGRAPHICLIB_PRECISION 2 +#endif + +#include +#include +#include + +#if GEOGRAPHICLIB_PRECISION == 4 +#include +#include +#include +#elif GEOGRAPHICLIB_PRECISION == 5 +#include +#endif + +#if GEOGRAPHICLIB_PRECISION > 3 +// volatile keyword makes no sense for multiprec types +#define GEOGRAPHICLIB_VOLATILE +// Signal a convergence failure with multiprec types by throwing an exception +// at loop exit. +#define GEOGRAPHICLIB_PANIC \ + (throw GeographicLib::GeographicErr("Convergence failure"), false) +#else +#define GEOGRAPHICLIB_VOLATILE volatile +// Ignore convergence failures with standard floating points types by allowing +// loop to exit cleanly. +#define GEOGRAPHICLIB_PANIC false +#endif + +namespace GeographicLib { + + /** + * \brief Mathematical functions needed by %GeographicLib + * + * Define mathematical functions in order to localize system dependencies and + * to provide generic versions of the functions. In addition define a real + * type to be used by %GeographicLib. + * + * Example of use: + * \include example-Math.cpp + **********************************************************************/ + class Math { + private: + void dummy(); // Static check for GEOGRAPHICLIB_PRECISION + Math(); // Disable constructor + public: + +#if GEOGRAPHICLIB_HAVE_LONG_DOUBLE + /** + * The extended precision type for real numbers, used for some testing. + * This is long double on computers with this type; otherwise it is double. + **********************************************************************/ + typedef long double extended; +#else + typedef double extended; +#endif + +#if GEOGRAPHICLIB_PRECISION == 2 + /** + * The real type for %GeographicLib. Nearly all the testing has been done + * with \e real = double. However, the algorithms should also work with + * float and long double (where available). (CAUTION: reasonable + * accuracy typically cannot be obtained using floats.) + **********************************************************************/ + typedef double real; +#elif GEOGRAPHICLIB_PRECISION == 1 + typedef float real; +#elif GEOGRAPHICLIB_PRECISION == 3 + typedef extended real; +#elif GEOGRAPHICLIB_PRECISION == 4 + typedef boost::multiprecision::float128 real; +#elif GEOGRAPHICLIB_PRECISION == 5 + typedef mpfr::mpreal real; +#else + typedef double real; +#endif + + /** + * @return the number of bits of precision in a real number. + **********************************************************************/ + static int digits(); + + /** + * Set the binary precision of a real number. + * + * @param[in] ndigits the number of bits of precision. + * @return the resulting number of bits of precision. + * + * This only has an effect when GEOGRAPHICLIB_PRECISION = 5. See also + * Utility::set_digits for caveats about when this routine should be + * called. + **********************************************************************/ + static int set_digits(int ndigits); + + /** + * @return the number of decimal digits of precision in a real number. + **********************************************************************/ + static int digits10(); + + /** + * Number of additional decimal digits of precision for real relative to + * double (0 for float). + **********************************************************************/ + static int extra_digits(); + + /** + * true if the machine is big-endian. + **********************************************************************/ + static const bool bigendian = GEOGRAPHICLIB_WORDS_BIGENDIAN; + + /** + * @tparam T the type of the returned value. + * @return π. + **********************************************************************/ + template static T pi() { + using std::atan2; + static const T pi = atan2(T(0), T(-1)); + return pi; + } + /** + * A synonym for pi(). + **********************************************************************/ + static real pi() { return pi(); } + + /** + * @tparam T the type of the returned value. + * @return the number of radians in a degree. + **********************************************************************/ + template static T degree() { + static const T degree = pi() / 180; + return degree; + } + /** + * A synonym for degree(). + **********************************************************************/ + static real degree() { return degree(); } + + /** + * Square a number. + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @return x2. + **********************************************************************/ + template static T sq(T x) + { return x * x; } + + /** + * The hypotenuse function avoiding underflow and overflow. + * + * @tparam T the type of the arguments and the returned value. + * @param[in] x + * @param[in] y + * @return sqrt(x2 + y2). + **********************************************************************/ + template static T hypot(T x, T y); + + /** + * exp(\e x) − 1 accurate near \e x = 0. + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @return exp(\e x) − 1. + **********************************************************************/ + template static T expm1(T x); + + /** + * log(1 + \e x) accurate near \e x = 0. + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @return log(1 + \e x). + **********************************************************************/ + template static T log1p(T x); + + /** + * The inverse hyperbolic sine function. + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @return asinh(\e x). + **********************************************************************/ + template static T asinh(T x); + + /** + * The inverse hyperbolic tangent function. + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @return atanh(\e x). + **********************************************************************/ + template static T atanh(T x); + + /** + * Copy the sign. + * + * @tparam T the type of the argument. + * @param[in] x gives the magitude of the result. + * @param[in] y gives the sign of the result. + * @return value with the magnitude of \e x and with the sign of \e y. + * + * This routine correctly handles the case \e y = −0, returning + * &minus|x|. + **********************************************************************/ + template static T copysign(T x, T y); + + /** + * The cube root function. + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @return the real cube root of \e x. + **********************************************************************/ + template static T cbrt(T x); + + /** + * The remainder function. + * + * @tparam T the type of the arguments and the returned value. + * @param[in] x + * @param[in] y + * @return the remainder of \e x/\e y in the range [−\e y/2, \e y/2]. + **********************************************************************/ + template static T remainder(T x, T y); + + /** + * The remquo function. + * + * @tparam T the type of the arguments and the returned value. + * @param[in] x + * @param[in] y + * @param[out] n the low 3 bits of the quotient + * @return the remainder of \e x/\e y in the range [−\e y/2, \e y/2]. + **********************************************************************/ + template static T remquo(T x, T y, int* n); + + /** + * The round function. + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @return \e x round to the nearest integer (ties round away from 0). + **********************************************************************/ + template static T round(T x); + + /** + * The lround function. + * + * @tparam T the type of the argument. + * @param[in] x + * @return \e x round to the nearest integer as a long int (ties round away + * from 0). + * + * If the result does not fit in a long int, the return value is undefined. + **********************************************************************/ + template static long lround(T x); + + /** + * Fused multiply and add. + * + * @tparam T the type of the arguments and the returned value. + * @param[in] x + * @param[in] y + * @param[in] z + * @return xy + z, correctly rounded (on those platforms with + * support for the fma instruction). + * + * On platforms without the fma instruction, no attempt is + * made to improve on the result of a rounded multiplication followed by a + * rounded addition. + **********************************************************************/ + template static T fma(T x, T y, T z); + + /** + * Normalize a two-vector. + * + * @tparam T the type of the argument and the returned value. + * @param[in,out] x on output set to x/hypot(x, y). + * @param[in,out] y on output set to y/hypot(x, y). + **********************************************************************/ + template static void norm(T& x, T& y) + { T h = hypot(x, y); x /= h; y /= h; } + + /** + * The error-free sum of two numbers. + * + * @tparam T the type of the argument and the returned value. + * @param[in] u + * @param[in] v + * @param[out] t the exact error given by (\e u + \e v) - \e s. + * @return \e s = round(\e u + \e v). + * + * See D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B. (Note that \e t can be + * the same as one of the first two arguments.) + **********************************************************************/ + template static T sum(T u, T v, T& t); + + /** + * Evaluate a polynomial. + * + * @tparam T the type of the arguments and returned value. + * @param[in] N the order of the polynomial. + * @param[in] p the coefficient array (of size \e N + 1). + * @param[in] x the variable. + * @return the value of the polynomial. + * + * Evaluate y = ∑n=0..N + * pn xN−n. + * Return 0 if \e N < 0. Return p0, if \e N = 0 (even + * if \e x is infinite or a nan). The evaluation uses Horner's method. + **********************************************************************/ + template static T polyval(int N, const T p[], T x) + // This used to employ Math::fma; but that's too slow and it seemed not to + // improve the accuracy noticeably. This might change when there's direct + // hardware support for fma. + { T y = N < 0 ? 0 : *p++; while (--N >= 0) y = y * x + *p++; return y; } + + /** + * Normalize an angle. + * + * @tparam T the type of the argument and returned value. + * @param[in] x the angle in degrees. + * @return the angle reduced to the range (−180°, 180°]. + * + * The range of \e x is unrestricted. + **********************************************************************/ + template static T AngNormalize(T x) { + x = remainder(x, T(360)); return x != -180 ? x : 180; + } + + /** + * Normalize a latitude. + * + * @tparam T the type of the argument and returned value. + * @param[in] x the angle in degrees. + * @return x if it is in the range [−90°, 90°], otherwise + * return NaN. + **********************************************************************/ + template static T LatFix(T x) + { using std::abs; return abs(x) > 90 ? NaN() : x; } + + /** + * The exact difference of two angles reduced to + * (−180°, 180°]. + * + * @tparam T the type of the arguments and returned value. + * @param[in] x the first angle in degrees. + * @param[in] y the second angle in degrees. + * @param[out] e the error term in degrees. + * @return \e d, the truncated value of \e y − \e x. + * + * This computes \e z = \e y − \e x exactly, reduced to + * (−180°, 180°]; and then sets \e z = \e d + \e e where \e d + * is the nearest representable number to \e z and \e e is the truncation + * error. If \e d = −180, then \e e > 0; If \e d = 180, then \e e + * ≤ 0. + **********************************************************************/ + template static T AngDiff(T x, T y, T& e) { + T t, d = AngNormalize(sum(remainder(-x, T(360)), + remainder( y, T(360)), t)); + // Here y - x = d + t (mod 360), exactly, where d is in (-180,180] and + // abs(t) <= eps (eps = 2^-45 for doubles). The only case where the + // addition of t takes the result outside the range (-180,180] is d = 180 + // and t > 0. The case, d = -180 + eps, t = -eps, can't happen, since + // sum would have returned the exact result in such a case (i.e., given t + // = 0). + return sum(d == 180 && t > 0 ? -180 : d, t, e); + } + + /** + * Difference of two angles reduced to [−180°, 180°] + * + * @tparam T the type of the arguments and returned value. + * @param[in] x the first angle in degrees. + * @param[in] y the second angle in degrees. + * @return \e y − \e x, reduced to the range [−180°, + * 180°]. + * + * The result is equivalent to computing the difference exactly, reducing + * it to (−180°, 180°] and rounding the result. Note that + * this prescription allows −180° to be returned (e.g., if \e x + * is tiny and negative and \e y = 180°). + **********************************************************************/ + template static T AngDiff(T x, T y) + { T e; return AngDiff(x, y, e); } + + /** + * Coarsen a value close to zero. + * + * @tparam T the type of the argument and returned value. + * @param[in] x + * @return the coarsened value. + * + * The makes the smallest gap in \e x = 1/16 − nextafter(1/16, 0) = + * 1/257 for reals = 0.7 pm on the earth if \e x is an angle in + * degrees. (This is about 1000 times more resolution than we get with + * angles around 90°.) We use this to avoid having to deal with near + * singular cases when \e x is non-zero but tiny (e.g., + * 10−200). This converts −0 to +0; however tiny + * negative numbers get converted to −0. + **********************************************************************/ + template static T AngRound(T x); + + /** + * Evaluate the sine and cosine function with the argument in degrees + * + * @tparam T the type of the arguments. + * @param[in] x in degrees. + * @param[out] sinx sin(x). + * @param[out] cosx cos(x). + * + * The results obey exactly the elementary properties of the trigonometric + * functions, e.g., sin 9° = cos 81° = − sin 123456789°. + * If x = −0, then \e sinx = −0; this is the only case where + * −0 is returned. + **********************************************************************/ + template static void sincosd(T x, T& sinx, T& cosx); + + /** + * Evaluate the sine function with the argument in degrees + * + * @tparam T the type of the argument and the returned value. + * @param[in] x in degrees. + * @return sin(x). + **********************************************************************/ + template static T sind(T x); + + /** + * Evaluate the cosine function with the argument in degrees + * + * @tparam T the type of the argument and the returned value. + * @param[in] x in degrees. + * @return cos(x). + **********************************************************************/ + template static T cosd(T x); + + /** + * Evaluate the tangent function with the argument in degrees + * + * @tparam T the type of the argument and the returned value. + * @param[in] x in degrees. + * @return tan(x). + * + * If \e x = ±90°, then a suitably large (but finite) value is + * returned. + **********************************************************************/ + template static T tand(T x); + + /** + * Evaluate the atan2 function with the result in degrees + * + * @tparam T the type of the arguments and the returned value. + * @param[in] y + * @param[in] x + * @return atan2(y, x) in degrees. + * + * The result is in the range (−180° 180°]. N.B., + * atan2d(±0, −1) = +180°; atan2d(−ε, + * −1) = −180°, for ε positive and tiny; + * atan2d(±0, +1) = ±0°. + **********************************************************************/ + template static T atan2d(T y, T x); + + /** + * Evaluate the atan function with the result in degrees + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @return atan(x) in degrees. + **********************************************************************/ + template static T atand(T x); + + /** + * Evaluate e atanh(e x) + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @param[in] es the signed eccentricity = sign(e2) + * sqrt(|e2|) + * @return e atanh(e x) + * + * If e2 is negative (e is imaginary), the + * expression is evaluated in terms of atan. + **********************************************************************/ + template static T eatanhe(T x, T es); + + /** + * tanχ in terms of tanφ + * + * @tparam T the type of the argument and the returned value. + * @param[in] tau τ = tanφ + * @param[in] es the signed eccentricity = sign(e2) + * sqrt(|e2|) + * @return τ′ = tanχ + * + * See Eqs. (7--9) of + * C. F. F. Karney, + * + * Transverse Mercator with an accuracy of a few nanometers, + * J. Geodesy 85(8), 475--485 (Aug. 2011) + * (preprint + * arXiv:1002.1417). + **********************************************************************/ + template static T taupf(T tau, T es); + + /** + * tanφ in terms of tanχ + * + * @tparam T the type of the argument and the returned value. + * @param[in] taup τ′ = tanχ + * @param[in] es the signed eccentricity = sign(e2) + * sqrt(|e2|) + * @return τ = tanφ + * + * See Eqs. (19--21) of + * C. F. F. Karney, + * + * Transverse Mercator with an accuracy of a few nanometers, + * J. Geodesy 85(8), 475--485 (Aug. 2011) + * (preprint + * arXiv:1002.1417). + **********************************************************************/ + template static T tauf(T taup, T es); + + /** + * Test for finiteness. + * + * @tparam T the type of the argument. + * @param[in] x + * @return true if number is finite, false if NaN or infinite. + **********************************************************************/ + template static bool isfinite(T x); + + /** + * The NaN (not a number) + * + * @tparam T the type of the returned value. + * @return NaN if available, otherwise return the max real of type T. + **********************************************************************/ + template static T NaN(); + + /** + * A synonym for NaN(). + **********************************************************************/ + static real NaN() { return NaN(); } + + /** + * Test for NaN. + * + * @tparam T the type of the argument. + * @param[in] x + * @return true if argument is a NaN. + **********************************************************************/ + template static bool isnan(T x); + + /** + * Infinity + * + * @tparam T the type of the returned value. + * @return infinity if available, otherwise return the max real. + **********************************************************************/ + template static T infinity(); + + /** + * A synonym for infinity(). + **********************************************************************/ + static real infinity() { return infinity(); } + + /** + * Swap the bytes of a quantity + * + * @tparam T the type of the argument and the returned value. + * @param[in] x + * @return x with its bytes swapped. + **********************************************************************/ + template static T swab(T x) { + union { + T r; + unsigned char c[sizeof(T)]; + } b; + b.r = x; + for (int i = sizeof(T)/2; i--; ) + std::swap(b.c[i], b.c[sizeof(T) - 1 - i]); + return b.r; + } + + }; + +} // namespace GeographicLib + +#endif // GEOGRAPHICLIB_MATH_HPP diff --git a/src/Geo/PolarStereographic.cpp b/src/Geo/PolarStereographic.cpp new file mode 100644 index 000000000..f110730d3 --- /dev/null +++ b/src/Geo/PolarStereographic.cpp @@ -0,0 +1,109 @@ +/** + * \file PolarStereographic.cpp + * \brief Implementation for GeographicLib::PolarStereographic class + * + * Copyright (c) Charles Karney (2008-2017) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#include "PolarStereographic.hpp" + +namespace GeographicLib { + + using namespace std; + + PolarStereographic::PolarStereographic(real a, real f, real k0) + : _a(a) + , _f(f) + , _e2(_f * (2 - _f)) + , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2))) + , _e2m(1 - _e2) + , _c( (1 - _f) * exp(Math::eatanhe(real(1), _es)) ) + , _k0(k0) + { + if (!(Math::isfinite(_a) && _a > 0)) + throw GeographicErr("Equatorial radius is not positive"); + if (!(Math::isfinite(_f) && _f < 1)) + throw GeographicErr("Polar semi-axis is not positive"); + if (!(Math::isfinite(_k0) && _k0 > 0)) + throw GeographicErr("Scale is not positive"); + } + + const PolarStereographic& PolarStereographic::UPS() { + static const PolarStereographic ups(Constants::WGS84_a(), + Constants::WGS84_f(), + Constants::UPS_k0()); + return ups; + } + + // This formulation converts to conformal coordinates by tau = tan(phi) and + // tau' = tan(phi') where phi' is the conformal latitude. The formulas are: + // tau = tan(phi) + // secphi = hypot(1, tau) + // sig = sinh(e * atanh(e * tau / secphi)) + // taup = tan(phip) = tau * hypot(1, sig) - sig * hypot(1, tau) + // c = (1 - f) * exp(e * atanh(e)) + // + // Forward: + // rho = (2*k0*a/c) / (hypot(1, taup) + taup) (taup >= 0) + // = (2*k0*a/c) * (hypot(1, taup) - taup) (taup < 0) + // + // Reverse: + // taup = ((2*k0*a/c) / rho - rho / (2*k0*a/c))/2 + // + // Scale: + // k = (rho/a) * secphi * sqrt((1-e2) + e2 / secphi^2) + // + // In limit rho -> 0, tau -> inf, taup -> inf, secphi -> inf, secphip -> inf + // secphip = taup = exp(-e * atanh(e)) * tau = exp(-e * atanh(e)) * secphi + + void PolarStereographic::Forward(bool northp, real lat, real lon, + real& x, real& y, + real& gamma, real& k) const { + lat = Math::LatFix(lat); + lat *= northp ? 1 : -1; + real + tau = Math::tand(lat), + secphi = Math::hypot(real(1), tau), + taup = Math::taupf(tau, _es), + rho = Math::hypot(real(1), taup) + abs(taup); + rho = taup >= 0 ? (lat != 90 ? 1/rho : 0) : rho; + rho *= 2 * _k0 * _a / _c; + k = lat != 90 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) : + _k0; + Math::sincosd(lon, x, y); + x *= rho; + y *= (northp ? -rho : rho); + gamma = Math::AngNormalize(northp ? lon : -lon); + } + + void PolarStereographic::Reverse(bool northp, real x, real y, + real& lat, real& lon, + real& gamma, real& k) const { + real + rho = Math::hypot(x, y), + t = rho != 0 ? rho / (2 * _k0 * _a / _c) : + Math::sq(numeric_limits::epsilon()), + taup = (1 / t - t) / 2, + tau = Math::tauf(taup, _es), + secphi = Math::hypot(real(1), tau); + k = rho != 0 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) : + _k0; + lat = (northp ? 1 : -1) * Math::atand(tau); + lon = Math::atan2d(x, northp ? -y : y ); + gamma = Math::AngNormalize(northp ? lon : -lon); + } + + void PolarStereographic::SetScale(real lat, real k) { + if (!(Math::isfinite(k) && k > 0)) + throw GeographicErr("Scale is not positive"); + if (!(-90 < lat && lat <= 90)) + throw GeographicErr("Latitude must be in (-90d, 90d]"); + real x, y, gamma, kold; + _k0 = 1; + Forward(true, lat, 0, x, y, gamma, kold); + _k0 *= k/kold; + } + +} // namespace GeographicLib diff --git a/src/Geo/PolarStereographic.hpp b/src/Geo/PolarStereographic.hpp new file mode 100644 index 000000000..823edd644 --- /dev/null +++ b/src/Geo/PolarStereographic.hpp @@ -0,0 +1,160 @@ +/** + * \file PolarStereographic.hpp + * \brief Header for GeographicLib::PolarStereographic class + * + * Copyright (c) Charles Karney (2008-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#if !defined(GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP) +#define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP 1 + +#include "Constants.hpp" + +namespace GeographicLib { + + /** + * \brief Polar stereographic projection + * + * Implementation taken from the report, + * - J. P. Snyder, + * Map Projections: A + * Working Manual, USGS Professional Paper 1395 (1987), + * pp. 160--163. + * + * This is a straightforward implementation of the equations in Snyder except + * that Newton's method is used to invert the projection. + * + * This class also returns the meridian convergence \e gamma and scale \e k. + * The meridian convergence is the bearing of grid north (the \e y axis) + * measured clockwise from true north. + * + * Example of use: + * \include example-PolarStereographic.cpp + **********************************************************************/ + class GEOGRAPHICLIB_EXPORT PolarStereographic { + private: + typedef Math::real real; + real _a, _f, _e2, _es, _e2m, _c; + real _k0; + public: + + /** + * Constructor for a ellipsoid with + * + * @param[in] a equatorial radius (meters). + * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. + * Negative \e f gives a prolate ellipsoid. + * @param[in] k0 central scale factor. + * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k0 is + * not positive. + **********************************************************************/ + PolarStereographic(real a, real f, real k0); + + /** + * Set the scale for the projection. + * + * @param[in] lat (degrees) assuming \e northp = true. + * @param[in] k scale at latitude \e lat (default 1). + * @exception GeographicErr \e k is not positive. + * @exception GeographicErr if \e lat is not in (−90°, + * 90°]. + **********************************************************************/ + void SetScale(real lat, real k = real(1)); + + /** + * Forward projection, from geographic to polar stereographic. + * + * @param[in] northp the pole which is the center of projection (true means + * north, false means south). + * @param[in] lat latitude of point (degrees). + * @param[in] lon longitude of point (degrees). + * @param[out] x easting of point (meters). + * @param[out] y northing of point (meters). + * @param[out] gamma meridian convergence at point (degrees). + * @param[out] k scale of projection at point. + * + * No false easting or northing is added. \e lat should be in the range + * (−90°, 90°] for \e northp = true and in the range + * [−90°, 90°) for \e northp = false. + **********************************************************************/ + void Forward(bool northp, real lat, real lon, + real& x, real& y, real& gamma, real& k) const; + + /** + * Reverse projection, from polar stereographic to geographic. + * + * @param[in] northp the pole which is the center of projection (true means + * north, false means south). + * @param[in] x easting of point (meters). + * @param[in] y northing of point (meters). + * @param[out] lat latitude of point (degrees). + * @param[out] lon longitude of point (degrees). + * @param[out] gamma meridian convergence at point (degrees). + * @param[out] k scale of projection at point. + * + * No false easting or northing is added. The value of \e lon returned is + * in the range [−180°, 180°]. + **********************************************************************/ + void Reverse(bool northp, real x, real y, + real& lat, real& lon, real& gamma, real& k) const; + + /** + * PolarStereographic::Forward without returning the convergence and scale. + **********************************************************************/ + void Forward(bool northp, real lat, real lon, + real& x, real& y) const { + real gamma, k; + Forward(northp, lat, lon, x, y, gamma, k); + } + + /** + * PolarStereographic::Reverse without returning the convergence and scale. + **********************************************************************/ + void Reverse(bool northp, real x, real y, + real& lat, real& lon) const { + real gamma, k; + Reverse(northp, x, y, lat, lon, gamma, k); + } + + /** \name Inspector functions + **********************************************************************/ + ///@{ + /** + * @return \e a the equatorial radius of the ellipsoid (meters). This is + * the value used in the constructor. + **********************************************************************/ + Math::real EquatorialRadius() const { return _a; } + + /** + * @return \e f the flattening of the ellipsoid. This is the value used in + * the constructor. + **********************************************************************/ + Math::real Flattening() const { return _f; } + + /** + * The central scale for the projection. This is the value of \e k0 used + * in the constructor and is the scale at the pole unless overridden by + * PolarStereographic::SetScale. + **********************************************************************/ + Math::real CentralScale() const { return _k0; } + + /** + * \deprecated An old name for EquatorialRadius(). + **********************************************************************/ + // GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()") + Math::real MajorRadius() const { return EquatorialRadius(); } + ///@} + + /** + * A global instantiation of PolarStereographic with the WGS84 ellipsoid + * and the UPS scale factor. However, unlike UPS, no false easting or + * northing is added. + **********************************************************************/ + static const PolarStereographic& UPS(); + }; + +} // namespace GeographicLib + +#endif // GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP diff --git a/src/QGCGeo.cc b/src/Geo/QGCGeo.cc similarity index 100% rename from src/QGCGeo.cc rename to src/Geo/QGCGeo.cc diff --git a/src/QGCGeo.h b/src/Geo/QGCGeo.h similarity index 100% rename from src/QGCGeo.h rename to src/Geo/QGCGeo.h diff --git a/src/Geo/TransverseMercator.cpp b/src/Geo/TransverseMercator.cpp new file mode 100644 index 000000000..5f2ecf383 --- /dev/null +++ b/src/Geo/TransverseMercator.cpp @@ -0,0 +1,596 @@ +/** + * \file TransverseMercator.cpp + * \brief Implementation for GeographicLib::TransverseMercator class + * + * Copyright (c) Charles Karney (2008-2017) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + * + * This implementation follows closely JHS 154, ETRS89 - + * järjestelmään liittyvät karttaprojektiot, + * tasokoordinaatistot ja karttalehtijako (Map projections, plane + * coordinates, and map sheet index for ETRS89), published by JUHTA, Finnish + * Geodetic Institute, and the National Land Survey of Finland (2006). + * + * The relevant section is available as the 2008 PDF file + * http://docs.jhs-suositukset.fi/jhs-suositukset/JHS154/JHS154_liite1.pdf + * + * This is a straight transcription of the formulas in this paper with the + * following exceptions: + * - use of 6th order series instead of 4th order series. This reduces the + * error to about 5nm for the UTM range of coordinates (instead of 200nm), + * with a speed penalty of only 1%; + * - use Newton's method instead of plain iteration to solve for latitude in + * terms of isometric latitude in the Reverse method; + * - use of Horner's representation for evaluating polynomials and Clenshaw's + * method for summing trigonometric series; + * - several modifications of the formulas to improve the numerical accuracy; + * - evaluating the convergence and scale using the expression for the + * projection or its inverse. + * + * If the preprocessor variable GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER is set + * to an integer between 4 and 8, then this specifies the order of the series + * used for the forward and reverse transformations. The default value is 6. + * (The series accurate to 12th order is given in \ref tmseries.) + **********************************************************************/ + +#include +#include +#include "TransverseMercator.hpp" + +namespace GeographicLib { + + using namespace std; + + TransverseMercator::TransverseMercator(real a, real f, real k0) + : _a(a) + , _f(f) + , _k0(k0) + , _e2(_f * (2 - _f)) + , _es((f < 0 ? -1 : 1) * sqrt(abs(_e2))) + , _e2m(1 - _e2) + // _c = sqrt( pow(1 + _e, 1 + _e) * pow(1 - _e, 1 - _e) ) ) + // See, for example, Lee (1976), p 100. + , _c( sqrt(_e2m) * exp(Math::eatanhe(real(1), _es)) ) + , _n(_f / (2 - _f)) + { + if (!(Math::isfinite(_a) && _a > 0)) + throw GeographicErr("Equatorial radius is not positive"); + if (!(Math::isfinite(_f) && _f < 1)) + throw GeographicErr("Polar semi-axis is not positive"); + if (!(Math::isfinite(_k0) && _k0 > 0)) + throw GeographicErr("Scale is not positive"); + + // Generated by Maxima on 2015-05-14 22:55:13-04:00 +#if GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2 + static const real b1coeff[] = { + // b1*(n+1), polynomial in n2 of order 2 + 1, 16, 64, 64, + }; // count = 4 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3 + static const real b1coeff[] = { + // b1*(n+1), polynomial in n2 of order 3 + 1, 4, 64, 256, 256, + }; // count = 5 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4 + static const real b1coeff[] = { + // b1*(n+1), polynomial in n2 of order 4 + 25, 64, 256, 4096, 16384, 16384, + }; // count = 6 +#else +#error "Bad value for GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER" +#endif + +#if GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 + static const real alpcoeff[] = { + // alp[1]/n^1, polynomial in n of order 3 + 164, 225, -480, 360, 720, + // alp[2]/n^2, polynomial in n of order 2 + 557, -864, 390, 1440, + // alp[3]/n^3, polynomial in n of order 1 + -1236, 427, 1680, + // alp[4]/n^4, polynomial in n of order 0 + 49561, 161280, + }; // count = 14 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 + static const real alpcoeff[] = { + // alp[1]/n^1, polynomial in n of order 4 + -635, 328, 450, -960, 720, 1440, + // alp[2]/n^2, polynomial in n of order 3 + 4496, 3899, -6048, 2730, 10080, + // alp[3]/n^3, polynomial in n of order 2 + 15061, -19776, 6832, 26880, + // alp[4]/n^4, polynomial in n of order 1 + -171840, 49561, 161280, + // alp[5]/n^5, polynomial in n of order 0 + 34729, 80640, + }; // count = 20 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 + static const real alpcoeff[] = { + // alp[1]/n^1, polynomial in n of order 5 + 31564, -66675, 34440, 47250, -100800, 75600, 151200, + // alp[2]/n^2, polynomial in n of order 4 + -1983433, 863232, 748608, -1161216, 524160, 1935360, + // alp[3]/n^3, polynomial in n of order 3 + 670412, 406647, -533952, 184464, 725760, + // alp[4]/n^4, polynomial in n of order 2 + 6601661, -7732800, 2230245, 7257600, + // alp[5]/n^5, polynomial in n of order 1 + -13675556, 3438171, 7983360, + // alp[6]/n^6, polynomial in n of order 0 + 212378941, 319334400, + }; // count = 27 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 + static const real alpcoeff[] = { + // alp[1]/n^1, polynomial in n of order 6 + 1804025, 2020096, -4267200, 2204160, 3024000, -6451200, 4838400, 9676800, + // alp[2]/n^2, polynomial in n of order 5 + 4626384, -9917165, 4316160, 3743040, -5806080, 2620800, 9676800, + // alp[3]/n^3, polynomial in n of order 4 + -67102379, 26816480, 16265880, -21358080, 7378560, 29030400, + // alp[4]/n^4, polynomial in n of order 3 + 155912000, 72618271, -85060800, 24532695, 79833600, + // alp[5]/n^5, polynomial in n of order 2 + 102508609, -109404448, 27505368, 63866880, + // alp[6]/n^6, polynomial in n of order 1 + -12282192400LL, 2760926233LL, 4151347200LL, + // alp[7]/n^7, polynomial in n of order 0 + 1522256789, 1383782400, + }; // count = 35 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 + static const real alpcoeff[] = { + // alp[1]/n^1, polynomial in n of order 7 + -75900428, 37884525, 42422016, -89611200, 46287360, 63504000, -135475200, + 101606400, 203212800, + // alp[2]/n^2, polynomial in n of order 6 + 148003883, 83274912, -178508970, 77690880, 67374720, -104509440, + 47174400, 174182400, + // alp[3]/n^3, polynomial in n of order 5 + 318729724, -738126169, 294981280, 178924680, -234938880, 81164160, + 319334400, + // alp[4]/n^4, polynomial in n of order 4 + -40176129013LL, 14967552000LL, 6971354016LL, -8165836800LL, 2355138720LL, + 7664025600LL, + // alp[5]/n^5, polynomial in n of order 3 + 10421654396LL, 3997835751LL, -4266773472LL, 1072709352, 2490808320LL, + // alp[6]/n^6, polynomial in n of order 2 + 175214326799LL, -171950693600LL, 38652967262LL, 58118860800LL, + // alp[7]/n^7, polynomial in n of order 1 + -67039739596LL, 13700311101LL, 12454041600LL, + // alp[8]/n^8, polynomial in n of order 0 + 1424729850961LL, 743921418240LL, + }; // count = 44 +#else +#error "Bad value for GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER" +#endif + +#if GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4 + static const real betcoeff[] = { + // bet[1]/n^1, polynomial in n of order 3 + -4, 555, -960, 720, 1440, + // bet[2]/n^2, polynomial in n of order 2 + -437, 96, 30, 1440, + // bet[3]/n^3, polynomial in n of order 1 + -148, 119, 3360, + // bet[4]/n^4, polynomial in n of order 0 + 4397, 161280, + }; // count = 14 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5 + static const real betcoeff[] = { + // bet[1]/n^1, polynomial in n of order 4 + -3645, -64, 8880, -15360, 11520, 23040, + // bet[2]/n^2, polynomial in n of order 3 + 4416, -3059, 672, 210, 10080, + // bet[3]/n^3, polynomial in n of order 2 + -627, -592, 476, 13440, + // bet[4]/n^4, polynomial in n of order 1 + -3520, 4397, 161280, + // bet[5]/n^5, polynomial in n of order 0 + 4583, 161280, + }; // count = 20 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6 + static const real betcoeff[] = { + // bet[1]/n^1, polynomial in n of order 5 + 384796, -382725, -6720, 932400, -1612800, 1209600, 2419200, + // bet[2]/n^2, polynomial in n of order 4 + -1118711, 1695744, -1174656, 258048, 80640, 3870720, + // bet[3]/n^3, polynomial in n of order 3 + 22276, -16929, -15984, 12852, 362880, + // bet[4]/n^4, polynomial in n of order 2 + -830251, -158400, 197865, 7257600, + // bet[5]/n^5, polynomial in n of order 1 + -435388, 453717, 15966720, + // bet[6]/n^6, polynomial in n of order 0 + 20648693, 638668800, + }; // count = 27 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7 + static const real betcoeff[] = { + // bet[1]/n^1, polynomial in n of order 6 + -5406467, 6156736, -6123600, -107520, 14918400, -25804800, 19353600, + 38707200, + // bet[2]/n^2, polynomial in n of order 5 + 829456, -5593555, 8478720, -5873280, 1290240, 403200, 19353600, + // bet[3]/n^3, polynomial in n of order 4 + 9261899, 3564160, -2708640, -2557440, 2056320, 58060800, + // bet[4]/n^4, polynomial in n of order 3 + 14928352, -9132761, -1742400, 2176515, 79833600, + // bet[5]/n^5, polynomial in n of order 2 + -8005831, -1741552, 1814868, 63866880, + // bet[6]/n^6, polynomial in n of order 1 + -261810608, 268433009, 8302694400LL, + // bet[7]/n^7, polynomial in n of order 0 + 219941297, 5535129600LL, + }; // count = 35 +#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8 + static const real betcoeff[] = { + // bet[1]/n^1, polynomial in n of order 7 + 31777436, -37845269, 43097152, -42865200, -752640, 104428800, -180633600, + 135475200, 270950400, + // bet[2]/n^2, polynomial in n of order 6 + 24749483, 14930208, -100683990, 152616960, -105719040, 23224320, 7257600, + 348364800, + // bet[3]/n^3, polynomial in n of order 5 + -232468668, 101880889, 39205760, -29795040, -28131840, 22619520, + 638668800, + // bet[4]/n^4, polynomial in n of order 4 + 324154477, 1433121792, -876745056, -167270400, 208945440, 7664025600LL, + // bet[5]/n^5, polynomial in n of order 3 + 457888660, -312227409, -67920528, 70779852, 2490808320LL, + // bet[6]/n^6, polynomial in n of order 2 + -19841813847LL, -3665348512LL, 3758062126LL, 116237721600LL, + // bet[7]/n^7, polynomial in n of order 1 + -1989295244, 1979471673, 49816166400LL, + // bet[8]/n^8, polynomial in n of order 0 + 191773887257LL, 3719607091200LL, + }; // count = 44 +#else +#error "Bad value for GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER" +#endif + + GEOGRAPHICLIB_STATIC_ASSERT(sizeof(b1coeff) / sizeof(real) == + maxpow_/2 + 2, + "Coefficient array size mismatch for b1"); + GEOGRAPHICLIB_STATIC_ASSERT(sizeof(alpcoeff) / sizeof(real) == + (maxpow_ * (maxpow_ + 3))/2, + "Coefficient array size mismatch for alp"); + GEOGRAPHICLIB_STATIC_ASSERT(sizeof(betcoeff) / sizeof(real) == + (maxpow_ * (maxpow_ + 3))/2, + "Coefficient array size mismatch for bet"); + int m = maxpow_/2; + _b1 = Math::polyval(m, b1coeff, Math::sq(_n)) / (b1coeff[m + 1] * (1+_n)); + // _a1 is the equivalent radius for computing the circumference of + // ellipse. + _a1 = _b1 * _a; + int o = 0; + real d = _n; + for (int l = 1; l <= maxpow_; ++l) { + m = maxpow_ - l; + _alp[l] = d * Math::polyval(m, alpcoeff + o, _n) / alpcoeff[o + m + 1]; + _bet[l] = d * Math::polyval(m, betcoeff + o, _n) / betcoeff[o + m + 1]; + o += m + 2; + d *= _n; + } + // Post condition: o == sizeof(alpcoeff) / sizeof(real) && + // o == sizeof(betcoeff) / sizeof(real) + } + + const TransverseMercator& TransverseMercator::UTM() { + static const TransverseMercator utm(Constants::WGS84_a(), + Constants::WGS84_f(), + Constants::UTM_k0()); + return utm; + } + + // Engsager and Poder (2007) use trigonometric series to convert between phi + // and phip. Here are the series... + // + // Conversion from phi to phip: + // + // phip = phi + sum(c[j] * sin(2*j*phi), j, 1, 6) + // + // c[1] = - 2 * n + // + 2/3 * n^2 + // + 4/3 * n^3 + // - 82/45 * n^4 + // + 32/45 * n^5 + // + 4642/4725 * n^6; + // c[2] = 5/3 * n^2 + // - 16/15 * n^3 + // - 13/9 * n^4 + // + 904/315 * n^5 + // - 1522/945 * n^6; + // c[3] = - 26/15 * n^3 + // + 34/21 * n^4 + // + 8/5 * n^5 + // - 12686/2835 * n^6; + // c[4] = 1237/630 * n^4 + // - 12/5 * n^5 + // - 24832/14175 * n^6; + // c[5] = - 734/315 * n^5 + // + 109598/31185 * n^6; + // c[6] = 444337/155925 * n^6; + // + // Conversion from phip to phi: + // + // phi = phip + sum(d[j] * sin(2*j*phip), j, 1, 6) + // + // d[1] = 2 * n + // - 2/3 * n^2 + // - 2 * n^3 + // + 116/45 * n^4 + // + 26/45 * n^5 + // - 2854/675 * n^6; + // d[2] = 7/3 * n^2 + // - 8/5 * n^3 + // - 227/45 * n^4 + // + 2704/315 * n^5 + // + 2323/945 * n^6; + // d[3] = 56/15 * n^3 + // - 136/35 * n^4 + // - 1262/105 * n^5 + // + 73814/2835 * n^6; + // d[4] = 4279/630 * n^4 + // - 332/35 * n^5 + // - 399572/14175 * n^6; + // d[5] = 4174/315 * n^5 + // - 144838/6237 * n^6; + // d[6] = 601676/22275 * n^6; + // + // In order to maintain sufficient relative accuracy close to the pole use + // + // S = sum(c[i]*sin(2*i*phi),i,1,6) + // taup = (tau + tan(S)) / (1 - tau * tan(S)) + + // In Math::taupf and Math::tauf we evaluate the forward transform explicitly + // and solve the reverse one by Newton's method. + // + // There are adapted from TransverseMercatorExact (taup and taupinv). tau = + // tan(phi), taup = sinh(psi) + + void TransverseMercator::Forward(real lon0, real lat, real lon, + real& x, real& y, + real& gamma, real& k) const { + lat = Math::LatFix(lat); + lon = Math::AngDiff(lon0, lon); + // Explicitly enforce the parity + int + latsign = (lat < 0) ? -1 : 1, + lonsign = (lon < 0) ? -1 : 1; + lon *= lonsign; + lat *= latsign; + bool backside = lon > 90; + if (backside) { + if (lat == 0) + latsign = -1; + lon = 180 - lon; + } + real sphi, cphi, slam, clam; + Math::sincosd(lat, sphi, cphi); + Math::sincosd(lon, slam, clam); + // phi = latitude + // phi' = conformal latitude + // psi = isometric latitude + // tau = tan(phi) + // tau' = tan(phi') + // [xi', eta'] = Gauss-Schreiber TM coordinates + // [xi, eta] = Gauss-Krueger TM coordinates + // + // We use + // tan(phi') = sinh(psi) + // sin(phi') = tanh(psi) + // cos(phi') = sech(psi) + // denom^2 = 1-cos(phi')^2*sin(lam)^2 = 1-sech(psi)^2*sin(lam)^2 + // sin(xip) = sin(phi')/denom = tanh(psi)/denom + // cos(xip) = cos(phi')*cos(lam)/denom = sech(psi)*cos(lam)/denom + // cosh(etap) = 1/denom = 1/denom + // sinh(etap) = cos(phi')*sin(lam)/denom = sech(psi)*sin(lam)/denom + real etap, xip; + if (lat != 90) { + real + tau = sphi / cphi, + taup = Math::taupf(tau, _es); + xip = atan2(taup, clam); + // Used to be + // etap = Math::atanh(sin(lam) / cosh(psi)); + etap = Math::asinh(slam / Math::hypot(taup, clam)); + // convergence and scale for Gauss-Schreiber TM (xip, etap) -- gamma0 = + // atan(tan(xip) * tanh(etap)) = atan(tan(lam) * sin(phi')); + // sin(phi') = tau'/sqrt(1 + tau'^2) + // Krueger p 22 (44) + gamma = Math::atan2d(slam * taup, clam * Math::hypot(real(1), taup)); + // k0 = sqrt(1 - _e2 * sin(phi)^2) * (cos(phi') / cos(phi)) * cosh(etap) + // Note 1/cos(phi) = cosh(psip); + // and cos(phi') * cosh(etap) = 1/hypot(sinh(psi), cos(lam)) + // + // This form has cancelling errors. This property is lost if cosh(psip) + // is replaced by 1/cos(phi), even though it's using "primary" data (phi + // instead of psip). + k = sqrt(_e2m + _e2 * Math::sq(cphi)) * Math::hypot(real(1), tau) + / Math::hypot(taup, clam); + } else { + xip = Math::pi()/2; + etap = 0; + gamma = lon; + k = _c; + } + // {xi',eta'} is {northing,easting} for Gauss-Schreiber transverse Mercator + // (for eta' = 0, xi' = bet). {xi,eta} is {northing,easting} for transverse + // Mercator with constant scale on the central meridian (for eta = 0, xip = + // rectifying latitude). Define + // + // zeta = xi + i*eta + // zeta' = xi' + i*eta' + // + // The conversion from conformal to rectifying latitude can be expressed as + // a series in _n: + // + // zeta = zeta' + sum(h[j-1]' * sin(2 * j * zeta'), j = 1..maxpow_) + // + // where h[j]' = O(_n^j). The reversion of this series gives + // + // zeta' = zeta - sum(h[j-1] * sin(2 * j * zeta), j = 1..maxpow_) + // + // which is used in Reverse. + // + // Evaluate sums via Clenshaw method. See + // https://en.wikipedia.org/wiki/Clenshaw_algorithm + // + // Let + // + // S = sum(a[k] * phi[k](x), k = 0..n) + // phi[k+1](x) = alpha[k](x) * phi[k](x) + beta[k](x) * phi[k-1](x) + // + // Evaluate S with + // + // b[n+2] = b[n+1] = 0 + // b[k] = alpha[k](x) * b[k+1] + beta[k+1](x) * b[k+2] + a[k] + // S = (a[0] + beta[1](x) * b[2]) * phi[0](x) + b[1] * phi[1](x) + // + // Here we have + // + // x = 2 * zeta' + // phi[k](x) = sin(k * x) + // alpha[k](x) = 2 * cos(x) + // beta[k](x) = -1 + // [ sin(A+B) - 2*cos(B)*sin(A) + sin(A-B) = 0, A = k*x, B = x ] + // n = maxpow_ + // a[k] = _alp[k] + // S = b[1] * sin(x) + // + // For the derivative we have + // + // x = 2 * zeta' + // phi[k](x) = cos(k * x) + // alpha[k](x) = 2 * cos(x) + // beta[k](x) = -1 + // [ cos(A+B) - 2*cos(B)*cos(A) + cos(A-B) = 0, A = k*x, B = x ] + // a[0] = 1; a[k] = 2*k*_alp[k] + // S = (a[0] - b[2]) + b[1] * cos(x) + // + // Matrix formulation (not used here): + // phi[k](x) = [sin(k * x); k * cos(k * x)] + // alpha[k](x) = 2 * [cos(x), 0; -sin(x), cos(x)] + // beta[k](x) = -1 * [1, 0; 0, 1] + // a[k] = _alp[k] * [1, 0; 0, 1] + // b[n+2] = b[n+1] = [0, 0; 0, 0] + // b[k] = alpha[k](x) * b[k+1] + beta[k+1](x) * b[k+2] + a[k] + // N.B., for all k: b[k](1,2) = 0; b[k](1,1) = b[k](2,2) + // S = (a[0] + beta[1](x) * b[2]) * phi[0](x) + b[1] * phi[1](x) + // phi[0](x) = [0; 0] + // phi[1](x) = [sin(x); cos(x)] + real + c0 = cos(2 * xip), ch0 = cosh(2 * etap), + s0 = sin(2 * xip), sh0 = sinh(2 * etap); + complex a(2 * c0 * ch0, -2 * s0 * sh0); // 2 * cos(2*zeta') + int n = maxpow_; + complex + y0(n & 1 ? _alp[n] : 0), y1, // default initializer is 0+i0 + z0(n & 1 ? 2*n * _alp[n] : 0), z1; + if (n & 1) --n; + while (n) { + y1 = a * y0 - y1 + _alp[n]; + z1 = a * z0 - z1 + 2*n * _alp[n]; + --n; + y0 = a * y1 - y0 + _alp[n]; + z0 = a * z1 - z0 + 2*n * _alp[n]; + --n; + } + a /= real(2); // cos(2*zeta') + z1 = real(1) - z1 + a * z0; + a = complex(s0 * ch0, c0 * sh0); // sin(2*zeta') + y1 = complex(xip, etap) + a * y0; + // Fold in change in convergence and scale for Gauss-Schreiber TM to + // Gauss-Krueger TM. + gamma -= Math::atan2d(z1.imag(), z1.real()); + k *= _b1 * abs(z1); + real xi = y1.real(), eta = y1.imag(); + y = _a1 * _k0 * (backside ? Math::pi() - xi : xi) * latsign; + x = _a1 * _k0 * eta * lonsign; + if (backside) + gamma = 180 - gamma; + gamma *= latsign * lonsign; + gamma = Math::AngNormalize(gamma); + k *= _k0; + } + + void TransverseMercator::Reverse(real lon0, real x, real y, + real& lat, real& lon, + real& gamma, real& k) const { + // This undoes the steps in Forward. The wrinkles are: (1) Use of the + // reverted series to express zeta' in terms of zeta. (2) Newton's method + // to solve for phi in terms of tan(phi). + real + xi = y / (_a1 * _k0), + eta = x / (_a1 * _k0); + // Explicitly enforce the parity + int + xisign = (xi < 0) ? -1 : 1, + etasign = (eta < 0) ? -1 : 1; + xi *= xisign; + eta *= etasign; + bool backside = xi > Math::pi()/2; + if (backside) + xi = Math::pi() - xi; + real + c0 = cos(2 * xi), ch0 = cosh(2 * eta), + s0 = sin(2 * xi), sh0 = sinh(2 * eta); + complex a(2 * c0 * ch0, -2 * s0 * sh0); // 2 * cos(2*zeta) + int n = maxpow_; + complex + y0(n & 1 ? -_bet[n] : 0), y1, // default initializer is 0+i0 + z0(n & 1 ? -2*n * _bet[n] : 0), z1; + if (n & 1) --n; + while (n) { + y1 = a * y0 - y1 - _bet[n]; + z1 = a * z0 - z1 - 2*n * _bet[n]; + --n; + y0 = a * y1 - y0 - _bet[n]; + z0 = a * z1 - z0 - 2*n * _bet[n]; + --n; + } + a /= real(2); // cos(2*zeta) + z1 = real(1) - z1 + a * z0; + a = complex(s0 * ch0, c0 * sh0); // sin(2*zeta) + y1 = complex(xi, eta) + a * y0; + // Convergence and scale for Gauss-Schreiber TM to Gauss-Krueger TM. + gamma = Math::atan2d(z1.imag(), z1.real()); + k = _b1 / abs(z1); + // JHS 154 has + // + // phi' = asin(sin(xi') / cosh(eta')) (Krueger p 17 (25)) + // lam = asin(tanh(eta') / cos(phi') + // psi = asinh(tan(phi')) + real + xip = y1.real(), etap = y1.imag(), + s = sinh(etap), + c = max(real(0), cos(xip)), // cos(pi/2) might be negative + r = Math::hypot(s, c); + if (r != 0) { + lon = Math::atan2d(s, c); // Krueger p 17 (25) + // Use Newton's method to solve for tau + real + sxip = sin(xip), + tau = Math::tauf(sxip/r, _es); + gamma += Math::atan2d(sxip * tanh(etap), c); // Krueger p 19 (31) + lat = Math::atand(tau); + // Note cos(phi') * cosh(eta') = r + k *= sqrt(_e2m + _e2 / (1 + Math::sq(tau))) * + Math::hypot(real(1), tau) * r; + } else { + lat = 90; + lon = 0; + k *= _c; + } + lat *= xisign; + if (backside) + lon = 180 - lon; + lon *= etasign; + lon = Math::AngNormalize(lon + lon0); + if (backside) + gamma = 180 - gamma; + gamma *= xisign * etasign; + gamma = Math::AngNormalize(gamma); + k *= _k0; + } + +} // namespace GeographicLib diff --git a/src/Geo/TransverseMercator.hpp b/src/Geo/TransverseMercator.hpp new file mode 100644 index 000000000..42ad91149 --- /dev/null +++ b/src/Geo/TransverseMercator.hpp @@ -0,0 +1,206 @@ +/** + * \file TransverseMercator.hpp + * \brief Header for GeographicLib::TransverseMercator class + * + * Copyright (c) Charles Karney (2008-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP) +#define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP 1 + +#include "Constants.hpp" + +#if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER) +/** + * The order of the series approximation used in TransverseMercator. + * GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER can be set to any integer in [4, 8]. + **********************************************************************/ +# define GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER \ + (GEOGRAPHICLIB_PRECISION == 2 ? 6 : \ + (GEOGRAPHICLIB_PRECISION == 1 ? 4 : 8)) +#endif + +namespace GeographicLib { + + /** + * \brief Transverse Mercator projection + * + * This uses Krüger's method which evaluates the projection and its + * inverse in terms of a series. See + * - L. Krüger, + * Konforme + * Abbildung des Erdellipsoids in der Ebene (Conformal mapping of the + * ellipsoidal earth to the plane), Royal Prussian Geodetic Institute, New + * Series 52, 172 pp. (1912). + * - C. F. F. Karney, + * + * Transverse Mercator with an accuracy of a few nanometers, + * J. Geodesy 85(8), 475--485 (Aug. 2011); + * preprint + * arXiv:1002.1417. + * + * Krüger's method has been extended from 4th to 6th order. The maximum + * error is 5 nm (5 nanometers), ground distance, for all positions within 35 + * degrees of the central meridian. The error in the convergence is 2 + * × 10−15" and the relative error in the scale + * is 6 × 10−12%%. See Sec. 4 of + * arXiv:1002.1417 for details. + * The speed penalty in going to 6th order is only about 1%. + * + * There's a singularity in the projection at φ = 0°, λ + * − λ0 = ±(1 − \e e)90° (≈ + * ±82.6° for the WGS84 ellipsoid), where \e e is the + * eccentricity. Beyond this point, the series ceases to converge and the + * results from this method will be garbage. To be on the safe side, don't + * use this method if the angular distance from the central meridian exceeds + * (1 − 2e)90° (≈ 75° for the WGS84 ellipsoid) + * + * TransverseMercatorExact is an alternative implementation of the projection + * using exact formulas which yield accurate (to 8 nm) results over the + * entire ellipsoid. + * + * The ellipsoid parameters and the central scale are set in the constructor. + * The central meridian (which is a trivial shift of the longitude) is + * specified as the \e lon0 argument of the TransverseMercator::Forward and + * TransverseMercator::Reverse functions. The latitude of origin is taken to + * be the equator. There is no provision in this class for specifying a + * false easting or false northing or a different latitude of origin. + * However these are can be simply included by the calling function. For + * example, the UTMUPS class applies the false easting and false northing for + * the UTM projections. A more complicated example is the British National + * Grid ( + * EPSG:7405) which requires the use of a latitude of origin. This is + * implemented by the GeographicLib::OSGB class. + * + * This class also returns the meridian convergence \e gamma and scale \e k. + * The meridian convergence is the bearing of grid north (the \e y axis) + * measured clockwise from true north. + * + * See TransverseMercator.cpp for more information on the implementation. + * + * See \ref transversemercator for a discussion of this projection. + * + * Example of use: + * \include example-TransverseMercator.cpp + * + * TransverseMercatorProj is a + * command-line utility providing access to the functionality of + * TransverseMercator and TransverseMercatorExact. + **********************************************************************/ + + class GEOGRAPHICLIB_EXPORT TransverseMercator { + private: + typedef Math::real real; + static const int maxpow_ = GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER; + static const int numit_ = 5; + real _a, _f, _k0, _e2, _es, _e2m, _c, _n; + // _alp[0] and _bet[0] unused + real _a1, _b1, _alp[maxpow_ + 1], _bet[maxpow_ + 1]; + friend class Ellipsoid; // For access to taupf, tauf. + public: + + /** + * Constructor for a ellipsoid with + * + * @param[in] a equatorial radius (meters). + * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. + * Negative \e f gives a prolate ellipsoid. + * @param[in] k0 central scale factor. + * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k0 is + * not positive. + **********************************************************************/ + TransverseMercator(real a, real f, real k0); + + /** + * Forward projection, from geographic to transverse Mercator. + * + * @param[in] lon0 central meridian of the projection (degrees). + * @param[in] lat latitude of point (degrees). + * @param[in] lon longitude of point (degrees). + * @param[out] x easting of point (meters). + * @param[out] y northing of point (meters). + * @param[out] gamma meridian convergence at point (degrees). + * @param[out] k scale of projection at point. + * + * No false easting or northing is added. \e lat should be in the range + * [−90°, 90°]. + **********************************************************************/ + void Forward(real lon0, real lat, real lon, + real& x, real& y, real& gamma, real& k) const; + + /** + * Reverse projection, from transverse Mercator to geographic. + * + * @param[in] lon0 central meridian of the projection (degrees). + * @param[in] x easting of point (meters). + * @param[in] y northing of point (meters). + * @param[out] lat latitude of point (degrees). + * @param[out] lon longitude of point (degrees). + * @param[out] gamma meridian convergence at point (degrees). + * @param[out] k scale of projection at point. + * + * No false easting or northing is added. The value of \e lon returned is + * in the range [−180°, 180°]. + **********************************************************************/ + void Reverse(real lon0, real x, real y, + real& lat, real& lon, real& gamma, real& k) const; + + /** + * TransverseMercator::Forward without returning the convergence and scale. + **********************************************************************/ + void Forward(real lon0, real lat, real lon, + real& x, real& y) const { + real gamma, k; + Forward(lon0, lat, lon, x, y, gamma, k); + } + + /** + * TransverseMercator::Reverse without returning the convergence and scale. + **********************************************************************/ + void Reverse(real lon0, real x, real y, + real& lat, real& lon) const { + real gamma, k; + Reverse(lon0, x, y, lat, lon, gamma, k); + } + + /** \name Inspector functions + **********************************************************************/ + ///@{ + /** + * @return \e a the equatorial radius of the ellipsoid (meters). This is + * the value used in the constructor. + **********************************************************************/ + Math::real EquatorialRadius() const { return _a; } + + /** + * @return \e f the flattening of the ellipsoid. This is the value used in + * the constructor. + **********************************************************************/ + Math::real Flattening() const { return _f; } + + /** + * @return \e k0 central scale for the projection. This is the value of \e + * k0 used in the constructor and is the scale on the central meridian. + **********************************************************************/ + Math::real CentralScale() const { return _k0; } + + /** + * \deprecated An old name for EquatorialRadius(). + **********************************************************************/ + // GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()") + Math::real MajorRadius() const { return EquatorialRadius(); } + ///@} + + /** + * A global instantiation of TransverseMercator with the WGS84 ellipsoid + * and the UTM scale factor. However, unlike UTM, no false easting or + * northing is added. + **********************************************************************/ + static const TransverseMercator& UTM(); + }; + +} // namespace GeographicLib + +#endif // GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP diff --git a/src/UTM.cpp b/src/Geo/UTM.cpp similarity index 100% rename from src/UTM.cpp rename to src/Geo/UTM.cpp diff --git a/src/UTM.h b/src/Geo/UTM.h similarity index 100% rename from src/UTM.h rename to src/Geo/UTM.h diff --git a/src/Geo/UTMUPS.cpp b/src/Geo/UTMUPS.cpp new file mode 100644 index 000000000..5fde1a636 --- /dev/null +++ b/src/Geo/UTMUPS.cpp @@ -0,0 +1,297 @@ +/** + * \file UTMUPS.cpp + * \brief Implementation for GeographicLib::UTMUPS class + * + * Copyright (c) Charles Karney (2008-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#include "UTMUPS.hpp" +#include "MGRS.hpp" +#include "PolarStereographic.hpp" +#include "TransverseMercator.hpp" +#include "Utility.hpp" + +namespace GeographicLib { + + using namespace std; + + const int UTMUPS::falseeasting_[] = + { MGRS::upseasting_ * MGRS::tile_, MGRS::upseasting_ * MGRS::tile_, + MGRS::utmeasting_ * MGRS::tile_, MGRS::utmeasting_ * MGRS::tile_ }; + const int UTMUPS::falsenorthing_[] = + { MGRS::upseasting_ * MGRS::tile_, MGRS::upseasting_ * MGRS::tile_, + MGRS::maxutmSrow_ * MGRS::tile_, MGRS::minutmNrow_ * MGRS::tile_ }; + const int UTMUPS::mineasting_[] = + { MGRS::minupsSind_ * MGRS::tile_, MGRS::minupsNind_ * MGRS::tile_, + MGRS::minutmcol_ * MGRS::tile_, MGRS::minutmcol_ * MGRS::tile_ }; + const int UTMUPS::maxeasting_[] = + { MGRS::maxupsSind_ * MGRS::tile_, MGRS::maxupsNind_ * MGRS::tile_, + MGRS::maxutmcol_ * MGRS::tile_, MGRS::maxutmcol_ * MGRS::tile_ }; + const int UTMUPS::minnorthing_[] = + { MGRS::minupsSind_ * MGRS::tile_, MGRS::minupsNind_ * MGRS::tile_, + MGRS::minutmSrow_ * MGRS::tile_, + (MGRS::minutmNrow_ + MGRS::minutmSrow_ - MGRS::maxutmSrow_) + * MGRS::tile_ }; + const int UTMUPS::maxnorthing_[] = + { MGRS::maxupsSind_ * MGRS::tile_, MGRS::maxupsNind_ * MGRS::tile_, + (MGRS::maxutmSrow_ + MGRS::maxutmNrow_ - MGRS::minutmNrow_) * + MGRS::tile_, + MGRS::maxutmNrow_ * MGRS::tile_ }; + + int UTMUPS::StandardZone(real lat, real lon, int setzone) { + if (!(setzone >= MINPSEUDOZONE && setzone <= MAXZONE)) + throw GeographicErr("Illegal zone requested " + Utility::str(setzone)); + if (setzone >= MINZONE || setzone == INVALID) + return setzone; + if (Math::isnan(lat) || Math::isnan(lon)) // Check if lat or lon is a NaN + return INVALID; + if (setzone == UTM || (lat >= -80 && lat < 84)) { + int ilon = int(floor(Math::AngNormalize(lon))); + if (ilon == 180) ilon = -180; // ilon now in [-180,180) + int zone = (ilon + 186)/6; + int band = MGRS::LatitudeBand(lat); + if (band == 7 && zone == 31 && ilon >= 3) // The Norway exception + zone = 32; + else if (band == 9 && ilon >= 0 && ilon < 42) // The Svalbard exception + zone = 2 * ((ilon + 183)/12) + 1; + return zone; + } else + return UPS; + } + + void UTMUPS::Forward(real lat, real lon, + int& zone, bool& northp, real& x, real& y, + real& gamma, real& k, + int setzone, bool mgrslimits) { + if (abs(lat) > 90) + throw GeographicErr("Latitude " + Utility::str(lat) + + "d not in [-90d, 90d]"); + bool northp1 = lat >= 0; + int zone1 = StandardZone(lat, lon, setzone); + if (zone1 == INVALID) { + zone = zone1; + northp = northp1; + x = y = gamma = k = Math::NaN(); + return; + } + real x1, y1, gamma1, k1; + bool utmp = zone1 != UPS; + if (utmp) { + real + lon0 = CentralMeridian(zone1), + dlon = lon - lon0; + dlon = abs(dlon - 360 * floor((dlon + 180)/360)); + if (!(dlon <= 60)) + // Check isn't really necessary because CheckCoords catches this case. + // But this allows a more meaningful error message to be given. + throw GeographicErr("Longitude " + Utility::str(lon) + + "d more than 60d from center of UTM zone " + + Utility::str(zone1)); + TransverseMercator::UTM().Forward(lon0, lat, lon, x1, y1, gamma1, k1); + } else { + if (abs(lat) < 70) + // Check isn't really necessary ... (see above). + throw GeographicErr("Latitude " + Utility::str(lat) + + "d more than 20d from " + + (northp1 ? "N" : "S") + " pole"); + PolarStereographic::UPS().Forward(northp1, lat, lon, x1, y1, gamma1, k1); + } + int ind = (utmp ? 2 : 0) + (northp1 ? 1 : 0); + x1 += falseeasting_[ind]; + y1 += falsenorthing_[ind]; + if (! CheckCoords(zone1 != UPS, northp1, x1, y1, mgrslimits, false) ) + throw GeographicErr("Latitude " + Utility::str(lat) + + ", longitude " + Utility::str(lon) + + " out of legal range for " + + (utmp ? "UTM zone " + Utility::str(zone1) : + "UPS")); + zone = zone1; + northp = northp1; + x = x1; + y = y1; + gamma = gamma1; + k = k1; + } + + void UTMUPS::Reverse(int zone, bool northp, real x, real y, + real& lat, real& lon, real& gamma, real& k, + bool mgrslimits) { + if (zone == INVALID || Math::isnan(x) || Math::isnan(y)) { + lat = lon = gamma = k = Math::NaN(); + return; + } + if (!(zone >= MINZONE && zone <= MAXZONE)) + throw GeographicErr("Zone " + Utility::str(zone) + + " not in range [0, 60]"); + bool utmp = zone != UPS; + CheckCoords(utmp, northp, x, y, mgrslimits); + int ind = (utmp ? 2 : 0) + (northp ? 1 : 0); + x -= falseeasting_[ind]; + y -= falsenorthing_[ind]; + if (utmp) + TransverseMercator::UTM().Reverse(CentralMeridian(zone), + x, y, lat, lon, gamma, k); + else + PolarStereographic::UPS().Reverse(northp, x, y, lat, lon, gamma, k); + } + + bool UTMUPS::CheckCoords(bool utmp, bool northp, real x, real y, + bool mgrslimits, bool throwp) { + // Limits are all multiples of 100km and are all closed on the both ends. + // Failure tests are such that NaNs succeed. + real slop = mgrslimits ? 0 : MGRS::tile_; + int ind = (utmp ? 2 : 0) + (northp ? 1 : 0); + if (x < mineasting_[ind] - slop || x > maxeasting_[ind] + slop) { + if (!throwp) return false; + throw GeographicErr("Easting " + Utility::str(x/1000) + "km not in " + + (mgrslimits ? "MGRS/" : "") + + (utmp ? "UTM" : "UPS") + " range for " + + (northp ? "N" : "S" ) + " hemisphere [" + + Utility::str((mineasting_[ind] - slop)/1000) + + "km, " + + Utility::str((maxeasting_[ind] + slop)/1000) + + "km]"); + } + if (y < minnorthing_[ind] - slop || y > maxnorthing_[ind] + slop) { + if (!throwp) return false; + throw GeographicErr("Northing " + Utility::str(y/1000) + "km not in " + + (mgrslimits ? "MGRS/" : "") + + (utmp ? "UTM" : "UPS") + " range for " + + (northp ? "N" : "S" ) + " hemisphere [" + + Utility::str((minnorthing_[ind] - slop)/1000) + + "km, " + + Utility::str((maxnorthing_[ind] + slop)/1000) + + "km]"); + } + return true; + } + + void UTMUPS::Transfer(int zonein, bool northpin, real xin, real yin, + int zoneout, bool northpout, real& xout, real& yout, + int& zone) { + bool northp = northpin; + if (zonein != zoneout) { + // Determine lat, lon + real lat, lon; + GeographicLib::UTMUPS::Reverse(zonein, northpin, xin, yin, lat, lon); + // Try converting to zoneout + real x, y; + int zone1; + GeographicLib::UTMUPS::Forward(lat, lon, zone1, northp, x, y, + zoneout == UTMUPS::MATCH + ? zonein : zoneout); + if (zone1 == 0 && northp != northpout) + throw GeographicErr + ("Attempt to transfer UPS coordinates between hemispheres"); + zone = zone1; + xout = x; + yout = y; + } else { + if (zoneout == 0 && northp != northpout) + throw GeographicErr + ("Attempt to transfer UPS coordinates between hemispheres"); + zone = zoneout; + xout = xin; + yout = yin; + } + if (northp != northpout) + // Can't get here if UPS + yout += (northpout ? -1 : 1) * MGRS::utmNshift_; + return; + } + + void UTMUPS::DecodeZone(const std::string& zonestr, int& zone, bool& northp) + { + unsigned zlen = unsigned(zonestr.size()); + if (zlen == 0) + throw GeographicErr("Empty zone specification"); + // Longest zone spec is 32north, 42south, invalid = 7 + if (zlen > 7) + throw GeographicErr("More than 7 characters in zone specification " + + zonestr); + + const char* c = zonestr.c_str(); + char* q; + int zone1 = strtol(c, &q, 10); + // if (zone1 == 0) zone1 = UPS; (not necessary) + + if (zone1 == UPS) { + if (!(q == c)) + // Don't allow 0n as an alternative to n for UPS coordinates + throw GeographicErr("Illegal zone 0 in " + zonestr + + ", use just the hemisphere for UPS"); + } else if (!(zone1 >= MINUTMZONE && zone1 <= MAXUTMZONE)) + throw GeographicErr("Zone " + Utility::str(zone1) + + " not in range [1, 60]"); + else if (!isdigit(zonestr[0])) + throw GeographicErr("Must use unsigned number for zone " + + Utility::str(zone1)); + else if (q - c > 2) + throw GeographicErr("More than 2 digits use to specify zone " + + Utility::str(zone1)); + + string hemi(zonestr, q - c); + for (std::string::iterator p = hemi.begin(); p != hemi.end(); ++p) + *p = char(std::tolower(*p)); + if (q == c && (hemi == "inv" || hemi == "invalid")) { + zone = INVALID; + northp = false; + return; + } + bool northp1 = hemi == "north" || hemi == "n"; + if (!(northp1 || hemi == "south" || hemi == "s")) + throw GeographicErr(string("Illegal hemisphere ") + hemi + " in " + + zonestr + ", specify north or south"); + zone = zone1; + northp = northp1; + } + + std::string UTMUPS::EncodeZone(int zone, bool northp, bool abbrev) { + if (zone == INVALID) + return string(abbrev ? "inv" : "invalid"); + if (!(zone >= MINZONE && zone <= MAXZONE)) + throw GeographicErr("Zone " + Utility::str(zone) + + " not in range [0, 60]"); + ostringstream os; + if (zone != UPS) + os << setfill('0') << setw(2) << zone; + if (abbrev) + os << (northp ? 'n' : 's'); + else + os << (northp ? "north" : "south"); + return os.str(); + } + + void UTMUPS::DecodeEPSG(int epsg, int& zone, bool& northp) { + northp = false; + if (epsg >= epsg01N && epsg <= epsg60N) { + zone = (epsg - epsg01N) + MINUTMZONE; + northp = true; + } else if (epsg == epsgN) { + zone = UPS; + northp = true; + } else if (epsg >= epsg01S && epsg <= epsg60S) { + zone = (epsg - epsg01S) + MINUTMZONE; + } else if (epsg == epsgS) { + zone = UPS; + } else { + zone = INVALID; + } + } + + int UTMUPS::EncodeEPSG(int zone, bool northp) { + int epsg = -1; + if (zone == UPS) + epsg = epsgS; + else if (zone >= MINUTMZONE && zone <= MAXUTMZONE) + epsg = (zone - MINUTMZONE) + epsg01S; + if (epsg >= 0 && northp) + epsg += epsgN - epsgS; + return epsg; + } + + Math::real UTMUPS::UTMShift() { return real(MGRS::utmNshift_); } + +} // namespace GeographicLib diff --git a/src/Geo/UTMUPS.hpp b/src/Geo/UTMUPS.hpp new file mode 100644 index 000000000..303386127 --- /dev/null +++ b/src/Geo/UTMUPS.hpp @@ -0,0 +1,428 @@ +/** + * \file UTMUPS.hpp + * \brief Header for GeographicLib::UTMUPS class + * + * Copyright (c) Charles Karney (2008-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#if !defined(GEOGRAPHICLIB_UTMUPS_HPP) +#define GEOGRAPHICLIB_UTMUPS_HPP 1 + +#include "Constants.hpp" + +namespace GeographicLib { + + /** + * \brief Convert between geographic coordinates and UTM/UPS + * + * UTM and UPS are defined + * - J. W. Hager, J. F. Behensky, and B. W. Drew, + * + * The Universal Grids: Universal Transverse Mercator (UTM) and Universal + * Polar Stereographic (UPS), Defense Mapping Agency, Technical Manual + * TM8358.2 (1989). + * . + * Section 2-3 defines UTM and section 3-2.4 defines UPS. This document also + * includes approximate algorithms for the computation of the underlying + * transverse Mercator and polar stereographic projections. Here we + * substitute much more accurate algorithms given by + * GeographicLib:TransverseMercator and GeographicLib:PolarStereographic. + * These are the algorithms recommended by the NGA document + * - + * The Universal Grids and the Transverse Mercator and Polar Stereographic + * Map Projections, NGA.SIG.0012_2.0.0_UTMUPS (2014). + * + * In this implementation, the conversions are closed, i.e., output from + * Forward is legal input for Reverse and vice versa. The error is about 5nm + * in each direction. However, the conversion from legal UTM/UPS coordinates + * to geographic coordinates and back might throw an error if the initial + * point is within 5nm of the edge of the allowed range for the UTM/UPS + * coordinates. + * + * The simplest way to guarantee the closed property is to define allowed + * ranges for the eastings and northings for UTM and UPS coordinates. The + * UTM boundaries are the same for all zones. (The only place the + * exceptional nature of the zone boundaries is evident is when converting to + * UTM/UPS coordinates requesting the standard zone.) The MGRS lettering + * scheme imposes natural limits on UTM/UPS coordinates which may be + * converted into MGRS coordinates. For the conversion to/from geographic + * coordinates these ranges have been extended by 100km in order to provide a + * generous overlap between UTM and UPS and between UTM zones. + * + * The NGA software package + * geotrans + * also provides conversions to and from UTM and UPS. Version 2.4.2 (and + * earlier) suffers from some drawbacks: + * - Inconsistent rules are used to determine the whether a particular UTM or + * UPS coordinate is legal. A more systematic approach is taken here. + * - The underlying projections are not very accurately implemented. + * + * The GeographicLib::UTMUPS::EncodeZone encodes the UTM zone and hemisphere + * to allow UTM/UPS coordinated to be displayed as, for example, "38N 444500 + * 3688500". According to NGA.SIG.0012_2.0.0_UTMUPS the use of "N" to denote + * "north" in the context is not allowed (since a upper case letter in this + * context denotes the MGRS latitude band). Consequently, as of version + * 1.36, EncodeZone uses the lower case letters "n" and "s" to denote the + * hemisphere. In addition EncodeZone accepts an optional final argument \e + * abbrev, which, if false, results in the hemisphere being spelled out as in + * "38north". + * + * Example of use: + * \include example-UTMUPS.cpp + **********************************************************************/ + class GEOGRAPHICLIB_EXPORT UTMUPS { + private: + typedef Math::real real; + static const int falseeasting_[4]; + static const int falsenorthing_[4]; + static const int mineasting_[4]; + static const int maxeasting_[4]; + static const int minnorthing_[4]; + static const int maxnorthing_[4]; + static const int epsg01N = 32601; // EPSG code for UTM 01N + static const int epsg60N = 32660; // EPSG code for UTM 60N + static const int epsgN = 32661; // EPSG code for UPS N + static const int epsg01S = 32701; // EPSG code for UTM 01S + static const int epsg60S = 32760; // EPSG code for UTM 60S + static const int epsgS = 32761; // EPSG code for UPS S + static real CentralMeridian(int zone) + { return real(6 * zone - 183); } + // Throw an error if easting or northing are outside standard ranges. If + // throwp = false, return bool instead. + static bool CheckCoords(bool utmp, bool northp, real x, real y, + bool msgrlimits = false, bool throwp = true); + UTMUPS(); // Disable constructor + + public: + + /** + * In this class we bring together the UTM and UPS coordinates systems. + * The UTM divides the earth between latitudes −80° and 84° + * into 60 zones numbered 1 thru 60. Zone assign zone number 0 to the UPS + * regions, covering the two poles. Within UTMUPS, non-negative zone + * numbers refer to one of the "physical" zones, 0 for UPS and [1, 60] for + * UTM. Negative "pseudo-zone" numbers are used to select one of the + * physical zones. + **********************************************************************/ + enum zonespec { + /** + * The smallest pseudo-zone number. + **********************************************************************/ + MINPSEUDOZONE = -4, + /** + * A marker for an undefined or invalid zone. Equivalent to NaN. + **********************************************************************/ + INVALID = -4, + /** + * If a coordinate already include zone information (e.g., it is an MGRS + * coordinate), use that, otherwise apply the UTMUPS::STANDARD rules. + **********************************************************************/ + MATCH = -3, + /** + * Apply the standard rules for UTM zone assigment extending the UTM zone + * to each pole to give a zone number in [1, 60]. For example, use UTM + * zone 38 for longitude in [42°, 48°). The rules include the + * Norway and Svalbard exceptions. + **********************************************************************/ + UTM = -2, + /** + * Apply the standard rules for zone assignment to give a zone number in + * [0, 60]. If the latitude is not in [−80°, 84°), then + * use UTMUPS::UPS = 0, otherwise apply the rules for UTMUPS::UTM. The + * tests on latitudes and longitudes are all closed on the lower end open + * on the upper. Thus for UTM zone 38, latitude is in [−80°, + * 84°) and longitude is in [42°, 48°). + **********************************************************************/ + STANDARD = -1, + /** + * The largest pseudo-zone number. + **********************************************************************/ + MAXPSEUDOZONE = -1, + /** + * The smallest physical zone number. + **********************************************************************/ + MINZONE = 0, + /** + * The zone number used for UPS + **********************************************************************/ + UPS = 0, + /** + * The smallest UTM zone number. + **********************************************************************/ + MINUTMZONE = 1, + /** + * The largest UTM zone number. + **********************************************************************/ + MAXUTMZONE = 60, + /** + * The largest physical zone number. + **********************************************************************/ + MAXZONE = 60, + }; + + /** + * The standard zone. + * + * @param[in] lat latitude (degrees). + * @param[in] lon longitude (degrees). + * @param[in] setzone zone override (optional). If omitted, use the + * standard rules for picking the zone. If \e setzone is given then use + * that zone if it is non-negative, otherwise apply the rules given in + * UTMUPS::zonespec. + * @exception GeographicErr if \e setzone is outside the range + * [UTMUPS::MINPSEUDOZONE, UTMUPS::MAXZONE] = [−4, 60]. + * + * This is exact. + **********************************************************************/ + static int StandardZone(real lat, real lon, int setzone = STANDARD); + + /** + * Forward projection, from geographic to UTM/UPS. + * + * @param[in] lat latitude of point (degrees). + * @param[in] lon longitude of point (degrees). + * @param[out] zone the UTM zone (zero means UPS). + * @param[out] northp hemisphere (true means north, false means south). + * @param[out] x easting of point (meters). + * @param[out] y northing of point (meters). + * @param[out] gamma meridian convergence at point (degrees). + * @param[out] k scale of projection at point. + * @param[in] setzone zone override (optional). + * @param[in] mgrslimits if true enforce the stricter MGRS limits on the + * coordinates (default = false). + * @exception GeographicErr if \e lat is not in [−90°, + * 90°]. + * @exception GeographicErr if the resulting \e x or \e y is out of allowed + * range (see Reverse); in this case, these arguments are unchanged. + * + * If \e setzone is omitted, use the standard rules for picking the zone. + * If \e setzone is given then use that zone if it is non-negative, + * otherwise apply the rules given in UTMUPS::zonespec. The accuracy of + * the conversion is about 5nm. + * + * The northing \e y jumps by UTMUPS::UTMShift() when crossing the equator + * in the southerly direction. Sometimes it is useful to remove this + * discontinuity in \e y by extending the "northern" hemisphere using + * UTMUPS::Transfer: + * \code + double lat = -1, lon = 123; + int zone; + bool northp; + double x, y, gamma, k; + GeographicLib::UTMUPS::Forward(lat, lon, zone, northp, x, y, gamma, k); + GeographicLib::UTMUPS::Transfer(zone, northp, x, y, + zone, true, x, y, zone); + northp = true; + \endcode + **********************************************************************/ + static void Forward(real lat, real lon, + int& zone, bool& northp, real& x, real& y, + real& gamma, real& k, + int setzone = STANDARD, bool mgrslimits = false); + + /** + * Reverse projection, from UTM/UPS to geographic. + * + * @param[in] zone the UTM zone (zero means UPS). + * @param[in] northp hemisphere (true means north, false means south). + * @param[in] x easting of point (meters). + * @param[in] y northing of point (meters). + * @param[out] lat latitude of point (degrees). + * @param[out] lon longitude of point (degrees). + * @param[out] gamma meridian convergence at point (degrees). + * @param[out] k scale of projection at point. + * @param[in] mgrslimits if true enforce the stricter MGRS limits on the + * coordinates (default = false). + * @exception GeographicErr if \e zone, \e x, or \e y is out of allowed + * range; this this case the arguments are unchanged. + * + * The accuracy of the conversion is about 5nm. + * + * UTM eastings are allowed to be in the range [0km, 1000km], northings are + * allowed to be in in [0km, 9600km] for the northern hemisphere and in + * [900km, 10000km] for the southern hemisphere. However UTM northings + * can be continued across the equator. So the actual limits on the + * northings are [-9100km, 9600km] for the "northern" hemisphere and + * [900km, 19600km] for the "southern" hemisphere. + * + * UPS eastings and northings are allowed to be in the range [1200km, + * 2800km] in the northern hemisphere and in [700km, 3300km] in the + * southern hemisphere. + * + * These ranges are 100km larger than allowed for the conversions to MGRS. + * (100km is the maximum extra padding consistent with eastings remaining + * non-negative.) This allows generous overlaps between zones and UTM and + * UPS. If \e mgrslimits = true, then all the ranges are shrunk by 100km + * so that they agree with the stricter MGRS ranges. No checks are + * performed besides these (e.g., to limit the distance outside the + * standard zone boundaries). + **********************************************************************/ + static void Reverse(int zone, bool northp, real x, real y, + real& lat, real& lon, real& gamma, real& k, + bool mgrslimits = false); + + /** + * UTMUPS::Forward without returning convergence and scale. + **********************************************************************/ + static void Forward(real lat, real lon, + int& zone, bool& northp, real& x, real& y, + int setzone = STANDARD, bool mgrslimits = false) { + real gamma, k; + Forward(lat, lon, zone, northp, x, y, gamma, k, setzone, mgrslimits); + } + + /** + * UTMUPS::Reverse without returning convergence and scale. + **********************************************************************/ + static void Reverse(int zone, bool northp, real x, real y, + real& lat, real& lon, bool mgrslimits = false) { + real gamma, k; + Reverse(zone, northp, x, y, lat, lon, gamma, k, mgrslimits); + } + + /** + * Transfer UTM/UPS coordinated from one zone to another. + * + * @param[in] zonein the UTM zone for \e xin and \e yin (or zero for UPS). + * @param[in] northpin hemisphere for \e xin and \e yin (true means north, + * false means south). + * @param[in] xin easting of point (meters) in \e zonein. + * @param[in] yin northing of point (meters) in \e zonein. + * @param[in] zoneout the requested UTM zone for \e xout and \e yout (or + * zero for UPS). + * @param[in] northpout hemisphere for \e xout output and \e yout. + * @param[out] xout easting of point (meters) in \e zoneout. + * @param[out] yout northing of point (meters) in \e zoneout. + * @param[out] zone the actual UTM zone for \e xout and \e yout (or zero + * for UPS); this equals \e zoneout if \e zoneout ≥ 0. + * @exception GeographicErr if \e zonein is out of range (see below). + * @exception GeographicErr if \e zoneout is out of range (see below). + * @exception GeographicErr if \e xin or \e yin fall outside their allowed + * ranges (see UTMUPS::Reverse). + * @exception GeographicErr if \e xout or \e yout fall outside their + * allowed ranges (see UTMUPS::Reverse). + * + * \e zonein must be in the range [UTMUPS::MINZONE, UTMUPS::MAXZONE] = [0, + * 60] with \e zonein = UTMUPS::UPS, 0, indicating UPS. \e zonein may + * also be UTMUPS::INVALID. + * + * \e zoneout must be in the range [UTMUPS::MINPSEUDOZONE, UTMUPS::MAXZONE] + * = [-4, 60]. If \e zoneout < UTMUPS::MINZONE then the rules give in + * the documentation of UTMUPS::zonespec are applied, and \e zone is set to + * the actual zone used for output. + * + * (\e xout, \e yout) can overlap with (\e xin, \e yin). + **********************************************************************/ + static void Transfer(int zonein, bool northpin, real xin, real yin, + int zoneout, bool northpout, real& xout, real& yout, + int& zone); + + /** + * Decode a UTM/UPS zone string. + * + * @param[in] zonestr string representation of zone and hemisphere. + * @param[out] zone the UTM zone (zero means UPS). + * @param[out] northp hemisphere (true means north, false means south). + * @exception GeographicErr if \e zonestr is malformed. + * + * For UTM, \e zonestr has the form of a zone number in the range + * [UTMUPS::MINUTMZONE, UTMUPS::MAXUTMZONE] = [1, 60] followed by a + * hemisphere letter, n or s (or "north" or "south" spelled out). For UPS, + * it consists just of the hemisphere letter (or the spelled out + * hemisphere). The returned value of \e zone is UTMUPS::UPS = 0 for UPS. + * Note well that "38s" indicates the southern hemisphere of zone 38 and + * not latitude band S, 32° ≤ \e lat < 40°. n, 01s, 2n, 38s, + * south, 3north are legal. 0n, 001s, +3n, 61n, 38P are illegal. INV is a + * special value for which the returned value of \e is UTMUPS::INVALID. + **********************************************************************/ + static void DecodeZone(const std::string& zonestr, + int& zone, bool& northp); + + /** + * Encode a UTM/UPS zone string. + * + * @param[in] zone the UTM zone (zero means UPS). + * @param[in] northp hemisphere (true means north, false means south). + * @param[in] abbrev if true (the default) use abbreviated (n/s) notation + * for hemisphere; otherwise spell out the hemisphere (north/south) + * @exception GeographicErr if \e zone is out of range (see below). + * @exception std::bad_alloc if memoy for the string can't be allocated. + * @return string representation of zone and hemisphere. + * + * \e zone must be in the range [UTMUPS::MINZONE, UTMUPS::MAXZONE] = [0, + * 60] with \e zone = UTMUPS::UPS, 0, indicating UPS (but the resulting + * string does not contain "0"). \e zone may also be UTMUPS::INVALID, in + * which case the returned string is "inv". This reverses + * UTMUPS::DecodeZone. + **********************************************************************/ + static std::string EncodeZone(int zone, bool northp, bool abbrev = true); + + /** + * Decode EPSG. + * + * @param[in] epsg the EPSG code. + * @param[out] zone the UTM zone (zero means UPS). + * @param[out] northp hemisphere (true means north, false means south). + * + * EPSG (European Petroleum Survery Group) codes are a way to refer to many + * different projections. DecodeEPSG decodes those referring to UTM or UPS + * projections for the WGS84 ellipsoid. If the code does not refer to one + * of these projections, \e zone is set to UTMUPS::INVALID. See + * https://www.spatialreference.org/ref/epsg/ + **********************************************************************/ + static void DecodeEPSG(int epsg, int& zone, bool& northp); + + /** + * Encode zone as EPSG. + * + * @param[in] zone the UTM zone (zero means UPS). + * @param[in] northp hemisphere (true means north, false means south). + * @return EPSG code (or -1 if \e zone is not in the range + * [UTMUPS::MINZONE, UTMUPS::MAXZONE] = [0, 60]) + * + * Convert \e zone and \e northp to the corresponding EPSG (European + * Petroleum Survery Group) codes + **********************************************************************/ + static int EncodeEPSG(int zone, bool northp); + + /** + * @return shift (meters) necessary to align north and south halves of a + * UTM zone (107). + **********************************************************************/ + static Math::real UTMShift(); + + /** \name Inspector functions + **********************************************************************/ + ///@{ + /** + * @return \e a the equatorial radius of the WGS84 ellipsoid (meters). + * + * (The WGS84 value is returned because the UTM and UPS projections are + * based on this ellipsoid.) + **********************************************************************/ + static Math::real EquatorialRadius() + { return Constants::WGS84_a(); } + + /** + * @return \e f the flattening of the WGS84 ellipsoid. + * + * (The WGS84 value is returned because the UTM and UPS projections are + * based on this ellipsoid.) + **********************************************************************/ + static Math::real Flattening() + { return Constants::WGS84_f(); } + + /** + * \deprecated An old name for EquatorialRadius(). + **********************************************************************/ + // GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()") + static Math::real MajorRadius() { return EquatorialRadius(); } + ///@} + + }; + +} // namespace GeographicLib + +#endif // GEOGRAPHICLIB_UTMUPS_HPP diff --git a/src/Geo/Utility.cpp b/src/Geo/Utility.cpp new file mode 100644 index 000000000..0877b6363 --- /dev/null +++ b/src/Geo/Utility.cpp @@ -0,0 +1,61 @@ +/** + * \file Utility.cpp + * \brief Implementation for GeographicLib::Utility class + * + * Copyright (c) Charles Karney (2011) and licensed under + * the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#include +#include "Utility.hpp" + +#if defined(_MSC_VER) +// Squelch warnings about unsafe use of getenv +# pragma warning (disable: 4996) +#endif + +namespace GeographicLib { + + using namespace std; + + bool Utility::ParseLine(const std::string& line, + std::string& key, std::string& val) { + const char* spaces = " \t\n\v\f\r"; + string::size_type n0 = line.find_first_not_of(spaces); + if (n0 == string::npos) + return false; // Blank line + string::size_type n1 = line.find_first_of('#', n0); + if (n0 == n1) + return false; // Only a comment + val = line.substr(n0, n1 == string::npos ? n1 : n1 - n0); + n0 = val.find_first_of(spaces); + key = val.substr(0, n0); + if (n0 == string::npos) { + val = ""; + return true; + } + n0 = val.find_first_not_of(spaces, n0); + if (n0 == string::npos) { + val = ""; + return true; + } + n1 = val.find_last_not_of(spaces); + val = val.substr(n0, n1 + 1 - n0); + return true; + } + + int Utility::set_digits(int ndigits) { +#if GEOGRAPHICLIB_PRECISION == 5 + if (ndigits <= 0) { + char* digitenv = getenv("GEOGRAPHICLIB_DIGITS"); + if (digitenv) + ndigits = strtol(digitenv, NULL, 0); + if (ndigits <= 0) + ndigits = 256; + } +#endif + return Math::set_digits(ndigits); + } + +} // namespace GeographicLib diff --git a/src/Geo/Utility.h b/src/Geo/Utility.h new file mode 100644 index 000000000..329fedf40 --- /dev/null +++ b/src/Geo/Utility.h @@ -0,0 +1,704 @@ +/** + * \file Utility.hpp + * \brief Header for GeographicLib::Utility class + * + * Copyright (c) Charles Karney (2011-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#if !defined(GEOGRAPHICLIB_UTILITY_HPP) +#define GEOGRAPHICLIB_UTILITY_HPP 1 + +#include +#include +#include +#include +#include +#include +#include + +#if defined(_MSC_VER) +// Squelch warnings about constant conditional expressions and unsafe gmtime +# pragma warning (push) +# pragma warning (disable: 4127 4996) +#endif + +namespace GeographicLib { + + /** + * \brief Some utility routines for %GeographicLib + * + * Example of use: + * \include example-Utility.cpp + **********************************************************************/ + class GEOGRAPHICLIB_EXPORT Utility { + private: + static bool gregorian(int y, int m, int d) { + // The original cut over to the Gregorian calendar in Pope Gregory XIII's + // time had 1582-10-04 followed by 1582-10-15. Here we implement the + // switch over used by the English-speaking world where 1752-09-02 was + // followed by 1752-09-14. We also assume that the year always begins + // with January 1, whereas in reality it often was reckoned to begin in + // March. + return 100 * (100 * y + m) + d >= 17520914; // or 15821015 + } + static bool gregorian(int s) { + return s >= 639799; // 1752-09-14 + } + public: + + /** + * Convert a date to the day numbering sequentially starting with + * 0001-01-01 as day 1. + * + * @param[in] y the year (must be positive). + * @param[in] m the month, Jan = 1, etc. (must be positive). Default = 1. + * @param[in] d the day of the month (must be positive). Default = 1. + * @return the sequential day number. + **********************************************************************/ + static int day(int y, int m = 1, int d = 1) { + // Convert from date to sequential day and vice versa + // + // Here is some code to convert a date to sequential day and vice + // versa. The sequential day is numbered so that January 1, 1 AD is day 1 + // (a Saturday). So this is offset from the "Julian" day which starts the + // numbering with 4713 BC. + // + // This is inspired by a talk by John Conway at the John von Neumann + // National Supercomputer Center when he described his Doomsday algorithm + // for figuring the day of the week. The code avoids explicitly doing ifs + // (except for the decision of whether to use the Julian or Gregorian + // calendar). Instead the equivalent result is achieved using integer + // arithmetic. I got this idea from the routine for the day of the week + // in MACLisp (I believe that that routine was written by Guy Steele). + // + // There are three issues to take care of + // + // 1. the rules for leap years, + // 2. the inconvenient placement of leap days at the end of February, + // 3. the irregular pattern of month lengths. + // + // We deal with these as follows: + // + // 1. Leap years are given by simple rules which are straightforward to + // accommodate. + // + // 2. We simplify the calculations by moving January and February to the + // previous year. Here we internally number the months March–December, + // January, February as 0–9, 10, 11. + // + // 3. The pattern of month lengths from March through January is regular + // with a 5-month period—31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31. The + // 5-month period is 153 days long. Since February is now at the end of + // the year, we don't need to include its length in this part of the + // calculation. + bool greg = gregorian(y, m, d); + y += (m + 9) / 12 - 1; // Move Jan and Feb to previous year, + m = (m + 9) % 12; // making March month 0. + return + (1461 * y) / 4 // Julian years converted to days. Julian year is 365 + + // 1/4 = 1461/4 days. + // Gregorian leap year corrections. The 2 offset with respect to the + // Julian calendar synchronizes the vernal equinox with that at the + // time of the Council of Nicea (325 AD). + + (greg ? (y / 100) / 4 - (y / 100) + 2 : 0) + + (153 * m + 2) / 5 // The zero-based start of the m'th month + + d - 1 // The zero-based day + - 305; // The number of days between March 1 and December 31. + // This makes 0001-01-01 day 1 + } + + /** + * Convert a date to the day numbering sequentially starting with + * 0001-01-01 as day 1. + * + * @param[in] y the year (must be positive). + * @param[in] m the month, Jan = 1, etc. (must be positive). Default = 1. + * @param[in] d the day of the month (must be positive). Default = 1. + * @param[in] check whether to check the date. + * @exception GeographicErr if the date is invalid and \e check is true. + * @return the sequential day number. + **********************************************************************/ + static int day(int y, int m, int d, bool check) { + int s = day(y, m, d); + if (!check) + return s; + int y1, m1, d1; + date(s, y1, m1, d1); + if (!(s > 0 && y == y1 && m == m1 && d == d1)) + throw GeographicErr("Invalid date " + + str(y) + "-" + str(m) + "-" + str(d) + + (s > 0 ? "; use " + + str(y1) + "-" + str(m1) + "-" + str(d1) : + " before 0001-01-01")); + return s; + } + + /** + * Given a day (counting from 0001-01-01 as day 1), return the date. + * + * @param[in] s the sequential day number (must be positive) + * @param[out] y the year. + * @param[out] m the month, Jan = 1, etc. + * @param[out] d the day of the month. + **********************************************************************/ + static void date(int s, int& y, int& m, int& d) { + int c = 0; + bool greg = gregorian(s); + s += 305; // s = 0 on March 1, 1BC + if (greg) { + s -= 2; // The 2 day Gregorian offset + // Determine century with the Gregorian rules for leap years. The + // Gregorian year is 365 + 1/4 - 1/100 + 1/400 = 146097/400 days. + c = (4 * s + 3) / 146097; + s -= (c * 146097) / 4; // s = 0 at beginning of century + } + y = (4 * s + 3) / 1461; // Determine the year using Julian rules. + s -= (1461 * y) / 4; // s = 0 at start of year, i.e., March 1 + y += c * 100; // Assemble full year + m = (5 * s + 2) / 153; // Determine the month + s -= (153 * m + 2) / 5; // s = 0 at beginning of month + d = s + 1; // Determine day of month + y += (m + 2) / 12; // Move Jan and Feb back to original year + m = (m + 2) % 12 + 1; // Renumber the months so January = 1 + } + + /** + * Given a date as a string in the format yyyy, yyyy-mm, or yyyy-mm-dd, + * return the numeric values for the year, month, and day. No checking is + * done on these values. The string "now" is interpreted as the present + * date (in UTC). + * + * @param[in] s the date in string format. + * @param[out] y the year. + * @param[out] m the month, Jan = 1, etc. + * @param[out] d the day of the month. + * @exception GeographicErr is \e s is malformed. + **********************************************************************/ + static void date(const std::string& s, int& y, int& m, int& d) { + if (s == "now") { + std::time_t t = std::time(0); + struct tm* now = gmtime(&t); + y = now->tm_year + 1900; + m = now->tm_mon + 1; + d = now->tm_mday; + return; + } + int y1, m1 = 1, d1 = 1; + const char* digits = "0123456789"; + std::string::size_type p1 = s.find_first_not_of(digits); + if (p1 == std::string::npos) + y1 = val(s); + else if (s[p1] != '-') + throw GeographicErr("Delimiter not hyphen in date " + s); + else if (p1 == 0) + throw GeographicErr("Empty year field in date " + s); + else { + y1 = val(s.substr(0, p1)); + if (++p1 == s.size()) + throw GeographicErr("Empty month field in date " + s); + std::string::size_type p2 = s.find_first_not_of(digits, p1); + if (p2 == std::string::npos) + m1 = val(s.substr(p1)); + else if (s[p2] != '-') + throw GeographicErr("Delimiter not hyphen in date " + s); + else if (p2 == p1) + throw GeographicErr("Empty month field in date " + s); + else { + m1 = val(s.substr(p1, p2 - p1)); + if (++p2 == s.size()) + throw GeographicErr("Empty day field in date " + s); + d1 = val(s.substr(p2)); + } + } + y = y1; m = m1; d = d1; + } + + /** + * Given the date, return the day of the week. + * + * @param[in] y the year (must be positive). + * @param[in] m the month, Jan = 1, etc. (must be positive). + * @param[in] d the day of the month (must be positive). + * @return the day of the week with Sunday, Monday--Saturday = 0, + * 1--6. + **********************************************************************/ + static int dow(int y, int m, int d) { return dow(day(y, m, d)); } + + /** + * Given the sequential day, return the day of the week. + * + * @param[in] s the sequential day (must be positive). + * @return the day of the week with Sunday, Monday--Saturday = 0, + * 1--6. + **********************************************************************/ + static int dow(int s) { + return (s + 5) % 7; // The 5 offset makes day 1 (0001-01-01) a Saturday. + } + + /** + * Convert a string representing a date to a fractional year. + * + * @tparam T the type of the argument. + * @param[in] s the string to be converted. + * @exception GeographicErr if \e s can't be interpreted as a date. + * @return the fractional year. + * + * The string is first read as an ordinary number (e.g., 2010 or 2012.5); + * if this is successful, the value is returned. Otherwise the string + * should be of the form yyyy-mm or yyyy-mm-dd and this is converted to a + * number with 2010-01-01 giving 2010.0 and 2012-07-03 giving 2012.5. + **********************************************************************/ + template static T fractionalyear(const std::string& s) { + try { + return val(s); + } + catch (const std::exception&) {} + int y, m, d; + date(s, y, m, d); + int t = day(y, m, d, true); + return T(y) + T(t - day(y)) / T(day(y + 1) - day(y)); + } + + /** + * Convert a object of type T to a string. + * + * @tparam T the type of the argument. + * @param[in] x the value to be converted. + * @param[in] p the precision used (default −1). + * @exception std::bad_alloc if memory for the string can't be allocated. + * @return the string representation. + * + * If \e p ≥ 0, then the number fixed format is used with p bits of + * precision. With p < 0, there is no manipulation of the format. + **********************************************************************/ + template static std::string str(T x, int p = -1) { + std::ostringstream s; + if (p >= 0) s << std::fixed << std::setprecision(p); + s << x; return s.str(); + } + + /** + * Convert a Math::real object to a string. + * + * @param[in] x the value to be converted. + * @param[in] p the precision used (default −1). + * @exception std::bad_alloc if memory for the string can't be allocated. + * @return the string representation. + * + * If \e p ≥ 0, then the number fixed format is used with p bits of + * precision. With p < 0, there is no manipulation of the format. This is + * an overload of str which deals with inf and nan. + **********************************************************************/ + static std::string str(Math::real x, int p = -1) { + if (!Math::isfinite(x)) + return x < 0 ? std::string("-inf") : + (x > 0 ? std::string("inf") : std::string("nan")); + std::ostringstream s; +#if GEOGRAPHICLIB_PRECISION == 4 + // boost-quadmath treats precision == 0 as "use as many digits as + // necessary" (see https://svn.boost.org/trac/boost/ticket/10103), so... + using std::floor; using std::fmod; + if (p == 0) { + x += Math::real(0.5); + Math::real ix = floor(x); + // Implement the "round ties to even" rule + x = (ix == x && fmod(ix, Math::real(2)) == 1) ? ix - 1 : ix; + s << std::fixed << std::setprecision(1) << x; + std::string r(s.str()); + // strip off trailing ".0" + return r.substr(0, (std::max)(int(r.size()) - 2, 0)); + } +#endif + if (p >= 0) s << std::fixed << std::setprecision(p); + s << x; return s.str(); + } + + /** + * Trim the white space from the beginning and end of a string. + * + * @param[in] s the string to be trimmed + * @return the trimmed string + **********************************************************************/ + static std::string trim(const std::string& s) { + unsigned + beg = 0, + end = unsigned(s.size()); + while (beg < end && isspace(s[beg])) + ++beg; + while (beg < end && isspace(s[end - 1])) + --end; + return std::string(s, beg, end-beg); + } + + /** + * Convert a string to type T. + * + * @tparam T the type of the return value. + * @param[in] s the string to be converted. + * @exception GeographicErr is \e s is not readable as a T. + * @return object of type T. + * + * White space at the beginning and end of \e s is ignored. + * + * Special handling is provided for some types. + * + * If T is a floating point type, then inf and nan are recognized. + * + * If T is bool, then \e s should either be string a representing 0 (false) + * or 1 (true) or one of the strings + * - "false", "f", "nil", "no", "n", "off", or "" meaning false, + * - "true", "t", "yes", "y", or "on" meaning true; + * . + * case is ignored. + * + * If T is std::string, then \e s is returned (with the white space at the + * beginning and end removed). + **********************************************************************/ + template static T val(const std::string& s) { + // If T is bool, then the specialization val() defined below is + // used. + T x; + std::string errmsg, t(trim(s)); + do { // Executed once (provides the ability to break) + std::istringstream is(t); + if (!(is >> x)) { + errmsg = "Cannot decode " + t; + break; + } + int pos = int(is.tellg()); // Returns -1 at end of string? + if (!(pos < 0 || pos == int(t.size()))) { + errmsg = "Extra text " + t.substr(pos) + " at end of " + t; + break; + } + return x; + } while (false); + x = std::numeric_limits::is_integer ? 0 : nummatch(t); + if (x == 0) + throw GeographicErr(errmsg); + return x; + } + /** + * \deprecated An old name for val(s). + **********************************************************************/ + template + GEOGRAPHICLIB_DEPRECATED("Use Utility::val(s)") + static T num(const std::string& s) { + return val(s); + } + + /** + * Match "nan" and "inf" (and variants thereof) in a string. + * + * @tparam T the type of the return value (this should be a floating point + * type). + * @param[in] s the string to be matched. + * @return appropriate special value (±∞, nan) or 0 if none is + * found. + * + * White space is not allowed at the beginning or end of \e s. + **********************************************************************/ + template static T nummatch(const std::string& s) { + if (s.length() < 3) + return 0; + std::string t(s); + for (std::string::iterator p = t.begin(); p != t.end(); ++p) + *p = char(std::toupper(*p)); + for (size_t i = s.length(); i--;) + t[i] = char(std::toupper(s[i])); + int sign = t[0] == '-' ? -1 : 1; + std::string::size_type p0 = t[0] == '-' || t[0] == '+' ? 1 : 0; + std::string::size_type p1 = t.find_last_not_of('0'); + if (p1 == std::string::npos || p1 + 1 < p0 + 3) + return 0; + // Strip off sign and trailing 0s + t = t.substr(p0, p1 + 1 - p0); // Length at least 3 + if (t == "NAN" || t == "1.#QNAN" || t == "1.#SNAN" || t == "1.#IND" || + t == "1.#R") + return Math::NaN(); + else if (t == "INF" || t == "1.#INF") + return sign * Math::infinity(); + return 0; + } + + /** + * Read a simple fraction, e.g., 3/4, from a string to an object of type T. + * + * @tparam T the type of the return value. + * @param[in] s the string to be converted. + * @exception GeographicErr is \e s is not readable as a fraction of type + * T. + * @return object of type T + * + * \note The msys shell under Windows converts arguments which look + * like pathnames into their Windows equivalents. As a result the argument + * "-1/300" gets mangled into something unrecognizable. A workaround is to + * use a floating point number in the numerator, i.e., "-1.0/300". + **********************************************************************/ + template static T fract(const std::string& s) { + std::string::size_type delim = s.find('/'); + return + !(delim != std::string::npos && delim >= 1 && delim + 2 <= s.size()) ? + val(s) : + // delim in [1, size() - 2] + val(s.substr(0, delim)) / val(s.substr(delim + 1)); + } + + /** + * Lookup up a character in a string. + * + * @param[in] s the string to be searched. + * @param[in] c the character to look for. + * @return the index of the first occurrence character in the string or + * −1 is the character is not present. + * + * \e c is converted to upper case before search \e s. Therefore, it is + * intended that \e s should not contain any lower case letters. + **********************************************************************/ + static int lookup(const std::string& s, char c) { + std::string::size_type r = s.find(char(toupper(c))); + return r == std::string::npos ? -1 : int(r); + } + + /** + * Lookup up a character in a char*. + * + * @param[in] s the char* string to be searched. + * @param[in] c the character to look for. + * @return the index of the first occurrence character in the string or + * −1 is the character is not present. + * + * \e c is converted to upper case before search \e s. Therefore, it is + * intended that \e s should not contain any lower case letters. + **********************************************************************/ + static int lookup(const char* s, char c) { + const char* p = std::strchr(s, toupper(c)); + return p != NULL ? int(p - s) : -1; + } + + /** + * Read data of type ExtT from a binary stream to an array of type IntT. + * The data in the file is in (bigendp ? big : little)-endian format. + * + * @tparam ExtT the type of the objects in the binary stream (external). + * @tparam IntT the type of the objects in the array (internal). + * @tparam bigendp true if the external storage format is big-endian. + * @param[in] str the input stream containing the data of type ExtT + * (external). + * @param[out] array the output array of type IntT (internal). + * @param[in] num the size of the array. + * @exception GeographicErr if the data cannot be read. + **********************************************************************/ + template + static void readarray(std::istream& str, IntT array[], size_t num) { +#if GEOGRAPHICLIB_PRECISION < 4 + if (sizeof(IntT) == sizeof(ExtT) && + std::numeric_limits::is_integer == + std::numeric_limits::is_integer) + { + // Data is compatible (aside from the issue of endian-ness). + str.read(reinterpret_cast(array), num * sizeof(ExtT)); + if (!str.good()) + throw GeographicErr("Failure reading data"); + if (bigendp != Math::bigendian) { // endian mismatch -> swap bytes + for (size_t i = num; i--;) + array[i] = Math::swab(array[i]); + } + } + else +#endif + { + const int bufsize = 1024; // read this many values at a time + ExtT buffer[bufsize]; // temporary buffer + int k = int(num); // data values left to read + int i = 0; // index into output array + while (k) { + int n = (std::min)(k, bufsize); + str.read(reinterpret_cast(buffer), n * sizeof(ExtT)); + if (!str.good()) + throw GeographicErr("Failure reading data"); + for (int j = 0; j < n; ++j) + // fix endian-ness and cast to IntT + array[i++] = IntT(bigendp == Math::bigendian ? buffer[j] : + Math::swab(buffer[j])); + k -= n; + } + } + return; + } + + /** + * Read data of type ExtT from a binary stream to a vector array of type + * IntT. The data in the file is in (bigendp ? big : little)-endian + * format. + * + * @tparam ExtT the type of the objects in the binary stream (external). + * @tparam IntT the type of the objects in the array (internal). + * @tparam bigendp true if the external storage format is big-endian. + * @param[in] str the input stream containing the data of type ExtT + * (external). + * @param[out] array the output vector of type IntT (internal). + * @exception GeographicErr if the data cannot be read. + **********************************************************************/ + template + static void readarray(std::istream& str, std::vector& array) { + if (array.size() > 0) + readarray(str, &array[0], array.size()); + } + + /** + * Write data in an array of type IntT as type ExtT to a binary stream. + * The data in the file is in (bigendp ? big : little)-endian format. + * + * @tparam ExtT the type of the objects in the binary stream (external). + * @tparam IntT the type of the objects in the array (internal). + * @tparam bigendp true if the external storage format is big-endian. + * @param[out] str the output stream for the data of type ExtT (external). + * @param[in] array the input array of type IntT (internal). + * @param[in] num the size of the array. + * @exception GeographicErr if the data cannot be written. + **********************************************************************/ + template + static void writearray(std::ostream& str, const IntT array[], size_t num) + { +#if GEOGRAPHICLIB_PRECISION < 4 + if (sizeof(IntT) == sizeof(ExtT) && + std::numeric_limits::is_integer == + std::numeric_limits::is_integer && + bigendp == Math::bigendian) + { + // Data is compatible (including endian-ness). + str.write(reinterpret_cast(array), num * sizeof(ExtT)); + if (!str.good()) + throw GeographicErr("Failure writing data"); + } + else +#endif + { + const int bufsize = 1024; // write this many values at a time + ExtT buffer[bufsize]; // temporary buffer + int k = int(num); // data values left to write + int i = 0; // index into output array + while (k) { + int n = (std::min)(k, bufsize); + for (int j = 0; j < n; ++j) + // cast to ExtT and fix endian-ness + buffer[j] = bigendp == Math::bigendian ? ExtT(array[i++]) : + Math::swab(ExtT(array[i++])); + str.write(reinterpret_cast(buffer), n * sizeof(ExtT)); + if (!str.good()) + throw GeographicErr("Failure writing data"); + k -= n; + } + } + return; + } + + /** + * Write data in an array of type IntT as type ExtT to a binary stream. + * The data in the file is in (bigendp ? big : little)-endian format. + * + * @tparam ExtT the type of the objects in the binary stream (external). + * @tparam IntT the type of the objects in the array (internal). + * @tparam bigendp true if the external storage format is big-endian. + * @param[out] str the output stream for the data of type ExtT (external). + * @param[in] array the input vector of type IntT (internal). + * @exception GeographicErr if the data cannot be written. + **********************************************************************/ + template + static void writearray(std::ostream& str, std::vector& array) { + if (array.size() > 0) + writearray(str, &array[0], array.size()); + } + + /** + * Parse a KEY VALUE line. + * + * @param[in] line the input line. + * @param[out] key the key. + * @param[out] val the value. + * @exception std::bad_alloc if memory for the internal strings can't be + * allocated. + * @return whether a key was found. + * + * A # character and everything after it are discarded. If the result is + * just white space, the routine returns false (and \e key and \e val are + * not set). Otherwise the first token is taken to be the key and the rest + * of the line (trimmed of leading and trailing white space) is the value. + **********************************************************************/ + static bool ParseLine(const std::string& line, + std::string& key, std::string& val); + + /** + * Set the binary precision of a real number. + * + * @param[in] ndigits the number of bits of precision. If ndigits is 0 + * (the default), then determine the precision from the environment + * variable GEOGRAPHICLIB_DIGITS. If this is undefined, use ndigits = + * 256 (i.e., about 77 decimal digits). + * @return the resulting number of bits of precision. + * + * This only has an effect when GEOGRAPHICLIB_PRECISION = 5. The + * precision should only be set once and before calls to any other + * GeographicLib functions. (Several functions, for example Math::pi(), + * cache the return value in a static local variable. The precision needs + * to be set before a call to any such functions.) In multi-threaded + * applications, it is necessary also to set the precision in each thread + * (see the example GeoidToGTX.cpp). + **********************************************************************/ + static int set_digits(int ndigits = 0); + + }; + + /** + * The specialization of Utility::val() for strings. + **********************************************************************/ + template<> inline std::string Utility::val(const std::string& s) + { return trim(s); } + + /** + * The specialization of Utility::val() for bools. + **********************************************************************/ + template<> inline bool Utility::val(const std::string& s) { + std::string t(trim(s)); + if (t.empty()) return false; + bool x; + std::istringstream is(t); + if (is >> x) { + int pos = int(is.tellg()); // Returns -1 at end of string? + if (!(pos < 0 || pos == int(t.size()))) + throw GeographicErr("Extra text " + t.substr(pos) + + " at end of " + t); + return x; + } + for (std::string::iterator p = t.begin(); p != t.end(); ++p) + *p = char(std::tolower(*p)); + switch (t[0]) { // already checked that t isn't empty + case 'f': + if (t == "f" || t == "false") return false; + break; + case 'n': + if (t == "n" || t == "nil" || t == "no") return false; + break; + case 'o': + if (t == "off") return false; + else if (t == "on") return true; + break; + case 't': + if (t == "t" || t == "true") return true; + break; + case 'y': + if (t == "y" || t == "yes") return true; + break; + } + throw GeographicErr("Cannot decode " + t + " as a bool"); + } + +} // namespace GeographicLib + +#if defined(_MSC_VER) +# pragma warning (pop) +#endif + +#endif // GEOGRAPHICLIB_UTILITY_HPP diff --git a/src/Geo/Utility.hpp b/src/Geo/Utility.hpp new file mode 100644 index 000000000..c0288ec2e --- /dev/null +++ b/src/Geo/Utility.hpp @@ -0,0 +1,704 @@ +/** + * \file Utility.hpp + * \brief Header for GeographicLib::Utility class + * + * Copyright (c) Charles Karney (2011-2019) and licensed + * under the MIT/X11 License. For more information, see + * https://geographiclib.sourceforge.io/ + **********************************************************************/ + +#if !defined(GEOGRAPHICLIB_UTILITY_HPP) +#define GEOGRAPHICLIB_UTILITY_HPP 1 + +#include "Constants.hpp" +#include +#include +#include +#include +#include +#include + +#if defined(_MSC_VER) +// Squelch warnings about constant conditional expressions and unsafe gmtime +# pragma warning (push) +# pragma warning (disable: 4127 4996) +#endif + +namespace GeographicLib { + + /** + * \brief Some utility routines for %GeographicLib + * + * Example of use: + * \include example-Utility.cpp + **********************************************************************/ + class GEOGRAPHICLIB_EXPORT Utility { + private: + static bool gregorian(int y, int m, int d) { + // The original cut over to the Gregorian calendar in Pope Gregory XIII's + // time had 1582-10-04 followed by 1582-10-15. Here we implement the + // switch over used by the English-speaking world where 1752-09-02 was + // followed by 1752-09-14. We also assume that the year always begins + // with January 1, whereas in reality it often was reckoned to begin in + // March. + return 100 * (100 * y + m) + d >= 17520914; // or 15821015 + } + static bool gregorian(int s) { + return s >= 639799; // 1752-09-14 + } + public: + + /** + * Convert a date to the day numbering sequentially starting with + * 0001-01-01 as day 1. + * + * @param[in] y the year (must be positive). + * @param[in] m the month, Jan = 1, etc. (must be positive). Default = 1. + * @param[in] d the day of the month (must be positive). Default = 1. + * @return the sequential day number. + **********************************************************************/ + static int day(int y, int m = 1, int d = 1) { + // Convert from date to sequential day and vice versa + // + // Here is some code to convert a date to sequential day and vice + // versa. The sequential day is numbered so that January 1, 1 AD is day 1 + // (a Saturday). So this is offset from the "Julian" day which starts the + // numbering with 4713 BC. + // + // This is inspired by a talk by John Conway at the John von Neumann + // National Supercomputer Center when he described his Doomsday algorithm + // for figuring the day of the week. The code avoids explicitly doing ifs + // (except for the decision of whether to use the Julian or Gregorian + // calendar). Instead the equivalent result is achieved using integer + // arithmetic. I got this idea from the routine for the day of the week + // in MACLisp (I believe that that routine was written by Guy Steele). + // + // There are three issues to take care of + // + // 1. the rules for leap years, + // 2. the inconvenient placement of leap days at the end of February, + // 3. the irregular pattern of month lengths. + // + // We deal with these as follows: + // + // 1. Leap years are given by simple rules which are straightforward to + // accommodate. + // + // 2. We simplify the calculations by moving January and February to the + // previous year. Here we internally number the months March–December, + // January, February as 0–9, 10, 11. + // + // 3. The pattern of month lengths from March through January is regular + // with a 5-month period—31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31. The + // 5-month period is 153 days long. Since February is now at the end of + // the year, we don't need to include its length in this part of the + // calculation. + bool greg = gregorian(y, m, d); + y += (m + 9) / 12 - 1; // Move Jan and Feb to previous year, + m = (m + 9) % 12; // making March month 0. + return + (1461 * y) / 4 // Julian years converted to days. Julian year is 365 + + // 1/4 = 1461/4 days. + // Gregorian leap year corrections. The 2 offset with respect to the + // Julian calendar synchronizes the vernal equinox with that at the + // time of the Council of Nicea (325 AD). + + (greg ? (y / 100) / 4 - (y / 100) + 2 : 0) + + (153 * m + 2) / 5 // The zero-based start of the m'th month + + d - 1 // The zero-based day + - 305; // The number of days between March 1 and December 31. + // This makes 0001-01-01 day 1 + } + + /** + * Convert a date to the day numbering sequentially starting with + * 0001-01-01 as day 1. + * + * @param[in] y the year (must be positive). + * @param[in] m the month, Jan = 1, etc. (must be positive). Default = 1. + * @param[in] d the day of the month (must be positive). Default = 1. + * @param[in] check whether to check the date. + * @exception GeographicErr if the date is invalid and \e check is true. + * @return the sequential day number. + **********************************************************************/ + static int day(int y, int m, int d, bool check) { + int s = day(y, m, d); + if (!check) + return s; + int y1, m1, d1; + date(s, y1, m1, d1); + if (!(s > 0 && y == y1 && m == m1 && d == d1)) + throw GeographicErr("Invalid date " + + str(y) + "-" + str(m) + "-" + str(d) + + (s > 0 ? "; use " + + str(y1) + "-" + str(m1) + "-" + str(d1) : + " before 0001-01-01")); + return s; + } + + /** + * Given a day (counting from 0001-01-01 as day 1), return the date. + * + * @param[in] s the sequential day number (must be positive) + * @param[out] y the year. + * @param[out] m the month, Jan = 1, etc. + * @param[out] d the day of the month. + **********************************************************************/ + static void date(int s, int& y, int& m, int& d) { + int c = 0; + bool greg = gregorian(s); + s += 305; // s = 0 on March 1, 1BC + if (greg) { + s -= 2; // The 2 day Gregorian offset + // Determine century with the Gregorian rules for leap years. The + // Gregorian year is 365 + 1/4 - 1/100 + 1/400 = 146097/400 days. + c = (4 * s + 3) / 146097; + s -= (c * 146097) / 4; // s = 0 at beginning of century + } + y = (4 * s + 3) / 1461; // Determine the year using Julian rules. + s -= (1461 * y) / 4; // s = 0 at start of year, i.e., March 1 + y += c * 100; // Assemble full year + m = (5 * s + 2) / 153; // Determine the month + s -= (153 * m + 2) / 5; // s = 0 at beginning of month + d = s + 1; // Determine day of month + y += (m + 2) / 12; // Move Jan and Feb back to original year + m = (m + 2) % 12 + 1; // Renumber the months so January = 1 + } + + /** + * Given a date as a string in the format yyyy, yyyy-mm, or yyyy-mm-dd, + * return the numeric values for the year, month, and day. No checking is + * done on these values. The string "now" is interpreted as the present + * date (in UTC). + * + * @param[in] s the date in string format. + * @param[out] y the year. + * @param[out] m the month, Jan = 1, etc. + * @param[out] d the day of the month. + * @exception GeographicErr is \e s is malformed. + **********************************************************************/ + static void date(const std::string& s, int& y, int& m, int& d) { + if (s == "now") { + std::time_t t = std::time(0); + struct tm* now = gmtime(&t); + y = now->tm_year + 1900; + m = now->tm_mon + 1; + d = now->tm_mday; + return; + } + int y1, m1 = 1, d1 = 1; + const char* digits = "0123456789"; + std::string::size_type p1 = s.find_first_not_of(digits); + if (p1 == std::string::npos) + y1 = val(s); + else if (s[p1] != '-') + throw GeographicErr("Delimiter not hyphen in date " + s); + else if (p1 == 0) + throw GeographicErr("Empty year field in date " + s); + else { + y1 = val(s.substr(0, p1)); + if (++p1 == s.size()) + throw GeographicErr("Empty month field in date " + s); + std::string::size_type p2 = s.find_first_not_of(digits, p1); + if (p2 == std::string::npos) + m1 = val(s.substr(p1)); + else if (s[p2] != '-') + throw GeographicErr("Delimiter not hyphen in date " + s); + else if (p2 == p1) + throw GeographicErr("Empty month field in date " + s); + else { + m1 = val(s.substr(p1, p2 - p1)); + if (++p2 == s.size()) + throw GeographicErr("Empty day field in date " + s); + d1 = val(s.substr(p2)); + } + } + y = y1; m = m1; d = d1; + } + + /** + * Given the date, return the day of the week. + * + * @param[in] y the year (must be positive). + * @param[in] m the month, Jan = 1, etc. (must be positive). + * @param[in] d the day of the month (must be positive). + * @return the day of the week with Sunday, Monday--Saturday = 0, + * 1--6. + **********************************************************************/ + static int dow(int y, int m, int d) { return dow(day(y, m, d)); } + + /** + * Given the sequential day, return the day of the week. + * + * @param[in] s the sequential day (must be positive). + * @return the day of the week with Sunday, Monday--Saturday = 0, + * 1--6. + **********************************************************************/ + static int dow(int s) { + return (s + 5) % 7; // The 5 offset makes day 1 (0001-01-01) a Saturday. + } + + /** + * Convert a string representing a date to a fractional year. + * + * @tparam T the type of the argument. + * @param[in] s the string to be converted. + * @exception GeographicErr if \e s can't be interpreted as a date. + * @return the fractional year. + * + * The string is first read as an ordinary number (e.g., 2010 or 2012.5); + * if this is successful, the value is returned. Otherwise the string + * should be of the form yyyy-mm or yyyy-mm-dd and this is converted to a + * number with 2010-01-01 giving 2010.0 and 2012-07-03 giving 2012.5. + **********************************************************************/ + template static T fractionalyear(const std::string& s) { + try { + return val(s); + } + catch (const std::exception&) {} + int y, m, d; + date(s, y, m, d); + int t = day(y, m, d, true); + return T(y) + T(t - day(y)) / T(day(y + 1) - day(y)); + } + + /** + * Convert a object of type T to a string. + * + * @tparam T the type of the argument. + * @param[in] x the value to be converted. + * @param[in] p the precision used (default −1). + * @exception std::bad_alloc if memory for the string can't be allocated. + * @return the string representation. + * + * If \e p ≥ 0, then the number fixed format is used with p bits of + * precision. With p < 0, there is no manipulation of the format. + **********************************************************************/ + template static std::string str(T x, int p = -1) { + std::ostringstream s; + if (p >= 0) s << std::fixed << std::setprecision(p); + s << x; return s.str(); + } + + /** + * Convert a Math::real object to a string. + * + * @param[in] x the value to be converted. + * @param[in] p the precision used (default −1). + * @exception std::bad_alloc if memory for the string can't be allocated. + * @return the string representation. + * + * If \e p ≥ 0, then the number fixed format is used with p bits of + * precision. With p < 0, there is no manipulation of the format. This is + * an overload of str which deals with inf and nan. + **********************************************************************/ + static std::string str(Math::real x, int p = -1) { + if (!Math::isfinite(x)) + return x < 0 ? std::string("-inf") : + (x > 0 ? std::string("inf") : std::string("nan")); + std::ostringstream s; +#if GEOGRAPHICLIB_PRECISION == 4 + // boost-quadmath treats precision == 0 as "use as many digits as + // necessary" (see https://svn.boost.org/trac/boost/ticket/10103), so... + using std::floor; using std::fmod; + if (p == 0) { + x += Math::real(0.5); + Math::real ix = floor(x); + // Implement the "round ties to even" rule + x = (ix == x && fmod(ix, Math::real(2)) == 1) ? ix - 1 : ix; + s << std::fixed << std::setprecision(1) << x; + std::string r(s.str()); + // strip off trailing ".0" + return r.substr(0, (std::max)(int(r.size()) - 2, 0)); + } +#endif + if (p >= 0) s << std::fixed << std::setprecision(p); + s << x; return s.str(); + } + + /** + * Trim the white space from the beginning and end of a string. + * + * @param[in] s the string to be trimmed + * @return the trimmed string + **********************************************************************/ + static std::string trim(const std::string& s) { + unsigned + beg = 0, + end = unsigned(s.size()); + while (beg < end && isspace(s[beg])) + ++beg; + while (beg < end && isspace(s[end - 1])) + --end; + return std::string(s, beg, end-beg); + } + + /** + * Convert a string to type T. + * + * @tparam T the type of the return value. + * @param[in] s the string to be converted. + * @exception GeographicErr is \e s is not readable as a T. + * @return object of type T. + * + * White space at the beginning and end of \e s is ignored. + * + * Special handling is provided for some types. + * + * If T is a floating point type, then inf and nan are recognized. + * + * If T is bool, then \e s should either be string a representing 0 (false) + * or 1 (true) or one of the strings + * - "false", "f", "nil", "no", "n", "off", or "" meaning false, + * - "true", "t", "yes", "y", or "on" meaning true; + * . + * case is ignored. + * + * If T is std::string, then \e s is returned (with the white space at the + * beginning and end removed). + **********************************************************************/ + template static T val(const std::string& s) { + // If T is bool, then the specialization val() defined below is + // used. + T x; + std::string errmsg, t(trim(s)); + do { // Executed once (provides the ability to break) + std::istringstream is(t); + if (!(is >> x)) { + errmsg = "Cannot decode " + t; + break; + } + int pos = int(is.tellg()); // Returns -1 at end of string? + if (!(pos < 0 || pos == int(t.size()))) { + errmsg = "Extra text " + t.substr(pos) + " at end of " + t; + break; + } + return x; + } while (false); + x = std::numeric_limits::is_integer ? 0 : nummatch(t); + if (x == 0) + throw GeographicErr(errmsg); + return x; + } + /** + * \deprecated An old name for val(s). + **********************************************************************/ + template + GEOGRAPHICLIB_DEPRECATED("Use Utility::val(s)") + static T num(const std::string& s) { + return val(s); + } + + /** + * Match "nan" and "inf" (and variants thereof) in a string. + * + * @tparam T the type of the return value (this should be a floating point + * type). + * @param[in] s the string to be matched. + * @return appropriate special value (±∞, nan) or 0 if none is + * found. + * + * White space is not allowed at the beginning or end of \e s. + **********************************************************************/ + template static T nummatch(const std::string& s) { + if (s.length() < 3) + return 0; + std::string t(s); + for (std::string::iterator p = t.begin(); p != t.end(); ++p) + *p = char(std::toupper(*p)); + for (size_t i = s.length(); i--;) + t[i] = char(std::toupper(s[i])); + int sign = t[0] == '-' ? -1 : 1; + std::string::size_type p0 = t[0] == '-' || t[0] == '+' ? 1 : 0; + std::string::size_type p1 = t.find_last_not_of('0'); + if (p1 == std::string::npos || p1 + 1 < p0 + 3) + return 0; + // Strip off sign and trailing 0s + t = t.substr(p0, p1 + 1 - p0); // Length at least 3 + if (t == "NAN" || t == "1.#QNAN" || t == "1.#SNAN" || t == "1.#IND" || + t == "1.#R") + return Math::NaN(); + else if (t == "INF" || t == "1.#INF") + return sign * Math::infinity(); + return 0; + } + + /** + * Read a simple fraction, e.g., 3/4, from a string to an object of type T. + * + * @tparam T the type of the return value. + * @param[in] s the string to be converted. + * @exception GeographicErr is \e s is not readable as a fraction of type + * T. + * @return object of type T + * + * \note The msys shell under Windows converts arguments which look + * like pathnames into their Windows equivalents. As a result the argument + * "-1/300" gets mangled into something unrecognizable. A workaround is to + * use a floating point number in the numerator, i.e., "-1.0/300". + **********************************************************************/ + template static T fract(const std::string& s) { + std::string::size_type delim = s.find('/'); + return + !(delim != std::string::npos && delim >= 1 && delim + 2 <= s.size()) ? + val(s) : + // delim in [1, size() - 2] + val(s.substr(0, delim)) / val(s.substr(delim + 1)); + } + + /** + * Lookup up a character in a string. + * + * @param[in] s the string to be searched. + * @param[in] c the character to look for. + * @return the index of the first occurrence character in the string or + * −1 is the character is not present. + * + * \e c is converted to upper case before search \e s. Therefore, it is + * intended that \e s should not contain any lower case letters. + **********************************************************************/ + static int lookup(const std::string& s, char c) { + std::string::size_type r = s.find(char(toupper(c))); + return r == std::string::npos ? -1 : int(r); + } + + /** + * Lookup up a character in a char*. + * + * @param[in] s the char* string to be searched. + * @param[in] c the character to look for. + * @return the index of the first occurrence character in the string or + * −1 is the character is not present. + * + * \e c is converted to upper case before search \e s. Therefore, it is + * intended that \e s should not contain any lower case letters. + **********************************************************************/ + static int lookup(const char* s, char c) { + const char* p = std::strchr(s, toupper(c)); + return p != NULL ? int(p - s) : -1; + } + + /** + * Read data of type ExtT from a binary stream to an array of type IntT. + * The data in the file is in (bigendp ? big : little)-endian format. + * + * @tparam ExtT the type of the objects in the binary stream (external). + * @tparam IntT the type of the objects in the array (internal). + * @tparam bigendp true if the external storage format is big-endian. + * @param[in] str the input stream containing the data of type ExtT + * (external). + * @param[out] array the output array of type IntT (internal). + * @param[in] num the size of the array. + * @exception GeographicErr if the data cannot be read. + **********************************************************************/ + template + static void readarray(std::istream& str, IntT array[], size_t num) { +#if GEOGRAPHICLIB_PRECISION < 4 + if (sizeof(IntT) == sizeof(ExtT) && + std::numeric_limits::is_integer == + std::numeric_limits::is_integer) + { + // Data is compatible (aside from the issue of endian-ness). + str.read(reinterpret_cast(array), num * sizeof(ExtT)); + if (!str.good()) + throw GeographicErr("Failure reading data"); + if (bigendp != Math::bigendian) { // endian mismatch -> swap bytes + for (size_t i = num; i--;) + array[i] = Math::swab(array[i]); + } + } + else +#endif + { + const int bufsize = 1024; // read this many values at a time + ExtT buffer[bufsize]; // temporary buffer + int k = int(num); // data values left to read + int i = 0; // index into output array + while (k) { + int n = (std::min)(k, bufsize); + str.read(reinterpret_cast(buffer), n * sizeof(ExtT)); + if (!str.good()) + throw GeographicErr("Failure reading data"); + for (int j = 0; j < n; ++j) + // fix endian-ness and cast to IntT + array[i++] = IntT(bigendp == Math::bigendian ? buffer[j] : + Math::swab(buffer[j])); + k -= n; + } + } + return; + } + + /** + * Read data of type ExtT from a binary stream to a vector array of type + * IntT. The data in the file is in (bigendp ? big : little)-endian + * format. + * + * @tparam ExtT the type of the objects in the binary stream (external). + * @tparam IntT the type of the objects in the array (internal). + * @tparam bigendp true if the external storage format is big-endian. + * @param[in] str the input stream containing the data of type ExtT + * (external). + * @param[out] array the output vector of type IntT (internal). + * @exception GeographicErr if the data cannot be read. + **********************************************************************/ + template + static void readarray(std::istream& str, std::vector& array) { + if (array.size() > 0) + readarray(str, &array[0], array.size()); + } + + /** + * Write data in an array of type IntT as type ExtT to a binary stream. + * The data in the file is in (bigendp ? big : little)-endian format. + * + * @tparam ExtT the type of the objects in the binary stream (external). + * @tparam IntT the type of the objects in the array (internal). + * @tparam bigendp true if the external storage format is big-endian. + * @param[out] str the output stream for the data of type ExtT (external). + * @param[in] array the input array of type IntT (internal). + * @param[in] num the size of the array. + * @exception GeographicErr if the data cannot be written. + **********************************************************************/ + template + static void writearray(std::ostream& str, const IntT array[], size_t num) + { +#if GEOGRAPHICLIB_PRECISION < 4 + if (sizeof(IntT) == sizeof(ExtT) && + std::numeric_limits::is_integer == + std::numeric_limits::is_integer && + bigendp == Math::bigendian) + { + // Data is compatible (including endian-ness). + str.write(reinterpret_cast(array), num * sizeof(ExtT)); + if (!str.good()) + throw GeographicErr("Failure writing data"); + } + else +#endif + { + const int bufsize = 1024; // write this many values at a time + ExtT buffer[bufsize]; // temporary buffer + int k = int(num); // data values left to write + int i = 0; // index into output array + while (k) { + int n = (std::min)(k, bufsize); + for (int j = 0; j < n; ++j) + // cast to ExtT and fix endian-ness + buffer[j] = bigendp == Math::bigendian ? ExtT(array[i++]) : + Math::swab(ExtT(array[i++])); + str.write(reinterpret_cast(buffer), n * sizeof(ExtT)); + if (!str.good()) + throw GeographicErr("Failure writing data"); + k -= n; + } + } + return; + } + + /** + * Write data in an array of type IntT as type ExtT to a binary stream. + * The data in the file is in (bigendp ? big : little)-endian format. + * + * @tparam ExtT the type of the objects in the binary stream (external). + * @tparam IntT the type of the objects in the array (internal). + * @tparam bigendp true if the external storage format is big-endian. + * @param[out] str the output stream for the data of type ExtT (external). + * @param[in] array the input vector of type IntT (internal). + * @exception GeographicErr if the data cannot be written. + **********************************************************************/ + template + static void writearray(std::ostream& str, std::vector& array) { + if (array.size() > 0) + writearray(str, &array[0], array.size()); + } + + /** + * Parse a KEY VALUE line. + * + * @param[in] line the input line. + * @param[out] key the key. + * @param[out] val the value. + * @exception std::bad_alloc if memory for the internal strings can't be + * allocated. + * @return whether a key was found. + * + * A # character and everything after it are discarded. If the result is + * just white space, the routine returns false (and \e key and \e val are + * not set). Otherwise the first token is taken to be the key and the rest + * of the line (trimmed of leading and trailing white space) is the value. + **********************************************************************/ + static bool ParseLine(const std::string& line, + std::string& key, std::string& val); + + /** + * Set the binary precision of a real number. + * + * @param[in] ndigits the number of bits of precision. If ndigits is 0 + * (the default), then determine the precision from the environment + * variable GEOGRAPHICLIB_DIGITS. If this is undefined, use ndigits = + * 256 (i.e., about 77 decimal digits). + * @return the resulting number of bits of precision. + * + * This only has an effect when GEOGRAPHICLIB_PRECISION = 5. The + * precision should only be set once and before calls to any other + * GeographicLib functions. (Several functions, for example Math::pi(), + * cache the return value in a static local variable. The precision needs + * to be set before a call to any such functions.) In multi-threaded + * applications, it is necessary also to set the precision in each thread + * (see the example GeoidToGTX.cpp). + **********************************************************************/ + static int set_digits(int ndigits = 0); + + }; + + /** + * The specialization of Utility::val() for strings. + **********************************************************************/ + template<> inline std::string Utility::val(const std::string& s) + { return trim(s); } + + /** + * The specialization of Utility::val() for bools. + **********************************************************************/ + template<> inline bool Utility::val(const std::string& s) { + std::string t(trim(s)); + if (t.empty()) return false; + bool x; + std::istringstream is(t); + if (is >> x) { + int pos = int(is.tellg()); // Returns -1 at end of string? + if (!(pos < 0 || pos == int(t.size()))) + throw GeographicErr("Extra text " + t.substr(pos) + + " at end of " + t); + return x; + } + for (std::string::iterator p = t.begin(); p != t.end(); ++p) + *p = char(std::tolower(*p)); + switch (t[0]) { // already checked that t isn't empty + case 'f': + if (t == "f" || t == "false") return false; + break; + case 'n': + if (t == "n" || t == "nil" || t == "no") return false; + break; + case 'o': + if (t == "off") return false; + else if (t == "on") return true; + break; + case 't': + if (t == "t" || t == "true") return true; + break; + case 'y': + if (t == "y" || t == "yes") return true; + break; + } + throw GeographicErr("Cannot decode " + t + " as a bool"); + } + +} // namespace GeographicLib + +#if defined(_MSC_VER) +# pragma warning (pop) +#endif + +#endif // GEOGRAPHICLIB_UTILITY_HPP -- 2.22.0