Merge pull request #8008 from mavlink/pr-mgrs

MGRS coordinates

qgroundcontrol.pro

@@ -407,6 +407,7 @@ INCLUDEPATH += \
     src/FlightMap \
     src/FlightMap/Widgets \
     src/FollowMe \
+    src/Geo \
     src/GPS \
     src/Joystick \
     src/PlanView \
@@ -631,12 +632,19 @@ HEADERS += \
     src/MissionManager/VisualMissionItem.h \
     src/PositionManager/PositionManager.h \
     src/PositionManager/SimulatedPosition.h \
+    src/Geo/QGCGeo.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 \
@@ -693,7 +701,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 \
@@ -864,11 +871,17 @@ SOURCES += \
     src/MissionManager/VisualMissionItem.cc \
     src/PositionManager/PositionManager.cpp \
     src/PositionManager/SimulatedPosition.cc \
+    src/Geo/QGCGeo.cc \
+    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 \
@@ -924,7 +937,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/PolarStereographic.cpp

+/**
+ * \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

+/**
+ * \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/QGCGeo.cc → src/Geo/QGCGeo.cc

@@ -8,25 +8,21 @@
  ****************************************************************************/
 
 #include <QDebug>
+#include <QString>
 
 #include <cmath>
 #include <limits>
 
 #include "QGCGeo.h"
-#include "UTM.h"
+#include "UTMUPS.hpp"
+#include "MGRS.hpp"
 
 // These defines are private
 #define M_DEG_TO_RAD (M_PI / 180.0)
-
 #define M_RAD_TO_DEG (180.0 / M_PI)
+#define CONSTANTS_RADIUS_OF_EARTH 6371000 // meters (m)
 
