Finished UTM2LatLon coordinate conversion.

radixgeo/coordinateconversion.cc

 
 #include "radixgeo/coordinateconversion.hh"
+#include "radixmath/constants.hh"
 
 #include <algorithm>
 
@@ -167,7 +168,7 @@ void CoordinateConversion::LatLon2UTM::init(double latitude, double longitude)
   }
   else
   {
-    var1 = ((int)(longitude / 6)) + 31;
+    var1 = ((int)(longitude / 6.)) + 31;
   }
   double var2 = (6 * var1) - 183;
   double var3 = longitude - var2;
@@ -225,6 +226,55 @@ void CoordinateConversion::validate(double latitude, double longitude)
   }
 }
 
+void CoordinateConversion::UTM2LatLon::init(
+    const CoordinateConversion::UTMCoordinate& utm)
+{
+  m_arc = utm.northing / m_k0;
+  m_mu  = m_arc /
+         (m_a * (1. - std::pow(m_e, 2.) / 4.0 - 3. * std::pow(m_e, 4.) / 64.0 -
+                 5. * std::pow(m_e, 6.) / 256.0));
+
+  m_ei = (1. - std::pow((1. - m_e * m_e), (1. / 2.0))) /
+         (1. + std::pow((1. - m_e * m_e), (1. / 2.0)));
+
+  m_ca = 3. * m_ei / 2. - 27. * std::pow(m_ei, 3.) / 32.0;
+
+  m_cb   = 21. * std::pow(m_ei, 2.) / 16. - 55. * std::pow(m_ei, 4.) / 32.;
+  m_cc   = 151. * std::pow(m_ei, 3.) / 96.;
+  m_cd   = 1097. * std::pow(m_ei, 4.) / 512.;
+  m_phi1 = m_mu + m_ca * std::sin(2. * m_mu) + m_cb * std::sin(4. * m_mu) +
+           m_cc * std::sin(6. * m_mu) + m_cd * std::sin(8. * m_mu);
+
+  m_n0 =
+      m_a / std::pow((1. - std::pow((m_e * std::sin(m_phi1)), 2.)), (1. / 2.0));
+
+  m_r0 = m_a * (1. - m_e * m_e) /
+         std::pow((1. - std::pow((m_e * std::sin(m_phi1)), 2.)), (3. / 2.0));
+  m_fact1 = m_n0 * std::tan(m_phi1) / m_r0;
+
+  m_a1    = 500000. - utm.easting;
+  m_dd0   = m_a1 / (m_n0 * m_k0);
+  m_fact2 = m_dd0 * m_dd0 / 2.;
+
+  m_t0    = std::pow(std::tan(m_phi1), 2.);
+  m_Q0    = m_e1sq * std::pow(std::cos(m_phi1), 2.);
+  m_fact3 = (5. + 3. * m_t0 + 10. * m_Q0 - 4. * m_Q0 * m_Q0 - 9. * m_e1sq) *
+            std::pow(m_dd0, 4.) / 24.;
+
+  m_fact4 = (61. + 90. * m_t0 + 298. * m_Q0 + 45. * m_t0 * m_t0 -
+             252. * m_e1sq - 3. * m_Q0 * m_Q0) *
+            std::pow(m_dd0, 6.) / 720.;
+
+  //
+  m_lof1 = m_a1 / (m_n0 * m_k0);
+  m_lof2 = (1. + 2. * m_t0 + m_Q0) * std::pow(m_dd0, 3.) / 6.0;
+  m_lof3 = (5. - 2. * m_Q0 + 28. * m_t0 - 3. * std::pow(m_Q0, 2.) +
+            8. * m_e1sq + 24. * std::pow(m_t0, 2.)) *
+           std::pow(m_dd0, 5.) / 120.;
+  m_a2 = (m_lof1 - m_lof2 + m_lof3) / std::cos(m_phi1);
+  m_a3 = m_a2 * 180. / radix::PI;
+}
+
 char CoordinateConversion::UTM2LatLon::hemisphere(char latitude_zone)
 {
   auto it = std::find(std::begin(southern_hemisphere),
@@ -239,5 +289,40 @@ char CoordinateConversion::UTM2LatLon::hemisphere(char latitude_zone)
   }
 }
 
