Starting on UTMZones.

radixgeo/coordinateconversion.cc

 
 #include "radixgeo/coordinateconversion.hh"
+#include <algorithm>
+
 namespace radix
 {
+const std::array<char, 22> CoordinateConversion::UTMZones::letters = {
+    'A', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K', 'L', 'M',
+    'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Z'};
+
+const std::array<short, 22> CoordinateConversion::UTMZones::degrees = {
+    -90, -84, -72, -64, -56, -48, -40, -32, -24, -16, -8,
+    0,   8,   16,  24,  32,  40,  48,  56,  64,  72,  84};
+
+const std::array<char, 11> CoordinateConversion::UTMZones::neg_letters = {
+    'A', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K', 'L', 'M'};
+
+const std::array<short, 11> CoordinateConversion::UTMZones::neg_degrees = {
+    -90, -84, -72, -64, -56, -48, -40, -32, -24, -16, -8};
+
+const std::array<char, 11> CoordinateConversion::UTMZones::pos_letters = {
+    'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Z'};
+
+const std::array<short, 11> CoordinateConversion::UTMZones::pos_degrees = {
+    0, 8, 16, 24, 32, 40, 48, 56, 64, 72, 84};
+
+short CoordinateConversion::UTMZones::latZoneDegree(char letter)
+{
+  auto itr = std::find(letters.begin(), letters.end(), letter);
+  if (itr == letters.end()) return -100;
+  return *itr;
+}
+
+short CoordinateConversion::UTMZones::lonZone(double longitude)
+{
+  double longZone = 0;
+  if (longitude < 0.0)
+  {
+    longZone = ((180.0 + longitude) / 6) + 1;
+  }
+  else
+  {
+    longZone = (longitude / 6) + 31;
+  }
+  return short(longZone);
+}
+
+char CoordinateConversion::UTMZones::latZone(double latitude)
+{
+  short latIndex = -2;
+  short lat      = short(latitude);
+
+  if (lat >= 0)
+  {
+    size_t len = pos_letters.size();
+    for (size_t i = 0; i < len; i++)
+    {
+      if (lat == pos_degrees[i])
+      {
+        latIndex = i;
+        break;
+      }
+
+      if (lat > pos_degrees[i])
+      {
+        continue;
+      }
+      else
+      {
+        latIndex = i - 1;
+        break;
+      }
+    }
+  }
+  else
+  {
+    size_t len = neg_letters.size();
+    for (size_t i = 0; i < len; i++)
+    {
+      if (lat == neg_degrees[i])
+      {
+        latIndex = i;
+        break;
+      }
+
+      if (lat < neg_degrees[i])
+      {
+        latIndex = i - 1;
+        break;
+      }
+      else
+      {
+        continue;
+      }
+    }
+  }
+
+  if (latIndex == -1)
+  {
+    latIndex = 0;
+  }
+  if (lat >= 0)
+  {
+    if (latIndex == -2)
+    {
+      latIndex = pos_letters.size() - 1;
+    }
+    return pos_letters[latIndex];
+  }
+  else
+  {
+    if (latIndex == -2)
+    {
+      latIndex = neg_letters.size() - 1;
+    }
+    return neg_letters[latIndex];
+  }
+}  // UTMZones::latZone
 }  // namespace radix
\ No newline at end of file

radixgeo/coordinateconversion.hh

@@ -2,13 +2,128 @@
 #define RADIX_RADIXGEO_COORDINATECONVERSION_HH_
 
 #include "radixcore/visibility.hh"
+#include "radixmath/util.hh"
 
+#include <cmath>
+#include <memory>
 #include <sstream>
 #include <stdexcept>
+
 namespace radix
 {
 class RADIX_PUBLIC CoordinateConversion
 {
+ public:
+  struct UTMCoordinate
+  {
+    // earth radius is ~40million meteres
+    // int will cover this
+    int easting;
+    int northing;
+    // 00 -> 60 longitude zones
+    short longitude_zone;
+    // Not officially part of UTM
+    // but frequently used with UTM
+    // C - > X, (A,B,Y,Z) are not used as they cover western and eastern sides
+    // of the Antarctic and Arctic regions
+    char lattitude_zone;
+  };
+
+  class UTMZones
+  {
+    static const std::array<char, 22> letters;
+
+    static const std::array<short, 22> degrees;
+
+    static const std::array<char, 11> neg_letters;
+
+    static const std::array<short, 11> neg_degrees;
+
+    static const std::array<char, 11> pos_letters;
+
+    static const std::array<short, 11> pos_degrees;
+
+   public:
+    static short latZoneDegree(char letter);
+    static short lonZone(double longitude);
+    static char latZone(double latitude);
+  };  // UTMZone
+
+ protected:
+  class LatLon2UTM
+  {
+   private:
+    // equatorial radius
+    double m_equatorialRadius = double(6378137);
+
+    // polar radius
+    double m_polarRadius = double(6356752.314);
+
+    // flattening
+    double m_flattening = double(
+        0.00335281066474748);  // (equatorialRadius-polarRadius)/equatorialRadius;
+
+    // inverse flattening 1/flattening
+    double m_inverseFlattening = double(298.257223563);  // 1/flattening;
+
+    // Mean radius
+    double m_rm =
+        std::pow(m_equatorialRadius * m_polarRadius, double(1.) / double(2.0));
+
+    // scale factor
+    double m_k0 = double(0.9996);
+
+    // eccentricity
+    double m_e = std::sqrt(
+        double(1) - std::pow(m_polarRadius / m_equatorialRadius, double(2.)));
+
+    double m_e1sq = m_e * m_e / (double(1) - m_e * m_e);
+
+    double m_n = (m_equatorialRadius - m_polarRadius) /
+                 (m_equatorialRadius + m_polarRadius);
+
+    // r curv 1
+    double m_rho = double(6368573.744);
+
+    // r curv 2
+    double m_nu = double(6389236.914);
+
+    // Calculate Meridional Arc Length
+    // Meridional Arc
+    double m_S = double(5103266.421);
+
+    double m_A0 = double(6367449.146);
+
+    double m_B0 = double(16038.42955);
+
+    double m_C0 = double(16.83261333);
+
+    double m_D0 = double(0.021984404);
+
+    double m_E0 = double(0.000312705);
+
+    // Calculation Constants
+    // Delta Long
+    double m_p = double(-0.483084);
+
+    double m_sin1 = double(4.84814E-06);
+
+    // Coefficients for UTM Coordinates
+    double m_K1 = double(5101225.115);
+
+    double m_K2 = double(3750.291596);
+
+    double m_K3 = double(1.397608151);
+
+    double m_K4 = double(214839.3105);
+
+    double m_K5 = double(-2.995382942);
+
+    double m_A6 = double(-1.00541E-07);
+
+    void init(double latitude, double longitude);
+  };  // class LatLon2UTM
+
  public:
   /**
    * @brief validate Validates coordinate
@@ -17,8 +132,7 @@ class RADIX_PUBLIC CoordinateConversion
    * Throws std::out_of_range exception if latitude is outside [-90.,90.] or
    * longitude must be [-180,180).
    */
-  template <typename data_type>
-  static void validate(data_type latitude, data_type longitude);
+  static void validate(double latitude, double longitude);
   /**
    * @brief validate Validates coordinate
    * @parm point where point.first is latitutde and point.second is longitude
@@ -26,23 +140,76 @@ class RADIX_PUBLIC CoordinateConversion
    * longitude must be [-180,180).
    */
 
-  template <typename data_type>
-  static void validate(const std::pair<data_type, data_type>& point);
+  static void validate(const std::pair<double, double>& point);
 
 };  // class CoordinateConversion
 
-template <typename data_type>
-void CoordinateConversion::validate(
-    const std::pair<data_type, data_type>& point)
+void CoordinateConversion::LatLon2UTM::init(double latitude, double longitude)
+{
+  latitude = radix::toRadians(latitude);
+  m_rho    = m_equatorialRadius * (double(1.) - m_e * m_e) /
+          std::pow(double(1.) - std::pow(m_e * std::sin(latitude), double(2)),
+                   double(3.) / double(2.0));
+
+  m_nu = m_equatorialRadius /
+         std::pow(double(1) - std::pow(m_e * std::sin(latitude), double(2)),
+                  (double(1) / double(2.0)));
+
+  double var1;
+  if (longitude < double(0.))
+  {
+    var1 = ((int)((180 + longitude) / double(6.))) + 1;
+  }
+  else
+  {
+    var1 = ((int)(longitude / 6)) + 31;
+  }
+  double var2 = (6 * var1) - 183;
+  double var3 = longitude - var2;
+  m_p         = var3 * 3600 / 10000;
+
+  m_S = m_A0 * latitude - m_B0 * std::sin(double(2) * latitude) +
+        m_C0 * std::sin(double(4) * latitude) -
+        m_D0 * std::sin(double(6) * latitude) +
+        m_E0 * std::sin(double(8) * latitude);
+
+  m_K1 = m_S * m_k0;
+  m_K2 = m_nu * std::sin(latitude) * std::cos(latitude) * std::pow(m_sin1, 2) *
+         m_k0 * (100000000) / 2;
+  m_K3 = ((std::pow(m_sin1, double(4)) * m_nu * std::sin(latitude) *
+           std::pow(std::cos(latitude), double(3))) /
+          double(24)) *
+         (double(5) - std::pow(std::tan(latitude), double(2)) +
+          double(9) * m_e1sq * std::pow(std::cos(latitude), double(2)) +
+          double(4) * std::pow(m_e1sq, double(2)) *
+              std::pow(std::cos(latitude), double(4))) *
+         m_k0 * (10000000000000000L);
+
+  m_K4 = m_nu * std::cos(latitude) * m_sin1 * m_k0 * double(10000);
+
+  m_K5 = std::pow(m_sin1 * std::cos(latitude), double(3)) * (m_nu / double(6)) *
+         (double(1) - std::pow(std::tan(latitude), double(2)) +
+          m_e1sq * std::pow(std::cos(latitude), double(2))) *
+         m_k0 * 1000000000000L;
+
+  m_A6 = (std::pow(m_p * m_sin1, double(6)) * m_nu * std::sin(latitude) *
+          std::pow(std::cos(latitude), double(5)) / double(720)) *
+         (double(61) - double(58) * std::pow(std::tan(latitude), double(2)) +
+          std::pow(std::tan(latitude), double(4)) +
+          double(270) * m_e1sq * std::pow(std::cos(latitude), double(2)) -
+          double(330) * m_e1sq * std::pow(std::sin(latitude), double(2))) *
+         m_k0 * (1E+24);
+}  // LatLon2UTM init
+
+void CoordinateConversion::validate(const std::pair<double, double>& point)
 {
   CoordinateConversion::validate(point.first, point.second);
 }
 
-template <typename data_type>
-void CoordinateConversion::validate(data_type latitude, data_type longitude)
+void CoordinateConversion::validate(double latitude, double longitude)
 {
-  if (latitude < data_type(-90.0) || latitude > data_type(90.0) ||
-      longitude < data_type(-180.0) || longitude >= data_type(180.0))
+  if (latitude < double(-90.0) || latitude > double(90.0) ||
+      longitude < double(-180.0) || longitude >= double(180.0))
   {
     std::ostringstream oss;
     oss << "Invalid coordinate [" << latitude << "," << longitude

radixgeo/tests/tstCoordinateConversion.cc

@@ -121,4 +121,29 @@ TEST(Radixgeo, CoordinateRange)
   }
 }
 
-TEST(Radixgeo, CoordinateConversion) {}
+TEST(Radixgeo, UTMZones)
+{
+  // latitude
+  EXPECT_EQ('N', CoordinateConversion::UTMZones::latZone(0.));
+  EXPECT_EQ('N', CoordinateConversion::UTMZones::latZone(0.13));
+  EXPECT_EQ('G', CoordinateConversion::UTMZones::latZone(-45.6456));
+  EXPECT_EQ('L', CoordinateConversion::UTMZones::latZone(-12.7650));
+  EXPECT_EQ('C', CoordinateConversion::UTMZones::latZone(-80.5434));
+  EXPECT_EQ('Z', CoordinateConversion::UTMZones::latZone(90.0000));
+  EXPECT_EQ('A', CoordinateConversion::UTMZones::latZone(-90.0000));
+  EXPECT_EQ('Q', CoordinateConversion::UTMZones::latZone(23.4578));
+  EXPECT_EQ('X', CoordinateConversion::UTMZones::latZone(77.3450));
+  EXPECT_EQ('A', CoordinateConversion::UTMZones::latZone(-89.3454));
+  // longitude
+  EXPECT_EQ(31, CoordinateConversion::UTMZones::lonZone(0.));
+  EXPECT_EQ(30, CoordinateConversion::UTMZones::lonZone(-0.2324));
+  EXPECT_EQ(34, CoordinateConversion::UTMZones::lonZone(23.3545));
+  EXPECT_EQ(25, CoordinateConversion::UTMZones::lonZone(-33.8765));
+  EXPECT_EQ(02, CoordinateConversion::UTMZones::lonZone(-170.6540));
+  EXPECT_EQ(60, CoordinateConversion::UTMZones::lonZone(177.0000));
+  EXPECT_EQ(01, CoordinateConversion::UTMZones::lonZone(-177.0000));
+  EXPECT_EQ(31, CoordinateConversion::UTMZones::lonZone(3.0000));
+  EXPECT_EQ(8, CoordinateConversion::UTMZones::lonZone(-135.4545));
+  EXPECT_EQ(57, CoordinateConversion::UTMZones::lonZone(156.9876));
+  EXPECT_EQ(22, CoordinateConversion::UTMZones::lonZone(-48.9306));
+}