Commit 6a35688d authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Merge branch 'coordinate-conversion' into 'master'

Coordinate conversion

See merge request !35
parents e675fb1c 58fc88a0
Pipeline #13217 passed with stages
in 13 minutes and 59 seconds
......@@ -76,7 +76,7 @@ linux_analysis:
- which cmake
- module load valgrind
- cmake -DDEBUG_OUTPUT=1 -D COVERAGE_EXTRA_FLAGS="-s ${CI_PROJECT_DIR}/googletest -d" -D MEMCHECK_COMMAND=$(which valgrind) -DBUILDNAME=$(uname -s)-GCC-4.8.5-Debug-${CI_BUILD_REF_NAME} -DCMAKE_BUILD_TYPE=DEBUG -Dradix_ENABLE_COVERAGE_TESTING=ON -Dradix_ENABLE_TESTS=ON -Dradix_ENABLE_SECONDARY_TESTED_CODE=ON -Dradix_ENABLE_TESTS=ON -Dradix_ENABLE_radixplot=OFF -Dradix_ENABLE_radixwidgets=OFF ..
- ctest -D ExperimentalStart -D ExperimentalBuild -D ExperimentalTest -D ExperimentalMemCheck -D ExperimentalSubmit
- ctest -D ExperimentalStart -D ExperimentalBuild -D ExperimentalTest -D ExperimentalMemCheck -D ExperimentalCoverage -D ExperimentalSubmit
allow_failure: true
linux_openmpi_testing:
......
......@@ -8,6 +8,7 @@ TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
para radixpara SS OPTIONAL
command radixcommand SS OPTIONAL
math radixmath SS OPTIONAL
geo radixgeo SS OPTIONAL
glls radixglls SS OPTIONAL
io radixio SS OPTIONAL
geometry radixgeometry SS OPTIONAL
......
TRIBITS_SUBPACKAGE(geo)
SET(SOURCE
coordinateconversion.cc
)
SET(HEADERS
coordinateconversion.hh
)
#
# Add library
TRIBITS_ADD_LIBRARY(radixgeo
SOURCES ${SOURCE}
NOINSTALLHEADERS ${HEADERS}
)
#
# Add testing directory
TRIBITS_ADD_TEST_DIRECTORIES(tests)
INSTALL(FILES ${HEADERS} DESTINATION "include/radixgeo/")
TRIBITS_SUBPACKAGE_POSTPROCESS()
TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
LIB_REQUIRED_PACKAGES radixcore radixbug
LIB_OPTIONAL_PACKAGES
TEST_REQUIRED_PACKAGES testframework googletest
TEST_OPTIONAL_PACKAGES
LIB_REQUIRED_TPLS
LIB_OPTIONAL_TPLS
TEST_REQUIRED_TPLS
TEST_OPTIONAL_TPLS
)
#include "radixgeo/coordinateconversion.hh"
#include "radixmath/constants.hh"
#include <algorithm>
namespace radix
{
const std::array<char, 22> 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> 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> UTMZones::neg_letters = {
{'A', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K', 'L', 'M'}};
const std::array<short, 11> UTMZones::neg_degrees = {
{-90, -84, -72, -64, -56, -48, -40, -32, -24, -16, -8}};
const std::array<char, 11> UTMZones::pos_letters = {
{'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Z'}};
const std::array<short, 11> UTMZones::pos_degrees = {
{0, 8, 16, 24, 32, 40, 48, 56, 64, 72, 84}};
const std::array<char, 11> UTM2LatLon::southern_hemisphere = {
{'A', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'L', 'L', 'M'}};
short UTMZones::latZoneDegree(char letter)
{
auto itr = std::find(letters.begin(), letters.end(), letter);
if (itr == letters.end()) return -100;
return *itr;
}
short 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 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];
}
}
double 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 LatLon2UTM::easting() const
{
return 500000. + (m_K4 * m_p + m_K5 * std::pow(m_p, 3.));
}
UTMCoordinate 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 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);
}
std::string CoordinateConversion::toString(const UTMCoordinate& utm)
{
char buff[120];
sprintf(buff, "%02d %c %f %f", utm.longitude_zone, utm.lattitude_zone,
utm.easting, utm.northing);
return std::string(buff);
}
void CoordinateConversion::validate(double latitude, double longitude)
{
radix_tagged_line("Validating (" << latitude << "," << 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());
}
}
void UTM2LatLon::init(const 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 UTM2LatLon::hemisphere(char latitude_zone)
{
auto it = std::find(std::begin(southern_hemisphere),
std::end(southern_hemisphere), latitude_zone);
if (it == southern_hemisphere.end())
{
return 'N';
}
else
{
return 'S';
}
}
std::pair<double, double> UTM2LatLon::toLatLon(const 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
#ifndef RADIX_RADIXGEO_COORDINATECONVERSION_HH_
#define RADIX_RADIXGEO_COORDINATECONVERSION_HH_
#include "radixcore/visibility.hh"
#include "radixmath/util.hh"
#include <array>
#include <cmath>
#include <memory>
#include <sstream>
#include <stdexcept>
namespace radix
{
struct RADIX_PUBLIC UTMCoordinate
{
// earth radius is ~40million meteres
// int will cover this
double easting;
double 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 RADIX_PUBLIC CoordinateConversion
{
public:
/**
* @brief validate Validates coordinate
* @param latitude latitude coordinate
* @param longitude longitude coordinate
* Throws std::out_of_range exception if latitude is outside [-90.,90.] or
* longitude must be [-180,180).
*/
static void validate(double latitude, double longitude);
/**
* @brief validate Validates coordinate
* @param point where point.first is latitutde and point.second is longitude
* Throws std::out_of_range exception if latitude is outside [-90.,90.] or
* longitude must be [-180,180).
*/
static void validate(const std::pair<double, double>& point);
/**
* @brief toString Converts UTM coordinate to string representation
* @param utm
* @return std::string "11 S 400000 5500"
*/
static std::string toString(const UTMCoordinate& utm);
}; // class CoordinateConversion
class RADIX_PUBLIC 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:
/**
* @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
class RADIX_PUBLIC 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);
/**
* @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
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;