Skip to content

GitLab

Commit 7dbe0e5c authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Initial version of LatLon2UTM conversion.

parent 48629f4a
Pipeline #12888 failed with stages
in 4 minutes and 13 seconds
Hide whitespace changes
Inline Side-by-side
radixgeo/coordinateconversion.cc
View file @ 7dbe0e5c
......@@ -115,5 +115,110 @@ char CoordinateConversion::UTMZones::latZone(double latitude)
}
return neg_letters[latIndex];
}
} // UTMZones::latZone
}
double CoordinateConversion::LatLon2UTM::northing(double latitude) const
{
double northing = m_K1 + m_K2 * m_p * m_p + m_K3 * std::pow(m_p, 4.);
if (latitude < 0.0)
{
northing = 10000000. + northing;
}
return northing;
}
double CoordinateConversion::LatLon2UTM::easting() const
{
return 500000. + (m_K4 * m_p + m_K5 * std::pow(m_p, 3.));
}
CoordinateConversion::UTMCoordinate CoordinateConversion::LatLon2UTM::toUTM(
double latitude, double longitude)
{
CoordinateConversion::validate(latitude, longitude);
UTMCoordinate utm;
init(latitude, longitude);
utm.longitude_zone = UTMZones::lonZone(longitude);
utm.lattitude_zone = UTMZones::latZone(latitude);
utm.easting = easting();
utm.northing = northing(latitude);
return utm;
}
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);
}
void CoordinateConversion::validate(double latitude, double longitude)
{
if (latitude < -90.0 || latitude > 90.0 || longitude < -180.0 ||
longitude >= 180.0)
{
std::ostringstream oss;
oss << "Invalid coordinate [" << latitude << "," << longitude
<< "].\nLatitude must be [-90,90] and longitude must be "
"[-180,180).";
throw std::out_of_range(oss.str());
}
}
// UTMZones::latZone
} // namespace radix
\ No newline at end of file
radixgeo/coordinateconversion.hh
View file @ 7dbe0e5c
......@@ -4,6 +4,7 @@
#include "radixcore/visibility.hh"
#include "radixmath/util.hh"
#include <array>
#include <cmath>
#include <memory>
#include <sstream>
......@@ -44,12 +45,27 @@ class RADIX_PUBLIC CoordinateConversion
static const std::array<short, 11> pos_degrees;
public:
/**
* @brief latZoneDegree Latitude
* @param letter 'A' thru 'Z' zone
* @return short -90 thru 90 latitude zone offset
*/
static short latZoneDegree(char letter);
/**
* @brief lonZone Longitude zone
* @param longitude degrees
* @return short 0 thru 60 zone
*/
static short lonZone(double longitude);
/**
* @brief latZone Latitude zone
* @param latitude degrees
* @return 'A' thru 'Z' zone
*/
static char latZone(double latitude);
}; // UTMZone
protected:
class LatLon2UTM
{
private:
......@@ -121,10 +137,19 @@ class RADIX_PUBLIC CoordinateConversion
double m_A6 = double(-1.00541E-07);
/**
* @brief init initializes internal data necessary for coordinate conversion
* @param latitude degrees
* @param longitude degrees
*/
void init(double latitude, double longitude);
public:
double northing(double latitude) const;
double easting() const;
UTMCoordinate toUTM(double latitude, double longitude);
}; // class LatLon2UTM
public:
/**
* @brief validate Validates coordinate
* @param latitude latitude coordinate
......@@ -144,79 +169,5 @@ class RADIX_PUBLIC CoordinateConversion
}; // class CoordinateConversion
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);
}
void CoordinateConversion::validate(double latitude, double longitude)
{
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
<< "].\nLatitude must be [-90,90] and longitude must be "
"[-180,180).";
throw std::out_of_range(oss.str());
}
}
} // namespace radix
#endif /** RADIX_RADIXGEO_COORDINATECONVERSION_HH_ */
\ No newline at end of file
radixgeo/tests/tstCoordinateConversion.cc
View file @ 7dbe0e5c
......@@ -147,3 +147,26 @@ TEST(Radixgeo, UTMZones)
EXPECT_EQ(57, CoordinateConversion::UTMZones::lonZone(156.9876));
EXPECT_EQ(22, CoordinateConversion::UTMZones::lonZone(-48.9306));
}
TEST(Radixgeo, LatLon2UTM)
{
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);
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);
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);
}
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