+std::pair<double, double> CoordinateConversion::UTM2LatLon::toLatLon(
+    const CoordinateConversion::UTMCoordinate& c)
+{
+  UTMCoordinate utm   = c;
+  double latitude     = 0.;
+  double longitude    = 0.;
+  char utm_hemisphere = hemisphere(utm.lattitude_zone);
+  // modify for Southern hemisphere
+  if ('S' == utm_hemisphere)
+  {
+    utm.northing = 10000000 - utm.northing;
+  }
+  // initialize internal memory
+  init(utm);
+  latitude =
+      180. * (m_phi1 - m_fact1 * (m_fact2 + m_fact3 + m_fact4)) / radix::PI;
+
+  if (utm.longitude_zone > 0)
+  {
+    m_zoneCM = 6 * utm.longitude_zone - 183.0;
+  }
+  else
+  {
+    m_zoneCM = 3.0;
+  }
+
+  longitude = m_zoneCM - m_a3;
+  if ('S' == utm_hemisphere)
+  {
+    latitude = -latitude;
+  }
+
+  return std::make_pair(latitude, longitude);
+}
+
 // UTMZones::latZone
 }  // namespace radix
\ No newline at end of file

radixgeo/coordinateconversion.hh

@@ -19,8 +19,8 @@ class RADIX_PUBLIC CoordinateConversion
   {
     // earth radius is ~40million meteres
     // int will cover this
-    int easting;
-    int northing;
+    double easting;
+    double northing;
     // 00 -> 60 longitude zones
     short longitude_zone;
     // Not officially part of UTM
@@ -153,9 +153,70 @@ class RADIX_PUBLIC CoordinateConversion
   class RADIX_PUBLIC UTM2LatLon
   {
    private:
+    double m_arc;
+
+    double m_mu;
+
+    double m_ei;
+
+    double m_ca;
+
+    double m_cb;
+
+    double m_cc;
+
+    double m_cd;
+
+    double m_n0;
+
+    double m_r0;
+
+    double m_a1;
+
+    double m_dd0;
+
+    double m_t0;
+
+    double m_Q0;
+
+    double m_lof1;
+
+    double m_lof2;
+
+    double m_lof3;
+
+    double m_a2;
+
+    double m_phi1;
+
+    double m_fact1;
+
+    double m_fact2;
+
+    double m_fact3;
+
+    double m_fact4;
+
+    double m_zoneCM;
+
+    double m_a3;
+
+    double m_b = 6356752.314;
+
+    double m_a = 6378137;
+
+    double m_e = 0.081819191;
+
+    double m_e1sq = 0.006739497;
+
+    double m_k0 = 0.9996;
+
+    void init(const UTMCoordinate& utm);
+
    public:
     static const std::array<char, 11> southern_hemisphere;
     static char hemisphere(char latitude_zone);
+    std::pair<double, double> toLatLon(const UTMCoordinate& c);
   };  // class UTM2LatLon
 
   /**

radixgeo/tests/tstCoordinateConversion.cc

@@ -150,86 +150,148 @@ TEST(Radixgeo, UTMZones)
 
 TEST(Radixgeo, LatLon2UTM)
 {
+  double fuzzy = 1;
   CoordinateConversion::LatLon2UTM conv;
   CoordinateConversion::UTMCoordinate utm = conv.toUTM(0.0000, 0.0000);
-  utm                                     = conv.toUTM(0.0000, 0.0000);
   EXPECT_EQ(31, utm.longitude_zone);
   EXPECT_EQ('N', utm.lattitude_zone);
-  EXPECT_EQ(166021, utm.easting);
-  EXPECT_EQ(0, utm.northing);
+  EXPECT_NEAR(166021.329, utm.easting, fuzzy);
+  EXPECT_NEAR(0., utm.northing, fuzzy);
 
   utm = conv.toUTM(0.1300, -0.2324);
   EXPECT_EQ(30, utm.longitude_zone);
   EXPECT_EQ('N', utm.lattitude_zone);
-  EXPECT_EQ(808084, utm.easting);
-  EXPECT_EQ(14385, utm.northing);
+  EXPECT_NEAR(808084.568, utm.easting, fuzzy);
+  EXPECT_NEAR(14385.8, utm.northing, fuzzy);
 
   utm = conv.toUTM(-45.6456, 23.3545);
   EXPECT_EQ(34, utm.longitude_zone);
   EXPECT_EQ('G', utm.lattitude_zone);
-  EXPECT_EQ(683473, utm.easting);
-  EXPECT_EQ(4942631, utm.northing);
+  EXPECT_NEAR(683473.881, utm.easting, fuzzy);
+  EXPECT_NEAR(4942631.186, utm.northing, fuzzy);
 
   utm = conv.toUTM(-12.7650, -33.8765);
   EXPECT_EQ(25, utm.longitude_zone);
   EXPECT_EQ('L', utm.lattitude_zone);
-  EXPECT_EQ(404859, utm.easting);
-  EXPECT_EQ(8588690, utm.northing);
+  EXPECT_NEAR(404859.077, utm.easting, fuzzy);
+  EXPECT_NEAR(8588690.973, utm.northing, fuzzy);
 
   utm = conv.toUTM(-80.5434, -170.6540);
   EXPECT_EQ(2, utm.longitude_zone);
   EXPECT_EQ('C', utm.lattitude_zone);
-  EXPECT_EQ(506346, utm.easting);
-  EXPECT_EQ(1057742, utm.northing);
+  EXPECT_NEAR(506346.425, utm.easting, fuzzy);
+  EXPECT_NEAR(1057742.538, utm.northing, fuzzy);
 
   utm = conv.toUTM(90.0000, 177.0000);
   EXPECT_EQ(60, utm.longitude_zone);
   EXPECT_EQ('Z', utm.lattitude_zone);
-  EXPECT_EQ(500000, utm.easting);
-  EXPECT_EQ(9997964, utm.northing);
+  EXPECT_NEAR(500000., utm.easting, fuzzy);
+  EXPECT_NEAR(9997964.943, utm.northing, fuzzy);
 
   utm = conv.toUTM(-90.0000, -177.0000);
   EXPECT_EQ(1, utm.longitude_zone);
   EXPECT_EQ('A', utm.lattitude_zone);
-  EXPECT_EQ(500000, utm.easting);
-  EXPECT_EQ(2035, utm.northing);
+  EXPECT_NEAR(500000, utm.easting, fuzzy);
+  EXPECT_NEAR(2035.06, utm.northing, fuzzy);
 
   utm = conv.toUTM(90.0000, 3.0000);
   EXPECT_EQ(31, utm.longitude_zone);
   EXPECT_EQ('Z', utm.lattitude_zone);
-  EXPECT_EQ(500000, utm.easting);
-  EXPECT_EQ(9997964, utm.northing);
+  EXPECT_NEAR(500000., utm.easting, fuzzy);
+  EXPECT_NEAR(9997964.943, utm.northing, fuzzy);
 
   utm = conv.toUTM(23.4578, -135.4545);
   EXPECT_EQ(8, utm.longitude_zone);
   EXPECT_EQ('Q', utm.lattitude_zone);
-  EXPECT_EQ(453580, utm.easting);
-  EXPECT_EQ(2594272, utm.northing);
+  EXPECT_NEAR(453580.338, utm.easting, fuzzy);
+  EXPECT_NEAR(2594272.912, utm.northing, fuzzy);
 
   utm = conv.toUTM(77.3450, 156.9876);
   EXPECT_EQ(57, utm.longitude_zone);
   EXPECT_EQ('X', utm.lattitude_zone);
-  EXPECT_EQ(450793, utm.easting);
-  EXPECT_EQ(8586116, utm.northing);
+  EXPECT_NEAR(450793.55, utm.easting, fuzzy);
+  EXPECT_NEAR(8586116.23, utm.northing, fuzzy);
 
   utm = conv.toUTM(-89.3454, -48.9306);
   EXPECT_EQ(22, utm.longitude_zone);
   EXPECT_EQ('A', utm.lattitude_zone);
-  EXPECT_EQ(502639, utm.easting);
-  EXPECT_EQ(75072, utm.northing);
+  EXPECT_NEAR(502639.064, utm.easting, fuzzy);
+  EXPECT_NEAR(75073., utm.northing, fuzzy);
 
   utm = conv.toUTM(37.131467, -116.039417);
   EXPECT_EQ(11, utm.longitude_zone);
   EXPECT_EQ('S', utm.lattitude_zone);
-  EXPECT_EQ(585322, utm.easting);
-  EXPECT_EQ(4109888, utm.northing);
-  CoordinateConversion::UTMZones::latZoneDegree(utm.lattitude_zone);
+  EXPECT_NEAR(585322.62, utm.easting, fuzzy);
+  EXPECT_NEAR(4109888.360, utm.northing, fuzzy);
 }
 
 TEST(Radixgeo, UTM2LatLon)
 {
+  double fuzzy = 1e-2;
   EXPECT_EQ('S', CoordinateConversion::UTM2LatLon::hemisphere('A'));
   EXPECT_EQ('S', CoordinateConversion::UTM2LatLon::hemisphere('M'));
   EXPECT_EQ('N', CoordinateConversion::UTM2LatLon::hemisphere('N'));
   EXPECT_EQ('N', CoordinateConversion::UTM2LatLon::hemisphere('Z'));
+
+  CoordinateConversion::LatLon2UTM latLonConv;
+  CoordinateConversion::UTM2LatLon utmConv;
+  CoordinateConversion::UTMCoordinate utm = latLonConv.toUTM(0.0000, 0.0000);
+  std::pair<double, double> latLong       = utmConv.toLatLon(utm);
+  EXPECT_NEAR(0., latLong.first, fuzzy);
+  EXPECT_NEAR(0., latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(0.1300, -0.2324);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(0.1300, latLong.first, fuzzy);
+  EXPECT_NEAR(-0.2324, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(-45.6456, 23.3545);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(-45.6456, latLong.first, fuzzy);
+  EXPECT_NEAR(23.3545, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(-12.7650, -33.8765);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(-12.7650, latLong.first, fuzzy);
+  EXPECT_NEAR(-33.8765, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(-80.5434, -170.6540);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(-80.5434, latLong.first, fuzzy);
+  EXPECT_NEAR(-170.6540, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(90.0000, 177.0000);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(90.0000, latLong.first, fuzzy);
+  EXPECT_NEAR(177.0000, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(-90.0000, -177.0000);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(-90.0000, latLong.first, fuzzy);
+  EXPECT_NEAR(-177.0000, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(90.0000, 3.0000);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(90.0000, latLong.first, fuzzy);
+  EXPECT_NEAR(3.0000, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(23.4578, -135.4545);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(23.4578, latLong.first, fuzzy);
+  EXPECT_NEAR(-135.4545, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(77.3450, 156.9876);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(77.3450, latLong.first, fuzzy);
+  EXPECT_NEAR(156.9876, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(-89.3454, -48.9306);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(-89.3454, latLong.first, fuzzy);
+  EXPECT_NEAR(-48.9306, latLong.second, fuzzy);
+
+  utm     = latLonConv.toUTM(37.131467, -116.039417);
+  latLong = utmConv.toLatLon(utm);
+  EXPECT_NEAR(37.131467, latLong.first, fuzzy);
+  EXPECT_NEAR(-116.039417, latLong.second, fuzzy);
 }