-#define CONSTANTS_ONE_G					9.80665f		/* m/s^2		*/
-#define CONSTANTS_AIR_DENSITY_SEA_LEVEL_15C		1.225f			/* kg/m^3		*/
-#define CONSTANTS_AIR_GAS_CONST				287.1f 			/* J/(kg * K)		*/
-#define CONSTANTS_ABSOLUTE_NULL_CELSIUS			-273.15f		/* °C			*/
-#define CONSTANTS_RADIUS_OF_EARTH			6371000			/* meters (m)		*/
-
-static const float epsilon = std::numeric_limits<double>::epsilon();
+static const double epsilon = std::numeric_limits<double>::epsilon();
 
 void convertGeoToNed(QGeoCoordinate coord, QGeoCoordinate origin, double* x, double* y, double* z)
 {
@@ -61,7 +57,7 @@ void convertGeoToNed(QGeoCoordinate coord, QGeoCoordinate origin, double* x, dou
 void convertNedToGeo(double x, double y, double z, QGeoCoordinate origin, QGeoCoordinate *coord) {
     double x_rad = x / CONSTANTS_RADIUS_OF_EARTH;
     double y_rad = y / CONSTANTS_RADIUS_OF_EARTH;
-    double c = sqrtf(x_rad * x_rad + y_rad * y_rad);
+    double c = sqrt(x_rad * x_rad + y_rad * y_rad);
     double sin_c = sin(c);
     double cos_c = cos(c);
 
@@ -91,14 +87,70 @@ void convertNedToGeo(double x, double y, double z, QGeoCoordinate origin, QGeoCo
 
 int convertGeoToUTM(const QGeoCoordinate& coord, double& easting, double& northing)
 {
-    return LatLonToUTMXY(coord.latitude(), coord.longitude(), -1 /* zone */, easting, northing);
+    try {
+        int zone;
+        bool northp;
+        GeographicLib::UTMUPS::Forward(coord.latitude(), coord.longitude(), zone, northp, easting, northing);
+        return zone;
+    } catch(...) {
+        return 0;
+    }
+}
+
+bool convertUTMToGeo(double easting, double northing, int zone, bool southhemi, QGeoCoordinate& coord)
+{
+    double lat, lon;
+    try {
+        GeographicLib::UTMUPS::Reverse(zone, !southhemi, easting, northing, lat, lon);
+    } catch(...) {
+        return false;
+    }
+    coord.setLatitude(lat);
+    coord.setLongitude(lon);
+
+    return true;
+}
+
+QString convertGeoToMGRS(const QGeoCoordinate& coord)
+{
+    int zone;
+    bool northp;
+    double x, y;
+    std::string mgrs;
+
+    try {
+        GeographicLib::UTMUPS::Forward(coord.latitude(), coord.longitude(), zone, northp, x, y);
+        GeographicLib::MGRS::Forward(zone, northp, x, y, coord.latitude(), 5, mgrs);
+    } catch(...) {
+        mgrs = "";
+    }
+
+    QString qstr = QString::fromStdString(mgrs);
+    for (int i = qstr.length() - 1; i >= 0; i--) {
+        if (!qstr.at(i).isDigit()) {
+            int l = (qstr.length() - i) / 2;
+            return qstr.left(i + 1) + " " + qstr.mid(i + 1, l) + " " + qstr.mid(i + 1 + l);
+        }
+    }
+
+    return qstr;
 }
 
-void convertUTMToGeo(double easting, double northing, int zone, bool southhemi, QGeoCoordinate& coord)
+bool convertMGRSToGeo(QString mgrs, QGeoCoordinate& coord)
 {
-    double latRadians, lonRadians;
+    int zone, prec;
+    bool northp;
+    double x, y;
+    double lat, lon;
+
+    try {
+        GeographicLib::MGRS::Reverse(mgrs.simplified().replace(" ", "").toStdString(), zone, northp, x, y, prec);
+        GeographicLib::UTMUPS::Reverse(zone, northp, x, y, lat, lon);
+    } catch(...) {
+        return false;
+    }
+    coord.setLatitude(lat);
+    coord.setLongitude(lon);
 
-    UTMXYToLatLon (easting, northing, zone, southhemi, latRadians, lonRadians);
-    coord.setLatitude(RadToDeg(latRadians));
-    coord.setLongitude(RadToDeg(lonRadians));
+    return true;
 }

src/QGCGeo.h → src/Geo/QGCGeo.h

@@ -58,6 +58,7 @@ void convertNedToGeo(double x, double y, double z, QGeoCoordinate origin, QGeoCo
 //
 // Returns:
 //   The UTM zone used for calculating the values of x and y.
+//   If conversion failed the function returns 0
 int convertGeoToUTM(const QGeoCoordinate& coord, double& easting, double& northing);
 
 // UTMXYToLatLon
@@ -78,7 +79,30 @@ int convertGeoToUTM(const QGeoCoordinate& coord, double& easting, double& northi
 // lon - The longitude of the point, in radians.
 //
 // Returns:
-// The function does not return a value.
-void convertUTMToGeo(double easting, double northing, int zone, bool southhemi, QGeoCoordinate& coord);
+// The function returns true if conversion succeeded.
+bool convertUTMToGeo(double easting, double northing, int zone, bool southhemi, QGeoCoordinate& coord);
+
+// Converts a latitude/longitude pair to MGRS string
+//
+// Inputs:
+//   coord - Latitude, Longiture coordinate
+//
+// Returns:
+//   The MGRS coordinate string
+//   If conversion fails the function returns empty string
+QString convertGeoToMGRS(const QGeoCoordinate& coord);
+
+// Converts MGRS string to a latitude/longitude pair.
+//
+// Inputs:
+// mgrs - MGRS coordinate string
+//
+// Outputs:
+// lat - The latitude of the point, in radians.
+// lon - The longitude of the point, in radians.
+//
+// Returns:
+// The function returns true if conversion succeeded.
+bool convertMGRSToGeo(QString mgrs, QGeoCoordinate& coord);
 
 #endif // QGCGEO_H

src/Geo/TransverseMercator.hpp

+/**
+ * \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/Utility.cpp

+/**
+ * \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/QmlControls/EditPositionDialog.FactMetaData.json

@@ -40,5 +40,10 @@
     "type":             "uint8",
     "enumStrings":      "North,South",
     "enumValues":       "0,1"
+},
+{
+    "name":             "MGRS",
+    "shortDescription": "MGRS coordinate",
+    "type":             "string"
 }
 ]

src/QmlControls/EditPositionDialog.qml

@@ -47,80 +47,108 @@ QGCViewDialog {
                 rowSpacing:     _margin
                 columns:        2
 
-                QGCLabel { text: qsTr("Latitude") }
+                QGCLabel {
+                    text: qsTr("Latitude")
+                }
                 FactTextField {
                     fact:               controller.latitude
                     Layout.fillWidth:   true
                 }
 
-                QGCLabel { text: qsTr("Longitude") }
+                QGCLabel {
+                    text: qsTr("Longitude")
+                }
                 FactTextField {
                     fact:               controller.longitude
                     Layout.fillWidth:   true
                 }
-            }
 
-            QGCButton {
-                anchors.right:  parent.right
-                text:           qsTr("Set Geographic")
-
-                onClicked: {
-                    controller.setFromGeo()
-                    reject()
+                QGCButton {
+                    text:               qsTr("Set Geographic")
+                    Layout.alignment:   Qt.AlignRight
+                    Layout.columnSpan:  2
+                    onClicked: {
+                        controller.setFromGeo()
+                        reject()
+                    }
                 }
-            }
 
-            GridLayout {
-                anchors.left:   parent.left
-                anchors.right:  parent.right
-                columnSpacing:  _margin
-                rowSpacing:     _margin
-                columns:        2
+                Item { width: 1; height: ScreenTools.defaultFontPixelHeight; Layout.columnSpan: 2}
 
-                QGCLabel { text: qsTr("Zone") }
+                QGCLabel {
+                    text: qsTr("Zone")
+                }
                 FactTextField {
                     fact:               controller.zone
                     Layout.fillWidth:   true
                 }
 
-                QGCLabel { text: qsTr("Hemisphere") }
+                QGCLabel {
+                    text: qsTr("Hemisphere")
+                }
                 FactComboBox {
                     fact:               controller.hemisphere
                     indexModel:         false
                     Layout.fillWidth:   true
                 }
 
-                QGCLabel { text: qsTr("Easting") }
+                QGCLabel {
+                    text: qsTr("Easting")
+                }
                 FactTextField {
                     fact:               controller.easting
                     Layout.fillWidth:   true
                 }
 
-                QGCLabel { text: qsTr("Northing") }
+                QGCLabel {
+                    text: qsTr("Northing")
+                }
                 FactTextField {
                     fact:               controller.northing
                     Layout.fillWidth:   true
                 }
-            }
 
-            QGCButton {
-                anchors.right:  parent.right
-                text:           qsTr("Set UTM")
+                QGCButton {
+                    text:               qsTr("Set UTM")
+                    Layout.alignment:   Qt.AlignRight
+                    Layout.columnSpan:  2
+                    onClicked: {
+                        controller.setFromUTM()
+                        reject()
+                    }
+                }
 
-                onClicked: {
-                    controller.setFromUTM()
-                    reject()
+                Item { width: 1; height: ScreenTools.defaultFontPixelHeight; Layout.columnSpan: 2}
+
+                QGCLabel {
+                    text:              qsTr("MGRS")
+                }
+                FactTextField {
+                    fact:              controller.mgrs
+                    Layout.fillWidth:  true
                 }
-            }
 
-            QGCButton {
-                anchors.right:  parent.right
-                text:           qsTr("Set From Vehicle Position")
-                visible:        QGroundControl.multiVehicleManager.activeVehicle && QGroundControl.multiVehicleManager.activeVehicle.coordinate.isValid
+                QGCButton {
+                    text:              qsTr("Set MGRS")
+                    Layout.alignment:  Qt.AlignRight
+                    Layout.columnSpan: 2
+                    onClicked: {
+                        controller.setFromMGRS()
+                        reject()
+                    }
+                }
 
-                onClicked: {
-                    controller.setFromVehicle()
-                    reject()
+                Item { width: 1; height: ScreenTools.defaultFontPixelHeight; Layout.columnSpan: 2}
+
+                QGCButton {
+                    text:              qsTr("Set From Vehicle Position")
+                    visible:           QGroundControl.multiVehicleManager.activeVehicle && QGroundControl.multiVehicleManager.activeVehicle.coordinate.isValid
+                    Layout.alignment:  Qt.AlignRight
+                    Layout.columnSpan: 2
+                    onClicked: {
+                        controller.setFromVehicle()
+                        reject()
+                    }
                 }
             }
         } // Column

src/QmlControls/EditPositionDialogController.cc

@@ -17,6 +17,7 @@ const char*  EditPositionDialogController::_zoneFactName =          "Zone";
 const char*  EditPositionDialogController::_hemisphereFactName =    "Hemisphere";
 const char*  EditPositionDialogController::_eastingFactName =       "Easting";
 const char*  EditPositionDialogController::_northingFactName =      "Northing";
+const char*  EditPositionDialogController::_mgrsFactName =          "MGRS";
 
 QMap<QString, FactMetaData*> EditPositionDialogController::_metaDataMap;
 
@@ -27,6 +28,7 @@ EditPositionDialogController::EditPositionDialogController(void)
     , _hemisphereFact   (0, _hemisphereFactName,    FactMetaData::valueTypeUint8)
     , _eastingFact      (0, _eastingFactName,       FactMetaData::valueTypeDouble)
     , _northingFact     (0, _northingFactName,      FactMetaData::valueTypeDouble)
+    , _mgrsFact         (0, _mgrsFactName,          FactMetaData::valueTypeString)
 {
     if (_metaDataMap.isEmpty()) {
         _metaDataMap = FactMetaData::createMapFromJsonFile(QStringLiteral(":/json/EditPositionDialog.FactMetaData.json"), nullptr /* QObject parent */);
@@ -38,6 +40,7 @@ EditPositionDialogController::EditPositionDialogController(void)
     _hemisphereFact.setMetaData (_metaDataMap[_hemisphereFactName]);
     _eastingFact.setMetaData    (_metaDataMap[_eastingFactName]);
     _northingFact.setMetaData   (_metaDataMap[_northingFactName]);
+    _mgrsFact.setMetaData       (_metaDataMap[_mgrsFactName]);
 }
 
 void EditPositionDialogController::setCoordinate(QGeoCoordinate coordinate)
@@ -55,10 +58,16 @@ void EditPositionDialogController::initValues(void)
 
     double easting, northing;
     int zone = convertGeoToUTM(_coordinate, easting, northing);
-    _zoneFact.setRawValue(zone);
-    _hemisphereFact.setRawValue(_coordinate.latitude() < 0);
-    _eastingFact.setRawValue(easting);
-    _northingFact.setRawValue(northing);
+    if (zone >= 1 && zone <= 60) {
+        _zoneFact.setRawValue(zone);
+        _hemisphereFact.setRawValue(_coordinate.latitude() < 0);
+        _eastingFact.setRawValue(easting);
+        _northingFact.setRawValue(northing);
+    }
+    QString mgrs = convertGeoToMGRS(_coordinate);
+    if (!mgrs.isEmpty()) {
+        _mgrsFact.setRawValue(mgrs);
+    }
 }
 
 void EditPositionDialogController::setFromGeo(void)
@@ -71,9 +80,21 @@ void EditPositionDialogController::setFromGeo(void)
 void EditPositionDialogController::setFromUTM(void)
 {
     qDebug() << _eastingFact.rawValue().toDouble() << _northingFact.rawValue().toDouble() << _zoneFact.rawValue().toInt() << (_hemisphereFact.rawValue().toInt() == 1);
-    convertUTMToGeo(_eastingFact.rawValue().toDouble(), _northingFact.rawValue().toDouble(), _zoneFact.rawValue().toInt(), _hemisphereFact.rawValue().toInt() == 1, _coordinate);
-    qDebug() << _eastingFact.rawValue().toDouble() << _northingFact.rawValue().toDouble() << _zoneFact.rawValue().toInt() << (_hemisphereFact.rawValue().toInt() == 1) << _coordinate;
-    emit coordinateChanged(_coordinate);
+    if (convertUTMToGeo(_eastingFact.rawValue().toDouble(), _northingFact.rawValue().toDouble(), _zoneFact.rawValue().toInt(), _hemisphereFact.rawValue().toInt() == 1, _coordinate)) {
+        qDebug() << _eastingFact.rawValue().toDouble() << _northingFact.rawValue().toDouble() << _zoneFact.rawValue().toInt() << (_hemisphereFact.rawValue().toInt() == 1) << _coordinate;
+        emit coordinateChanged(_coordinate);
+    } else {
+        initValues();
+    }
+}
+
+void EditPositionDialogController::setFromMGRS(void)
+{
+    if (convertMGRSToGeo(_mgrsFact.rawValue().toString(), _coordinate)) {
+        emit coordinateChanged(_coordinate);
+    } else {
+        initValues();
+    }
 }
 
 void EditPositionDialogController::setFromVehicle(void)

src/QmlControls/EditPositionDialogController.h

@@ -28,6 +28,7 @@ public:
     Q_PROPERTY(Fact*            hemisphere  READ hemisphere                     CONSTANT)
     Q_PROPERTY(Fact*            easting     READ easting                        CONSTANT)
     Q_PROPERTY(Fact*            northing    READ northing                       CONSTANT)
+    Q_PROPERTY(Fact*            mgrs        READ mgrs                           CONSTANT)
 
     QGeoCoordinate  coordinate(void) const { return _coordinate; }
     Fact* latitude  (void) { return &_latitudeFact; }
@@ -36,12 +37,14 @@ public:
     Fact* hemisphere(void) { return &_hemisphereFact; }
     Fact* easting   (void) { return &_eastingFact; }
     Fact* northing  (void) { return &_northingFact; }
+    Fact* mgrs      (void) { return &_mgrsFact; }
 
     void setCoordinate(QGeoCoordinate coordinate);
 
     Q_INVOKABLE void initValues(void);
     Q_INVOKABLE void setFromGeo(void);
     Q_INVOKABLE void setFromUTM(void);
+    Q_INVOKABLE void setFromMGRS(void);
     Q_INVOKABLE void setFromVehicle(void);
 
 signals:
@@ -58,6 +61,7 @@ private:
     Fact _hemisphereFact;
     Fact _eastingFact;
     Fact _northingFact;
+    Fact _mgrsFact;
 
     static const char*  _latitudeFactName;
     static const char*  _longitudeFactName;
@@ -65,4 +69,5 @@ private:
     static const char*  _hemisphereFactName;
     static const char*  _eastingFactName;
     static const char*  _northingFactName;
+    static const char*  _mgrsFactName;
 };

src/SHPFileHelper.cc

@@ -8,7 +8,7 @@
  ****************************************************************************/
 
 #include "SHPFileHelper.h"
-#include "UTM.h"
+#include "QGCGeo.h"
 
 #include <QFile>
 #include <QVariant>
@@ -141,14 +141,12 @@ bool SHPFileHelper::loadPolygonFromFile(const QString& shpFile, QList<QGeoCoordi
     }
 
     for (int i=0; i<shpObject->nVertices; i++) {
-        double lat, lon;
-        if (utmZone) {
-            UTMXYToLatLon(shpObject->padfX[i], shpObject->padfY[i], utmZone, utmSouthernHemisphere, lat, lon);
-        } else {
-            lat = shpObject->padfY[i];
-            lon = shpObject->padfX[i];
+        QGeoCoordinate coord;
+        if (!utmZone || !convertUTMToGeo(shpObject->padfX[i], shpObject->padfY[i], utmZone, utmSouthernHemisphere, coord)) {
+            coord.setLatitude(shpObject->padfY[i]);
+            coord.setLongitude(shpObject->padfX[i]);
         }
-        vertices.append(QGeoCoordinate(lat, lon));
+        vertices.append(coord);
     }
 
     // Filter last vertex such that it differs from first

src/UTM.h

-// UTM.h
-
-// Original Javascript by Chuck Taylor
-// Port to C++ by Alex Hajnal
-//
-// This is a simple port of the code on the Geographic/UTM Coordinate Converter (1) page from Javascript to C++.
-// Using this you can easily convert between UTM and WGS84 (latitude and longitude).
-// Accuracy seems to be around 50cm (I suspect rounding errors are limiting precision).
-// This code is provided as-is and has been minimally tested; enjoy but use at your own risk!
-// The license for UTM.cpp and UTM.h is the same as the original Javascript: 
-// "The C++ source code in UTM.cpp and UTM.h may be copied and reused without restriction."
-// 
-// 1) http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html
-
-// QGC Note: This file has been slightly modified to prevent possible conflicts with other parts of the system
-
-#ifndef UTM_H
-#define UTM_H
-
-// DegToRad
-// Converts degrees to radians.
-double DegToRad(double deg);
-
-// RadToDeg
-// Converts radians to degrees.
-double RadToDeg(double rad);
-
-// ArcLengthOfMeridian
-// Computes the ellipsoidal distance from the equator to a point at a
-// given latitude.
-// 
-// Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
-// GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
-// 
-// Inputs:
-//     phi - Latitude of the point, in radians.
-// 
-// Globals:
-//     sm_a - Ellipsoid model major axis.
-//     sm_b - Ellipsoid model minor axis.
-// 
-// Returns:
-//     The ellipsoidal distance of the point from the equator, in meters.
-double ArcLengthOfMeridian (double phi);
-
-// UTMCentralMeridian
-// Determines the central meridian for the given UTM zone.
-//
-// Inputs:
-//     zone - An integer value designating the UTM zone, range [1,60].
-//
-// Returns:
-//   The central meridian for the given UTM zone, in radians
-//   Range of the central meridian is the radian equivalent of [-177,+177].
-double UTMCentralMeridian(int zone);
-
-// FootpointLatitude
-//
-// Computes the footpoint latitude for use in converting transverse
-// Mercator coordinates to ellipsoidal coordinates.
-//
-// Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
-//   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
-//
-// Inputs:
-//   y - The UTM northing coordinate, in meters.
-//
-// Returns:
-//   The footpoint latitude, in radians.
-double FootpointLatitude(double y);
-
-// MapLatLonToXY
-// Converts a latitude/longitude pair to x and y coordinates in the
-// Transverse Mercator projection.  Note that Transverse Mercator is not
-// the same as UTM; a scale factor is required to convert between them.
-//
-// Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
-// GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
-//
-// Inputs:
-//    phi - Latitude of the point, in radians.
-//    lambda - Longitude of the point, in radians.
-//    lambda0 - Longitude of the central meridian to be used, in radians.
-//
-// Outputs:
-//    x - The x coordinate of the computed point.
-//    y - The y coordinate of the computed point.
-//
-// Returns:
-//    The function does not return a value.
-void MapLatLonToXY (double phi, double lambda, double lambda0, double &x, double &y);
-
-// MapXYToLatLon
-// Converts x and y coordinates in the Transverse Mercator projection to
-// a latitude/longitude pair.  Note that Transverse Mercator is not
-// the same as UTM; a scale factor is required to convert between them.
-//
-// Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
-//   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
-//
-// Inputs:
-//   x - The easting of the point, in meters.
-//   y - The northing of the point, in meters.
-//   lambda0 - Longitude of the central meridian to be used, in radians.
-//
-// Outputs:
-//   phi    - Latitude in radians.
-//   lambda - Longitude in radians.
-//
-// Returns:
-//   The function does not return a value.
-//
-// Remarks:
-//   The local variables Nf, nuf2, tf, and tf2 serve the same purpose as
-//   N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect
-//   to the footpoint latitude phif.
-//
-//   x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
-//   to optimize computations.
-void MapXYToLatLon (double x, double y, double lambda0, double& phi, double& lambda);
-
-// LatLonToUTMXY
-// Converts a latitude/longitude pair to x and y coordinates in the
-// Universal Transverse Mercator projection.
-//
-// Inputs:
-//   lat - Latitude of the point, in radians.
-//   lon - Longitude of the point, in radians.
-//   zone - UTM zone to be used for calculating values for x and y.
-//          If zone is less than 1 or greater than 60, the routine
-//          will determine the appropriate zone from the value of lon.
-//
-// Outputs:
-//   x - The x coordinate (easting) of the computed point. (in meters)
-//   y - The y coordinate (northing) of the computed point. (in meters)
-//
-// Returns:
-//   The UTM zone used for calculating the values of x and y.
-int LatLonToUTMXY (double lat, double lon, int zone, double& x, double& y);
-
-// UTMXYToLatLon
-//
-// Converts x and y coordinates in the Universal Transverse Mercator//   The UTM zone parameter should be in the range [1,60].
-
-// projection to a latitude/longitude pair.
-//
-// Inputs:
-// x - The easting of the point, in meters.
-// y - The northing of the point, in meters.
-// zone - The UTM zone in which the point lies.
-// southhemi - True if the point is in the southern hemisphere;
-//               false otherwise.
-//
-// Outputs:
-// lat - The latitude of the point, in radians.
-// lon - The longitude of the point, in radians.
-// 
-// Returns:
-// The function does not return a value.
-void UTMXYToLatLon (double x, double y, int zone, bool southhemi, double& lat, double& lon);
-
-#endif
-