Skip to content

Q
qgroundcontrol
Commit 1ee854e1 authored by Matej Frančeškin's avatar Matej Frančeškin
Browse files
Options

MGRS - Added GeographicLib

parent afcc3fb6
Hide whitespace changes
Inline Side-by-side
Showing with 5768 additions and 4 deletions
qgroundcontrol.pro
View file @ 1ee854e1
......@@ -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 \
......
src/Geo/Constants.hpp 0 → 100644
View file @ 1ee854e1
/**
* \file Constants.hpp
* \brief Header for GeographicLib::Constants class
*
* Copyright (c) Charles Karney (2008-2019) <charles@karney.com> 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 <stdexcept>
#include <string>
/**
* \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<real>().
**********************************************************************/
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<typename T> static T WGS84_a()
{ return 6378137 * meter<T>(); }
/**
* A synonym for WGS84_a<real>().
**********************************************************************/
static real WGS84_a() { return WGS84_a<real>(); }
/**
* @tparam T the type of the returned value.
* @return the flattening of WGS84 ellipsoid (1/298.257223563).
**********************************************************************/
template<typename T> 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<real>().
**********************************************************************/
static real WGS84_f() { return WGS84_f<real>(); }
/**
* @tparam T the type of the returned value.
* @return the gravitational constant of the WGS84 ellipsoid, \e GM, in
* m<sup>3</sup> s<sup>&minus;2</sup>.
**********************************************************************/
template<typename T> static T WGS84_GM()
{ return T(3986004) * 100000000 + 41800000; }
/**
* A synonym for WGS84_GM<real>().
**********************************************************************/
static real WGS84_GM() { return WGS84_GM<real>(); }
/**
* @tparam T the type of the returned value.
* @return the angular velocity of the WGS84 ellipsoid, &omega;, in rad
* s<sup>&minus;1</sup>.
**********************************************************************/
template<typename T> static T WGS84_omega()
{ return 7292115 / (T(1000000) * 100000); }
/**
* A synonym for WGS84_omega<real>().
**********************************************************************/
static real WGS84_omega() { return WGS84_omega<real>(); }
/**
* @tparam T the type of the returned value.
* @return the equatorial radius of GRS80 ellipsoid, \e a, in m.
**********************************************************************/
template<typename T> static T GRS80_a()
{ return 6378137 * meter<T>(); }
/**
* A synonym for GRS80_a<real>().
**********************************************************************/
static real GRS80_a() { return GRS80_a<real>(); }
/**
* @tparam T the type of the returned value.
* @return the gravitational constant of the GRS80 ellipsoid, \e GM, in
* m<sup>3</sup> s<sup>&minus;2</sup>.
**********************************************************************/
template<typename T> static T GRS80_GM()
{ return T(3986005) * 100000000; }
/**
* A synonym for GRS80_GM<real>().
**********************************************************************/
static real GRS80_GM() { return GRS80_GM<real>(); }
/**
* @tparam T the type of the returned value.
* @return the angular velocity of the GRS80 ellipsoid, &omega;, in rad
* s<sup>&minus;1</sup>.
*
* This is about 2 &pi; 366.25 / (365.25 &times; 24 &times; 3600) rad
* s<sup>&minus;1</sup>. 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<typename T> static T GRS80_omega()
{ return 7292115 / (T(1000000) * 100000); }
/**
* A synonym for GRS80_omega<real>().
**********************************************************************/
static real GRS80_omega() { return GRS80_omega<real>(); }
/**
* @tparam T the type of the returned value.
* @return the dynamical form factor of the GRS80 ellipsoid,
* <i>J</i><sub>2</sub>.
**********************************************************************/
template<typename T> static T GRS80_J2()
{ return T(108263) / 100000000; }
/**
* A synonym for GRS80_J2<real>().
**********************************************************************/
static real GRS80_J2() { return GRS80_J2<real>(); }
/**
* @tparam T the type of the returned value.
* @return the central scale factor for UTM (0.9996).
**********************************************************************/
template<typename T> static T UTM_k0()
{return T(9996) / 10000; }
/**
* A synonym for UTM_k0<real>().
**********************************************************************/
static real UTM_k0() { return UTM_k0<real>(); }
/**
* @tparam T the type of the returned value.
* @return the central scale factor for UPS (0.994).
**********************************************************************/
template<typename T> static T UPS_k0()
{ return T(994) / 1000; }
/**
* A synonym for UPS_k0<real>().
**********************************************************************/
static real UPS_k0() { return UPS_k0<real>(); }
///@}
/** \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<typename T> static T meter() { return T(1); }
/**
* A synonym for meter<real>().
**********************************************************************/
static real meter() { return meter<real>(); }
/**
* @return the number of meters in a kilometer.
**********************************************************************/
static real kilometer()
{ return 1000 * meter<real>(); }
/**
* @return the number of meters in a nautical mile (approximately 1 arc
* minute)
**********************************************************************/
static real nauticalmile()
{ return 1852 * meter<real>(); }
/**
* @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<typename T> static T square_meter()
{ return meter<real>() * meter<real>(); }
/**
* A synonym for square_meter<real>().
**********************************************************************/
static real square_meter()
{ return square_meter<real>(); }
/**
* @return the number of square meters in a hectare.
**********************************************************************/
static real hectare()
{ return 10000 * square_meter<real>(); }
/**
* @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<real>(); }
/**
* @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<real>(); }
///@}
};
/**
* \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
src/Geo/MGRS.cpp 0 → 100644
View file @ 1ee854e1
/**
* \file MGRS.cpp
* \brief Implementation for GeographicLib::MGRS class
*
* Copyright (c) Charles Karney (2008-2017) <charles@karney.com> 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<long long>::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
src/Geo/MGRS.hpp 0 → 100644
View file @ 1ee854e1
/**
* \file MGRS.hpp
* \brief Header for GeographicLib::MGRS class
*
* Copyright (c) Charles Karney (2008-2019) <charles@karney.com> 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,
* <a href="http://earth-info.nga.mil/GandG/publications/tm8358.1/pdf/TM8358_1.pdf">
* Datums, Ellipsoids, Grids, and Grid Reference Systems</a>,
* Defense Mapping Agency, Technical Manual TM8358.1 (1990).
* .
* This document has been updated by the two NGA documents
* - <a href="http://earth-info.nga.mil/GandG/publications/NGA_STND_0037_2_0_0_GRIDS/NGA.STND.0037_2.0.0_GRIDS.pdf">
* Universal Grids and Grid Reference Systems</a>,
* NGA.STND.0037_2.0.0_GRIDS (2014).
* - <a href="http://earth-info.nga.mil/GandG/publications/NGA_SIG_0012_2_0_0_UTMUPS/NGA.SIG.0012_2.0.0_UTMUPS.pdf">
* The Universal Grids and the Transverse Mercator and Polar Stereographic
* Map Projections</a>, 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 &lfloor;10<sup>6</sup> <i>x</i>&rfloor; 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 <a href="http://www.nga.mil">NGA</a> software package
* <a href="http://earth-info.nga.mil/GandG/geotrans/index.html">geotrans</a>
* 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 = &minus;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
* - &hellip;
* - \e prec = 11 (max), 1 &mu;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 [&minus;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 = &minus;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 &minus;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
* &minus;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
src/Geo/Math.cpp 0 → 100644
View file @ 1ee854e1
/**
* \file Math.cpp
* \brief Implementation for GeographicLib::Math class
*
* Copyright (c) Charles Karney (2015-2019) <charles@karney.com> 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<real>::digits;
#else
return std::numeric_limits<real>::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<real>::digits10;
#else
return std::numeric_limits<real>::digits10();
#endif
}
int Math::extra_digits() {
return
digits10() > std::numeric_limits<double>::digits10 ?
digits10() - std::numeric_limits<double>::digits10 : 0;
}
template<typename T> 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<typename T> 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<typename T> 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<typename T> 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<typename T> 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<typename T> 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<typename T> 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<typename T> 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<typename T> 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<typename T> 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<typename T> 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<long>::min();
x = round(x);
if (abs(x) < -T(r)) // Assume T(LONG_MIN) is exact
r = long(x);
return r;
#endif
}
template<typename T> 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<typename T> 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<typename T> 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<typename T> 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<T>();
// 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<typename T> T Math::sind(T x) {
// See sincosd
T r; int q;
r = remquo(x, T(90), &q); // now abs(r) <= 45
r *= degree<T>();
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<typename T> T Math::cosd(T x) {
// See sincosd
T r; int q;
r = remquo(x, T(90), &q); // now abs(r) <= 45
r *= degree<T>();
unsigned p = unsigned(q + 1);
r = p & 1U ? cos(r) : sin(r);
if (p & 2U) r = -r;
return T(0) + r;
}
template<typename T> T Math::tand(T x) {
static const T overflow = 1 / sq(std::numeric_limits<T>::epsilon());
T s, c;
sincosd(x, s, c);
return c != 0 ? s / c : (s < 0 ? -overflow : overflow);
}
template<typename T> 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<T>();
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<typename T> T Math::atand(T x)
{ return atan2d(x, T(1)); }
template<typename T> T Math::eatanhe(T x, T es) {
return es > T(0) ? es * atanh(es * x) : -es * atan(es * x);
}
template<typename T> 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<typename T> T Math::tauf(T taup, T es) {
const int numit = 5;
const T tol = sqrt(numeric_limits<T>::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<typename T> 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<T>::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<T>::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<T>::max();
#endif
#endif
}
template<typename T> T Math::NaN() {
#if defined(_MSC_VER)
return std::numeric_limits<T>::has_quiet_NaN ?
std::numeric_limits<T>::quiet_NaN() :
(std::numeric_limits<T>::max)();
#else
return std::numeric_limits<T>::has_quiet_NaN ?
std::numeric_limits<T>::quiet_NaN() :
std::numeric_limits<T>::max();
#endif
}
template<typename T> bool Math::isnan(T x) {
#if GEOGRAPHICLIB_CXX11_MATH
using std::isnan; return isnan(x);
#else
return x != x;
#endif
}
template<typename T> T Math::infinity() {
#if defined(_MSC_VER)
return std::numeric_limits<T>::has_infinity ?
std::numeric_limits<T>::infinity() :
(std::numeric_limits<T>::max)();
#else
return std::numeric_limits<T>::has_infinity ?
std::numeric_limits<T>::infinity() :
std::numeric_limits<T>::max();
#endif
}
/// \cond SKIP
// Instantiate
template Math::real GEOGRAPHICLIB_EXPORT
Math::hypot<Math::real>(Math::real, Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::expm1<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::log1p<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::asinh<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::atanh<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::cbrt<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::remainder<Math::real>(Math::real, Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::remquo<Math::real>(Math::real, Math::real, int*);
template Math::real GEOGRAPHICLIB_EXPORT
Math::round<Math::real>(Math::real);
template long GEOGRAPHICLIB_EXPORT
Math::lround<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::copysign<Math::real>(Math::real, Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::fma<Math::real>(Math::real, Math::real, Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::sum<Math::real>(Math::real, Math::real, Math::real&);
template Math::real GEOGRAPHICLIB_EXPORT
Math::AngRound<Math::real>(Math::real);
template void GEOGRAPHICLIB_EXPORT
Math::sincosd<Math::real>(Math::real, Math::real&, Math::real&);
template Math::real GEOGRAPHICLIB_EXPORT
Math::sind<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::cosd<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::tand<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::atan2d<Math::real>(Math::real, Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::atand<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::eatanhe<Math::real>(Math::real, Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::taupf<Math::real>(Math::real, Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::tauf<Math::real>(Math::real, Math::real);
template bool GEOGRAPHICLIB_EXPORT
Math::isfinite<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::NaN<Math::real>();
template bool GEOGRAPHICLIB_EXPORT
Math::isnan<Math::real>(Math::real);
template Math::real GEOGRAPHICLIB_EXPORT
Math::infinity<Math::real>();
#if GEOGRAPHICLIB_PRECISION != 2
// Always have double versions available
template double GEOGRAPHICLIB_EXPORT
Math::hypot<double>(double, double);
template double GEOGRAPHICLIB_EXPORT
Math::expm1<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::log1p<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::asinh<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::atanh<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::cbrt<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::remainder<double>(double, double);
template double GEOGRAPHICLIB_EXPORT
Math::remquo<double>(double, double, int*);
template double GEOGRAPHICLIB_EXPORT
Math::round<double>(double);
template long GEOGRAPHICLIB_EXPORT
Math::lround<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::copysign<double>(double, double);
template double GEOGRAPHICLIB_EXPORT
Math::fma<double>(double, double, double);
template double GEOGRAPHICLIB_EXPORT
Math::sum<double>(double, double, double&);
template double GEOGRAPHICLIB_EXPORT
Math::AngRound<double>(double);
template void GEOGRAPHICLIB_EXPORT
Math::sincosd<double>(double, double&, double&);
template double GEOGRAPHICLIB_EXPORT
Math::sind<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::cosd<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::tand<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::atan2d<double>(double, double);
template double GEOGRAPHICLIB_EXPORT
Math::atand<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::eatanhe<double>(double, double);
template double GEOGRAPHICLIB_EXPORT
Math::taupf<double>(double, double);
template double GEOGRAPHICLIB_EXPORT
Math::tauf<double>(double, double);
template bool GEOGRAPHICLIB_EXPORT
Math::isfinite<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::NaN<double>();
template bool GEOGRAPHICLIB_EXPORT
Math::isnan<double>(double);
template double GEOGRAPHICLIB_EXPORT
Math::infinity<double>();
#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, long double);
template long double GEOGRAPHICLIB_EXPORT
Math::expm1<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::log1p<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::asinh<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::atanh<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::cbrt<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::remainder<long double>(long double, long double);
template long double GEOGRAPHICLIB_EXPORT
Math::remquo<long double>(long double, long double, int*);
template long double GEOGRAPHICLIB_EXPORT
Math::round<long double>(long double);
template long GEOGRAPHICLIB_EXPORT
Math::lround<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::copysign<long double>(long double, long double);
template long double GEOGRAPHICLIB_EXPORT
Math::fma<long double>(long double, long double, long double);
template long double GEOGRAPHICLIB_EXPORT
Math::sum<long double>(long double, long double, long double&);
template long double GEOGRAPHICLIB_EXPORT
Math::AngRound<long double>(long double);
template void GEOGRAPHICLIB_EXPORT
Math::sincosd<long double>(long double, long double&, long double&);
template long double GEOGRAPHICLIB_EXPORT
Math::sind<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::cosd<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::tand<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::atan2d<long double>(long double, long double);
template long double GEOGRAPHICLIB_EXPORT
Math::atand<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::eatanhe<long double>(long double, long double);
template long double GEOGRAPHICLIB_EXPORT
Math::taupf<long double>(long double, long double);
template long double GEOGRAPHICLIB_EXPORT
Math::tauf<long double>(long double, long double);
template bool GEOGRAPHICLIB_EXPORT
Math::isfinite<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::NaN<long double>();
template bool GEOGRAPHICLIB_EXPORT
Math::isnan<long double>(long double);
template long double GEOGRAPHICLIB_EXPORT
Math::infinity<long double>();
#endif
// Also we need int versions for Utility::nummatch
template int GEOGRAPHICLIB_EXPORT Math::NaN<int>();
template int GEOGRAPHICLIB_EXPORT Math::infinity<int>();
/// \endcond
} // namespace GeographicLib
src/Geo/Math.hpp 0 → 100644
View file @ 1ee854e1
/**
* \file Math.hpp
* \brief Header for GeographicLib::Math class
*
* Copyright (c) Charles Karney (2008-2019) <charles@karney.com> 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 <cmath>
#include <algorithm>
#include <limits>
#if GEOGRAPHICLIB_PRECISION == 4
#include <boost/version.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/math/special_functions.hpp>
#elif GEOGRAPHICLIB_PRECISION == 5
#include <mpreal.h>
#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). (<b>CAUTION</b>: 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 &pi;.
**********************************************************************/
template<typename T> static T pi() {
using std::atan2;
static const T pi = atan2(T(0), T(-1));
return pi;
}
/**
* A synonym for pi<real>().
**********************************************************************/
static real pi() { return pi<real>(); }
/**
* @tparam T the type of the returned value.
* @return the number of radians in a degree.
**********************************************************************/
template<typename T> static T degree() {
static const T degree = pi<T>() / 180;
return degree;
}
/**
* A synonym for degree<real>().
**********************************************************************/
static real degree() { return degree<real>(); }
/**
* Square a number.
*
* @tparam T the type of the argument and the returned value.
* @param[in] x
* @return <i>x</i><sup>2</sup>.
**********************************************************************/
template<typename T> 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(<i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>).
**********************************************************************/
template<typename T> static T hypot(T x, T y);
/**
* exp(\e x) &minus; 1 accurate near \e x = 0.
*
* @tparam T the type of the argument and the returned value.
* @param[in] x
* @return exp(\e x) &minus; 1.
**********************************************************************/
template<typename T> 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<typename T> 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<typename T> 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<typename T> 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 = &minus;0, returning
* &minus|<i>x</i>|.
**********************************************************************/
template<typename T> 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<typename T> 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 [&minus;\e y/2, \e y/2].
**********************************************************************/
template<typename T> 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 [&minus;\e y/2, \e y/2].
**********************************************************************/
template<typename T> 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<typename T> 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<typename T> 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 <i>xy</i> + <i>z</i>, correctly rounded (on those platforms with
* support for the <code>fma</code> instruction).
*
* On platforms without the <code>fma</code> instruction, no attempt is
* made to improve on the result of a rounded multiplication followed by a
* rounded addition.
**********************************************************************/
template<typename T> 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 <i>x</i>/hypot(<i>x</i>, <i>y</i>).
* @param[in,out] y on output set to <i>y</i>/hypot(<i>x</i>, <i>y</i>).
**********************************************************************/
template<typename T> 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<typename T> 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 <i>y</i> = &sum;<sub><i>n</i>=0..<i>N</i></sub>
* <i>p</i><sub><i>n</i></sub> <i>x</i><sup><i>N</i>&minus;<i>n</i></sup>.
* Return 0 if \e N &lt; 0. Return <i>p</i><sub>0</sub>, if \e N = 0 (even
* if \e x is infinite or a nan). The evaluation uses Horner's method.
**********************************************************************/
template<typename T> 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 (&minus;180&deg;, 180&deg;].
*
* The range of \e x is unrestricted.
**********************************************************************/
template<typename T> 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 [&minus;90&deg;, 90&deg;], otherwise
* return NaN.
**********************************************************************/
template<typename T> static T LatFix(T x)
{ using std::abs; return abs(x) > 90 ? NaN<T>() : x; }
/**
* The exact difference of two angles reduced to
* (&minus;180&deg;, 180&deg;].
*
* @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 &minus; \e x.
*
* This computes \e z = \e y &minus; \e x exactly, reduced to
* (&minus;180&deg;, 180&deg;]; 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 = &minus;180, then \e e &gt; 0; If \e d = 180, then \e e
* &le; 0.
**********************************************************************/
template<typename T> 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 [&minus;180&deg;, 180&deg;]
*
* @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 &minus; \e x, reduced to the range [&minus;180&deg;,
* 180&deg;].
*
* The result is equivalent to computing the difference exactly, reducing
* it to (&minus;180&deg;, 180&deg;] and rounding the result. Note that
* this prescription allows &minus;180&deg; to be returned (e.g., if \e x
* is tiny and negative and \e y = 180&deg;).
**********************************************************************/
template<typename T> 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 &minus; nextafter(1/16, 0) =
* 1/2<sup>57</sup> 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&deg;.) We use this to avoid having to deal with near
* singular cases when \e x is non-zero but tiny (e.g.,
* 10<sup>&minus;200</sup>). This converts &minus;0 to +0; however tiny
* negative numbers get converted to &minus;0.
**********************************************************************/
template<typename T> 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(<i>x</i>).
* @param[out] cosx cos(<i>x</i>).
*
* The results obey exactly the elementary properties of the trigonometric
* functions, e.g., sin 9&deg; = cos 81&deg; = &minus; sin 123456789&deg;.
* If x = &minus;0, then \e sinx = &minus;0; this is the only case where
* &minus;0 is returned.
**********************************************************************/
template<typename T> 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(<i>x</i>).
**********************************************************************/
template<typename T> 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(<i>x</i>).
**********************************************************************/
template<typename T> 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(<i>x</i>).
*
* If \e x = &plusmn;90&deg;, then a suitably large (but finite) value is
* returned.
**********************************************************************/
template<typename T> 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(<i>y</i>, <i>x</i>) in degrees.
*
* The result is in the range (&minus;180&deg; 180&deg;]. N.B.,
* atan2d(&plusmn;0, &minus;1) = +180&deg;; atan2d(&minus;&epsilon;,
* &minus;1) = &minus;180&deg;, for &epsilon; positive and tiny;
* atan2d(&plusmn;0, +1) = &plusmn;0&deg;.
**********************************************************************/
template<typename T> 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(<i>x</i>) in degrees.
**********************************************************************/
template<typename T> static T atand(T x);
/**
* Evaluate <i>e</i> atanh(<i>e x</i>)
*
* @tparam T the type of the argument and the returned value.
* @param[in] x
* @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
* sqrt(|<i>e</i><sup>2</sup>|)
* @return <i>e</i> atanh(<i>e x</i>)
*
* If <i>e</i><sup>2</sup> is negative (<i>e</i> is imaginary), the
* expression is evaluated in terms of atan.
**********************************************************************/
template<typename T> static T eatanhe(T x, T es);
/**
* tan&chi; in terms of tan&phi;
*
* @tparam T the type of the argument and the returned value.
* @param[in] tau &tau; = tan&phi;
* @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
* sqrt(|<i>e</i><sup>2</sup>|)
* @return &tau;&prime; = tan&chi;
*
* See Eqs. (7--9) of
* C. F. F. Karney,
* <a href="https://doi.org/10.1007/s00190-011-0445-3">
* Transverse Mercator with an accuracy of a few nanometers,</a>
* J. Geodesy 85(8), 475--485 (Aug. 2011)
* (preprint
* <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
**********************************************************************/
template<typename T> static T taupf(T tau, T es);
/**
* tan&phi; in terms of tan&chi;
*
* @tparam T the type of the argument and the returned value.
* @param[in] taup &tau;&prime; = tan&chi;
* @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
* sqrt(|<i>e</i><sup>2</sup>|)
* @return &tau; = tan&phi;
*
* See Eqs. (19--21) of
* C. F. F. Karney,
* <a href="https://doi.org/10.1007/s00190-011-0445-3">
* Transverse Mercator with an accuracy of a few nanometers,</a>
* J. Geodesy 85(8), 475--485 (Aug. 2011)
* (preprint
* <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
**********************************************************************/
template<typename T> 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<typename T> 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<typename T> static T NaN();
/**
* A synonym for NaN<real>().
**********************************************************************/
static real NaN() { return NaN<real>(); }
/**
* Test for NaN.
*
* @tparam T the type of the argument.
* @param[in] x
* @return true if argument is a NaN.
**********************************************************************/
template<typename T> static bool isnan(T x);
/**
* Infinity
*
* @tparam T the type of the returned value.
* @return infinity if available, otherwise return the max real.
**********************************************************************/
template<typename T> static T infinity();
/**
* A synonym for infinity<real>().
**********************************************************************/
static real infinity() { return infinity<real>(); }
/**
* 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<typename T> 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
src/Geo/PolarStereographic.cpp 0 → 100644
View file @ 1ee854e1
/**
* \file PolarStereographic.cpp
* \brief Implementation for GeographicLib::PolarStereographic class
*
* Copyright (c) Charles Karney (2008-2017) <charles@karney.com> 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<real>::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
src/Geo/PolarStereographic.hpp 0 → 100644
View file @ 1ee854e1
/**
* \file PolarStereographic.hpp
* \brief Header for GeographicLib::PolarStereographic class
*
* Copyright (c) Charles Karney (2008-2019) <charles@karney.com> 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,
* <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projections: A
* Working Manual</a>, 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 &minus; \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 (&minus;90&deg;,
* 90&deg;].
**********************************************************************/
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
* (&minus;90&deg;, 90&deg;] for \e northp = true and in the range
* [&minus;90&deg;, 90&deg;) 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 [&minus;180&deg;, 180&deg;].
**********************************************************************/
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
src/Geo/TransverseMercator.cpp 0 → 100644
View file @ 1ee854e1
/**
* \file TransverseMercator.cpp
* \brief Implementation for GeographicLib::TransverseMercator class
*
* Copyright (c) Charles Karney (2008-2017) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
*
* This implementation follows closely JHS 154, ETRS89 -
* j&auml;rjestelm&auml;&auml;n liittyv&auml;t karttaprojektiot,
* tasokoordinaatistot ja karttalehtijako</a> (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 <iostream>
#include <complex>
#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<real> a(2 * c0 * ch0, -2 * s0 * sh0); // 2 * cos(2*zeta')
int n = maxpow_;
complex<real>
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<real>(s0 * ch0, c0 * sh0); // sin(2*zeta')
y1 = complex<real>(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<real> a(2 * c0 * ch0, -2 * s0 * sh0); // 2 * cos(2*zeta)
int n = maxpow_;
complex<real>
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<real>(s0 * ch0, c0 * sh0); // sin(2*zeta)
y1 = complex<real>(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
src/Geo/TransverseMercator.hpp 0 → 100644
View file @ 1ee854e1
/**
* \file TransverseMercator.hpp
* \brief Header for GeographicLib::TransverseMercator class
*
* Copyright (c) Charles Karney (2008-2019) <charles@karney.com> 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&uuml;ger's method which evaluates the projection and its
* inverse in terms of a series. See
* - L. Kr&uuml;ger,
* <a href="https://doi.org/10.2312/GFZ.b103-krueger28"> Konforme
* Abbildung des Erdellipsoids in der Ebene</a> (Conformal mapping of the
* ellipsoidal earth to the plane), Royal Prussian Geodetic Institute, New
* Series 52, 172 pp. (1912).
* - C. F. F. Karney,
* <a href="https://doi.org/10.1007/s00190-011-0445-3">
* Transverse Mercator with an accuracy of a few nanometers,</a>
* J. Geodesy 85(8), 475--485 (Aug. 2011);
* preprint
* <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>.
*
* Kr&uuml;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
* &times; 10<sup>&minus;15</sup>&quot; and the relative error in the scale
* is 6 &times; 10<sup>&minus;12</sup>%%. See Sec. 4 of
* <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a> for details.
* The speed penalty in going to 6th order is only about 1%.
*
* There's a singularity in the projection at &phi; = 0&deg;, &lambda;
* &minus; &lambda;<sub>0</sub> = &plusmn;(1 &minus; \e e)90&deg; (&asymp;
* &plusmn;82.6&deg; 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 &minus; 2e)90&deg; (&asymp; 75&deg; 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 (<a href="https://www.spatialreference.org/ref/epsg/7405/">
* EPSG:7405</a>) 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
*
* <a href="TransverseMercatorProj.1.html">TransverseMercatorProj</a> 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 &minus; \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
* [&minus;90&deg;, 90&deg;].
**********************************************************************/
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 [&minus;180&deg;, 180&deg;].
**********************************************************************/
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
src/Geo/UTMUPS.cpp 0 → 100644
View file @ 1ee854e1
/**
* \file UTMUPS.cpp
* \brief Implementation for GeographicLib::UTMUPS class
*
* Copyright (c) Charles Karney (2008-2019) <charles@karney.com> 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
src/Geo/UTMUPS.hpp 0 → 100644
View file @ 1ee854e1
/**
* \file UTMUPS.hpp
* \brief Header for GeographicLib::UTMUPS class
*
* Copyright (c) Charles Karney (2008-2019) <charles@karney.com> 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,
* <a href="http://earth-info.nga.mil/GandG/publications/tm8358.2/TM8358_2.pdf">
* The Universal Grids: Universal Transverse Mercator (UTM) and Universal
* Polar Stereographic (UPS)</a>, 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
* - <a href="http://earth-info.nga.mil/GandG/publications/NGA_SIG_0012_2_0_0_UTMUPS/NGA.SIG.0012_2.0.0_UTMUPS.pdf">
* The Universal Grids and the Transverse Mercator and Polar Stereographic
* Map Projections</a>, 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 <a href="http://www.nga.mil">NGA</a> software package
* <a href="http://earth-info.nga.mil/GandG/geotrans/index.html">geotrans</a>
* 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 &minus;80&deg; and 84&deg;
* 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&deg;, 48&deg;). 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 [&minus;80&deg;, 84&deg;), 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 [&minus;80&deg;,
* 84&deg;) and longitude is in [42&deg;, 48&deg;).
**********************************************************************/
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] = [&minus;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 [&minus;90&deg;,
* 90&deg;].
* @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 &ge; 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 &lt; 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&deg; &le; \e lat &lt; 40&deg;. 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 (10<sup>7</sup>).
**********************************************************************/
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
src/Geo/Utility.cpp 0 → 100644
View file @ 1ee854e1
/**
* \file Utility.cpp
* \brief Implementation for GeographicLib::Utility class
*
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#include <cstdlib>
#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
src/Geo/Utility.h 0 → 100644
View file @ 1ee854e1
/**
* \file Utility.hpp
* \brief Header for GeographicLib::Utility class
*
* Copyright (c) Charles Karney (2011-2019) <charles@karney.com> 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 <GeographicLib/Constants.hpp>
#include <iomanip>
#include <vector>
#include <sstream>
#include <cctype>
#include <ctime>
#include <cstring>
#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<int>(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<int>(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<int>(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<int>(s.substr(p1, p2 - p1));
if (++p2 == s.size())
throw GeographicErr("Empty day field in date " + s);
d1 = val<int>(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<typename T> static T fractionalyear(const std::string& s) {
try {
return val<T>(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 &minus;1).
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return the string representation.
*
* If \e p &ge; 0, then the number fixed format is used with p bits of
* precision. With p < 0, there is no manipulation of the format.
**********************************************************************/
template<typename T> 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 &minus;1).
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return the string representation.
*
* If \e p &ge; 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<T> 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<typename T> static T val(const std::string& s) {
// If T is bool, then the specialization val<bool>() 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<T>::is_integer ? 0 : nummatch<T>(t);
if (x == 0)
throw GeographicErr(errmsg);
return x;
}
/**
* \deprecated An old name for val<T>(s).
**********************************************************************/
template<typename T>
GEOGRAPHICLIB_DEPRECATED("Use Utility::val<T>(s)")
static T num(const std::string& s) {
return val<T>(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 (&plusmn;&infin;, nan) or 0 if none is
* found.
*
* White space is not allowed at the beginning or end of \e s.
**********************************************************************/
template<typename T> 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<T>();
else if (t == "INF" || t == "1.#INF")
return sign * Math::infinity<T>();
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<typename T> 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<T>(s) :
// delim in [1, size() - 2]
val<T>(s.substr(0, delim)) / val<T>(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
* &minus;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
* &minus;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<typename ExtT, typename IntT, bool bigendp>
static void readarray(std::istream& str, IntT array[], size_t num) {
#if GEOGRAPHICLIB_PRECISION < 4
if (sizeof(IntT) == sizeof(ExtT) &&
std::numeric_limits<IntT>::is_integer ==
std::numeric_limits<ExtT>::is_integer)
{
// Data is compatible (aside from the issue of endian-ness).
str.read(reinterpret_cast<char*>(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<IntT>(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<char*>(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<ExtT>(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<typename ExtT, typename IntT, bool bigendp>
static void readarray(std::istream& str, std::vector<IntT>& array) {
if (array.size() > 0)
readarray<ExtT, IntT, bigendp>(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<typename ExtT, typename IntT, bool bigendp>
static void writearray(std::ostream& str, const IntT array[], size_t num)
{
#if GEOGRAPHICLIB_PRECISION < 4
if (sizeof(IntT) == sizeof(ExtT) &&
std::numeric_limits<IntT>::is_integer ==
std::numeric_limits<ExtT>::is_integer &&
bigendp == Math::bigendian)
{
// Data is compatible (including endian-ness).
str.write(reinterpret_cast<const char*>(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>(ExtT(array[i++]));
str.write(reinterpret_cast<const char*>(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<typename ExtT, typename IntT, bool bigendp>
static void writearray(std::ostream& str, std::vector<IntT>& array) {
if (array.size() > 0)
writearray<ExtT, IntT, bigendp>(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<T>() for strings.
**********************************************************************/
template<> inline std::string Utility::val<std::string>(const std::string& s)
{ return trim(s); }
/**
* The specialization of Utility::val<T>() for bools.
**********************************************************************/
template<> inline bool Utility::val<bool>(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
src/Geo/Utility.hpp 0 → 100644
View file @ 1ee854e1
/**
* \file Utility.hpp
* \brief Header for GeographicLib::Utility class
*
* Copyright (c) Charles Karney (2011-2019) <charles@karney.com> 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 <iomanip>
#include <vector>
#include <sstream>
#include <cctype>
#include <ctime>
#include <cstring>
#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<int>(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<int>(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<int>(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<int>(s.substr(p1, p2 - p1));
if (++p2 == s.size())
throw GeographicErr("Empty day field in date " + s);
d1 = val<int>(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<typename T> static T fractionalyear(const std::string& s) {
try {
return val<T>(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 &minus;1).
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return the string representation.
*
* If \e p &ge; 0, then the number fixed format is used with p bits of
* precision. With p < 0, there is no manipulation of the format.
**********************************************************************/
template<typename T> 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 &minus;1).
* @exception std::bad_alloc if memory for the string can't be allocated.
* @return the string representation.
*
* If \e p &ge; 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<T> 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<typename T> static T val(const std::string& s) {
// If T is bool, then the specialization val<bool>() 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<T>::is_integer ? 0 : nummatch<T>(t);
if (x == 0)
throw GeographicErr(errmsg);
return x;
}
/**
* \deprecated An old name for val<T>(s).
**********************************************************************/
template<typename T>
GEOGRAPHICLIB_DEPRECATED("Use Utility::val<T>(s)")
static T num(const std::string& s) {
return val<T>(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 (&plusmn;&infin;, nan) or 0 if none is
* found.
*
* White space is not allowed at the beginning or end of \e s.
**********************************************************************/
template<typename T> 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<T>();
else if (t == "INF" || t == "1.#INF")
return sign * Math::infinity<T>();
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<typename T> 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<T>(s) :
// delim in [1, size() - 2]
val<T>(s.substr(0, delim)) / val<T>(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
* &minus;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
* &minus;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<typename ExtT, typename IntT, bool bigendp>
static void readarray(std::istream& str, IntT array[], size_t num) {
#if GEOGRAPHICLIB_PRECISION < 4
if (sizeof(IntT) == sizeof(ExtT) &&
std::numeric_limits<IntT>::is_integer ==
std::numeric_limits<ExtT>::is_integer)
{
// Data is compatible (aside from the issue of endian-ness).
str.read(reinterpret_cast<char*>(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<IntT>(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<char*>(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<ExtT>(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<typename ExtT, typename IntT, bool bigendp>
static void readarray(std::istream& str, std::vector<IntT>& array) {
if (array.size() > 0)
readarray<ExtT, IntT, bigendp>(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<typename ExtT, typename IntT, bool bigendp>
static void writearray(std::ostream& str, const IntT array[], size_t num)
{
#if GEOGRAPHICLIB_PRECISION < 4
if (sizeof(IntT) == sizeof(ExtT) &&
std::numeric_limits<IntT>::is_integer ==
std::numeric_limits<ExtT>::is_integer &&
bigendp == Math::bigendian)
{
// Data is compatible (including endian-ness).
str.write(reinterpret_cast<const char*>(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>(ExtT(array[i++]));
str.write(reinterpret_cast<const char*>(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<typename ExtT, typename IntT, bool bigendp>
static void writearray(std::ostream& str, std::vector<IntT>& array) {
if (array.size() > 0)
writearray<ExtT, IntT, bigendp>(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<T>() for strings.
**********************************************************************/
template<> inline std::string Utility::val<std::string>(const std::string& s)
{ return trim(s); }
/**
* The specialization of Utility::val<T>() for bools.
**********************************************************************/
template<> inline bool Utility::val<bool>(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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment