Unverified Commit 4f2b19dd authored by Don Gagne's avatar Don Gagne Committed by GitHub

Merge pull request #7033 from DonLakeFlyer/UTMSHPLoad

Support UTM SHP files
parents 3c45054e 687dde15
...@@ -8,15 +8,22 @@ ...@@ -8,15 +8,22 @@
****************************************************************************/ ****************************************************************************/
#include "SHPFileHelper.h" #include "SHPFileHelper.h"
#include "UTM.h"
#include <QFile> #include <QFile>
#include <QVariant> #include <QVariant>
#include <QtDebug> #include <QtDebug>
#include <QRegularExpression>
const char* SHPFileHelper::_errorPrefix = QT_TR_NOOP("SHP file load failed. %1"); const char* SHPFileHelper::_errorPrefix = QT_TR_NOOP("SHP file load failed. %1");
bool SHPFileHelper::_validateSHPFiles(const QString& shpFile, QString& errorString) /// Validates the specified SHP file is truly a SHP file and is in the format we understand.
/// @param utmZone[out] Zone for UTM shape, 0 for lat/lon shape
/// @param utmSouthernHemisphere[out] true/false for UTM hemisphere
/// @return true: Valid supported SHP file found, false: Invalid or unsupported file found
bool SHPFileHelper::_validateSHPFiles(const QString& shpFile, int* utmZone, bool* utmSouthernHemisphere, QString& errorString)
{ {
*utmZone = 0;
errorString.clear(); errorString.clear();
if (shpFile.endsWith(QStringLiteral(".shp"))) { if (shpFile.endsWith(QStringLiteral(".shp"))) {
...@@ -26,8 +33,25 @@ bool SHPFileHelper::_validateSHPFiles(const QString& shpFile, QString& errorStri ...@@ -26,8 +33,25 @@ bool SHPFileHelper::_validateSHPFiles(const QString& shpFile, QString& errorStri
if (prjFile.open(QIODevice::ReadOnly | QIODevice::Text)) { if (prjFile.open(QIODevice::ReadOnly | QIODevice::Text)) {
QTextStream strm(&prjFile); QTextStream strm(&prjFile);
QString line = strm.readLine(); QString line = strm.readLine();
if (!line.startsWith(QStringLiteral("GEOGCS[\"GCS_WGS_1984\","))) { if (line.startsWith(QStringLiteral("GEOGCS[\"GCS_WGS_1984\","))) {
errorString = QString(_errorPrefix).arg(tr("Only WGS84 projections are supported.")); *utmZone = 0;
*utmSouthernHemisphere = false;
} else if (line.startsWith(QStringLiteral("PROJCS[\"WGS_1984_UTM_Zone_"))) {
QRegularExpression regEx(QStringLiteral("^PROJCS\\[\"WGS_1984_UTM_Zone_(\\d+){1,2}([NS]{1})"));
QRegularExpressionMatch regExMatch = regEx.match(line);
QStringList rgCapture = regExMatch.capturedTexts();
if (rgCapture.count() == 3) {
int zone = rgCapture[1].toInt();
if (zone >= 1 && zone <= 60) {
*utmZone = zone;
*utmSouthernHemisphere = rgCapture[2] == QStringLiteral("S");
}
}
if (*utmZone == 0) {
errorString = QString(_errorPrefix).arg(tr("UTM projection is not in supported format. Must be PROJCS[\"WGS_1984_UTM_Zone_##N/S"));
}
} else {
errorString = QString(_errorPrefix).arg(tr("Only WGS84 or UTM projections are supported."));
} }
} else { } else {
errorString = QString(_errorPrefix).arg(tr("PRJ file open failed: %1").arg(prjFile.errorString())); errorString = QString(_errorPrefix).arg(tr("PRJ file open failed: %1").arg(prjFile.errorString()));
...@@ -42,13 +66,15 @@ bool SHPFileHelper::_validateSHPFiles(const QString& shpFile, QString& errorStri ...@@ -42,13 +66,15 @@ bool SHPFileHelper::_validateSHPFiles(const QString& shpFile, QString& errorStri
return errorString.isEmpty(); return errorString.isEmpty();
} }
SHPHandle SHPFileHelper::_loadShape(const QString& shpFile, QString& errorString) /// @param utmZone[out] Zone for UTM shape, 0 for lat/lon shape
/// @param utmSouthernHemisphere[out] true/false for UTM hemisphere
SHPHandle SHPFileHelper::_loadShape(const QString& shpFile, int* utmZone, bool* utmSouthernHemisphere, QString& errorString)
{ {
SHPHandle shpHandle = Q_NULLPTR; SHPHandle shpHandle = Q_NULLPTR;
errorString.clear(); errorString.clear();
if (_validateSHPFiles(shpFile, errorString)) { if (_validateSHPFiles(shpFile, utmZone, utmSouthernHemisphere, errorString)) {
if (!(shpHandle = SHPOpen(shpFile.toUtf8(), "rb"))) { if (!(shpHandle = SHPOpen(shpFile.toUtf8(), "rb"))) {
errorString = QString(_errorPrefix).arg(tr("SHPOpen failed.")); errorString = QString(_errorPrefix).arg(tr("SHPOpen failed."));
} }
...@@ -63,7 +89,9 @@ ShapeFileHelper::ShapeType SHPFileHelper::determineShapeType(const QString& shpF ...@@ -63,7 +89,9 @@ ShapeFileHelper::ShapeType SHPFileHelper::determineShapeType(const QString& shpF
errorString.clear(); errorString.clear();
SHPHandle shpHandle = SHPFileHelper::_loadShape(shpFile, errorString); int utmZone;
bool utmSouthernHemisphere;
SHPHandle shpHandle = SHPFileHelper::_loadShape(shpFile, &utmZone, &utmSouthernHemisphere, errorString);
if (errorString.isEmpty()) { if (errorString.isEmpty()) {
int cEntities, type; int cEntities, type;
...@@ -85,6 +113,8 @@ ShapeFileHelper::ShapeType SHPFileHelper::determineShapeType(const QString& shpF ...@@ -85,6 +113,8 @@ ShapeFileHelper::ShapeType SHPFileHelper::determineShapeType(const QString& shpF
bool SHPFileHelper::loadPolygonFromFile(const QString& shpFile, QList<QGeoCoordinate>& vertices, QString& errorString) bool SHPFileHelper::loadPolygonFromFile(const QString& shpFile, QList<QGeoCoordinate>& vertices, QString& errorString)
{ {
int utmZone = 0;
bool utmSouthernHemisphere;
double vertexFilterMeters = 5; double vertexFilterMeters = 5;
SHPHandle shpHandle = Q_NULLPTR; SHPHandle shpHandle = Q_NULLPTR;
SHPObject* shpObject = Q_NULLPTR; SHPObject* shpObject = Q_NULLPTR;
...@@ -92,7 +122,7 @@ bool SHPFileHelper::loadPolygonFromFile(const QString& shpFile, QList<QGeoCoordi ...@@ -92,7 +122,7 @@ bool SHPFileHelper::loadPolygonFromFile(const QString& shpFile, QList<QGeoCoordi
errorString.clear(); errorString.clear();
vertices.clear(); vertices.clear();
shpHandle = SHPFileHelper::_loadShape(shpFile, errorString); shpHandle = SHPFileHelper::_loadShape(shpFile, &utmZone, &utmSouthernHemisphere, errorString);
if (!errorString.isEmpty()) { if (!errorString.isEmpty()) {
goto Error; goto Error;
} }
...@@ -111,7 +141,14 @@ bool SHPFileHelper::loadPolygonFromFile(const QString& shpFile, QList<QGeoCoordi ...@@ -111,7 +141,14 @@ bool SHPFileHelper::loadPolygonFromFile(const QString& shpFile, QList<QGeoCoordi
} }
for (int i=0; i<shpObject->nVertices; i++) { for (int i=0; i<shpObject->nVertices; i++) {
vertices.append(QGeoCoordinate(shpObject->padfY[i], shpObject->padfX[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];
}
vertices.append(QGeoCoordinate(lat, lon));
} }
// Filter last vertex such that it differs from first // Filter last vertex such that it differs from first
......
...@@ -29,8 +29,8 @@ public: ...@@ -29,8 +29,8 @@ public:
static bool loadPolygonFromFile(const QString& shpFile, QList<QGeoCoordinate>& vertices, QString& errorString); static bool loadPolygonFromFile(const QString& shpFile, QList<QGeoCoordinate>& vertices, QString& errorString);
private: private:
static bool _validateSHPFiles(const QString& shpFile, QString& errorString); static bool _validateSHPFiles(const QString& shpFile, int* utmZone, bool* utmSouthernHemisphere, QString& errorString);
static SHPHandle _loadShape(const QString& shpFile, QString& errorString); static SHPHandle _loadShape(const QString& shpFile, int* utmZone, bool* utmSouthernHemisphere, QString& errorString);
static const char* _errorPrefix; static const char* _errorPrefix;
}; };
...@@ -15,18 +15,31 @@ ...@@ -15,18 +15,31 @@
// //
// 1) http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html // 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
#include "UTM.h" #include "UTM.h"
#include <math.h>
#define pi 3.14159265358979
/* Ellipsoid model constants (actual values here are for WGS84) */
#define sm_a 6378137.0
#define sm_b 6356752.314
#define sm_EccSquared 6.69437999013e-03
#define UTMScaleFactor 0.9996
// DegToRad // DegToRad
// Converts degrees to radians. // Converts degrees to radians.
FLOAT DegToRad(FLOAT deg) { double DegToRad(double deg) {
return (deg / 180.0 * pi); return (deg / 180.0 * pi);
} }
// RadToDeg // RadToDeg
// Converts radians to degrees. // Converts radians to degrees.
FLOAT RadToDeg(FLOAT rad) { double RadToDeg(double rad) {
return (rad / pi * 180.0); return (rad / pi * 180.0);
} }
...@@ -46,38 +59,38 @@ FLOAT RadToDeg(FLOAT rad) { ...@@ -46,38 +59,38 @@ FLOAT RadToDeg(FLOAT rad) {
// //
// Returns: // Returns:
// The ellipsoidal distance of the point from the equator, in meters. // The ellipsoidal distance of the point from the equator, in meters.
FLOAT ArcLengthOfMeridian (FLOAT phi) { double ArcLengthOfMeridian (double phi) {
FLOAT alpha, beta, gamma, delta, epsilon, n; double alpha, beta, gamma, delta, epsilon, n;
FLOAT result; double result;
/* Precalculate n */ /* Precalculate n */
n = (sm_a - sm_b) / (sm_a + sm_b); n = (sm_a - sm_b) / (sm_a + sm_b);
/* Precalculate alpha */ /* Precalculate alpha */
alpha = ((sm_a + sm_b) / 2.0) alpha = ((sm_a + sm_b) / 2.0)
* (1.0 + (POW(n, 2.0) / 4.0) + (POW(n, 4.0) / 64.0)); * (1.0 + (pow(n, 2.0) / 4.0) + (pow(n, 4.0) / 64.0));
/* Precalculate beta */ /* Precalculate beta */
beta = (-3.0 * n / 2.0) + (9.0 * POW(n, 3.0) / 16.0) beta = (-3.0 * n / 2.0) + (9.0 * pow(n, 3.0) / 16.0)
+ (-3.0 * POW(n, 5.0) / 32.0); + (-3.0 * pow(n, 5.0) / 32.0);
/* Precalculate gamma */ /* Precalculate gamma */
gamma = (15.0 * POW(n, 2.0) / 16.0) gamma = (15.0 * pow(n, 2.0) / 16.0)
+ (-15.0 * POW(n, 4.0) / 32.0); + (-15.0 * pow(n, 4.0) / 32.0);
/* Precalculate delta */ /* Precalculate delta */
delta = (-35.0 * POW(n, 3.0) / 48.0) delta = (-35.0 * pow(n, 3.0) / 48.0)
+ (105.0 * POW(n, 5.0) / 256.0); + (105.0 * pow(n, 5.0) / 256.0);
/* Precalculate epsilon */ /* Precalculate epsilon */
epsilon = (315.0 * POW(n, 4.0) / 512.0); epsilon = (315.0 * pow(n, 4.0) / 512.0);
/* Now calculate the sum of the series and return */ /* Now calculate the sum of the series and return */
result = alpha result = alpha
* (phi + (beta * SIN(2.0 * phi)) * (phi + (beta * sin(2.0 * phi))
+ (gamma * SIN(4.0 * phi)) + (gamma * sin(4.0 * phi))
+ (delta * SIN(6.0 * phi)) + (delta * sin(6.0 * phi))
+ (epsilon * SIN(8.0 * phi))); + (epsilon * sin(8.0 * phi)));
return result; return result;
} }
...@@ -93,9 +106,9 @@ FLOAT ArcLengthOfMeridian (FLOAT phi) { ...@@ -93,9 +106,9 @@ FLOAT ArcLengthOfMeridian (FLOAT phi) {
// Returns: // Returns:
// The central meridian for the given UTM zone, in radians // The central meridian for the given UTM zone, in radians
// Range of the central meridian is the radian equivalent of [-177,+177]. // Range of the central meridian is the radian equivalent of [-177,+177].
FLOAT UTMCentralMeridian(int zone) { double UTMCentralMeridian(int zone) {
FLOAT cmeridian; double cmeridian;
cmeridian = DegToRad(-183.0 + ((FLOAT)zone * 6.0)); cmeridian = DegToRad(-183.0 + ((double)zone * 6.0));
return cmeridian; return cmeridian;
} }
...@@ -115,9 +128,9 @@ FLOAT UTMCentralMeridian(int zone) { ...@@ -115,9 +128,9 @@ FLOAT UTMCentralMeridian(int zone) {
// //
// Returns: // Returns:
// The footpoint latitude, in radians. // The footpoint latitude, in radians.
FLOAT FootpointLatitude(FLOAT y) { double FootpointLatitude(double y) {
FLOAT y_, alpha_, beta_, gamma_, delta_, epsilon_, n; double y_, alpha_, beta_, gamma_, delta_, epsilon_, n;
FLOAT result; double result;
/* Precalculate n (Eq. 10.18) */ /* Precalculate n (Eq. 10.18) */
n = (sm_a - sm_b) / (sm_a + sm_b); n = (sm_a - sm_b) / (sm_a + sm_b);
...@@ -125,31 +138,31 @@ FLOAT FootpointLatitude(FLOAT y) { ...@@ -125,31 +138,31 @@ FLOAT FootpointLatitude(FLOAT y) {
/* Precalculate alpha_ (Eq. 10.22) */ /* Precalculate alpha_ (Eq. 10.22) */
/* (Same as alpha in Eq. 10.17) */ /* (Same as alpha in Eq. 10.17) */
alpha_ = ((sm_a + sm_b) / 2.0) alpha_ = ((sm_a + sm_b) / 2.0)
* (1 + (POW(n, 2.0) / 4) + (POW(n, 4.0) / 64)); * (1 + (pow(n, 2.0) / 4) + (pow(n, 4.0) / 64));
/* Precalculate y_ (Eq. 10.23) */ /* Precalculate y_ (Eq. 10.23) */
y_ = y / alpha_; y_ = y / alpha_;
/* Precalculate beta_ (Eq. 10.22) */ /* Precalculate beta_ (Eq. 10.22) */
beta_ = (3.0 * n / 2.0) + (-27.0 * POW(n, 3.0) / 32.0) beta_ = (3.0 * n / 2.0) + (-27.0 * pow(n, 3.0) / 32.0)
+ (269.0 * POW(n, 5.0) / 512.0); + (269.0 * pow(n, 5.0) / 512.0);
/* Precalculate gamma_ (Eq. 10.22) */ /* Precalculate gamma_ (Eq. 10.22) */
gamma_ = (21.0 * POW(n, 2.0) / 16.0) gamma_ = (21.0 * pow(n, 2.0) / 16.0)
+ (-55.0 * POW(n, 4.0) / 32.0); + (-55.0 * pow(n, 4.0) / 32.0);
/* Precalculate delta_ (Eq. 10.22) */ /* Precalculate delta_ (Eq. 10.22) */
delta_ = (151.0 * POW(n, 3.0) / 96.0) delta_ = (151.0 * pow(n, 3.0) / 96.0)
+ (-417.0 * POW(n, 5.0) / 128.0); + (-417.0 * pow(n, 5.0) / 128.0);
/* Precalculate epsilon_ (Eq. 10.22) */ /* Precalculate epsilon_ (Eq. 10.22) */
epsilon_ = (1097.0 * POW(n, 4.0) / 512.0); epsilon_ = (1097.0 * pow(n, 4.0) / 512.0);
/* Now calculate the sum of the series (Eq. 10.21) */ /* Now calculate the sum of the series (Eq. 10.21) */
result = y_ + (beta_ * SIN(2.0 * y_)) result = y_ + (beta_ * sin(2.0 * y_))
+ (gamma_ * SIN(4.0 * y_)) + (gamma_ * sin(4.0 * y_))
+ (delta_ * SIN(6.0 * y_)) + (delta_ * sin(6.0 * y_))
+ (epsilon_ * SIN(8.0 * y_)); + (epsilon_ * sin(8.0 * y_));
return result; return result;
} }
...@@ -175,24 +188,24 @@ FLOAT FootpointLatitude(FLOAT y) { ...@@ -175,24 +188,24 @@ FLOAT FootpointLatitude(FLOAT y) {
// //
// Returns: // Returns:
// The function does not return a value. // The function does not return a value.
void MapLatLonToXY (FLOAT phi, FLOAT lambda, FLOAT lambda0, FLOAT &x, FLOAT &y) { void MapLatLonToXY (double phi, double lambda, double lambda0, double &x, double &y) {
FLOAT N, nu2, ep2, t, t2, l; double N, nu2, ep2, t, t2, l;
FLOAT l3coef, l4coef, l5coef, l6coef, l7coef, l8coef; double l3coef, l4coef, l5coef, l6coef, l7coef, l8coef;
//FLOAT tmp; // Unused //double tmp; // Unused
/* Precalculate ep2 */ /* Precalculate ep2 */
ep2 = (POW(sm_a, 2.0) - POW(sm_b, 2.0)) / POW(sm_b, 2.0); ep2 = (pow(sm_a, 2.0) - pow(sm_b, 2.0)) / pow(sm_b, 2.0);
/* Precalculate nu2 */ /* Precalculate nu2 */
nu2 = ep2 * POW(COS(phi), 2.0); nu2 = ep2 * pow(cos(phi), 2.0);
/* Precalculate N */ /* Precalculate N */
N = POW(sm_a, 2.0) / (sm_b * SQRT(1 + nu2)); N = pow(sm_a, 2.0) / (sm_b * sqrt(1 + nu2));
/* Precalculate t */ /* Precalculate t */
t = TAN(phi); t = tan(phi);
t2 = t * t; t2 = t * t;
//tmp = (t2 * t2 * t2) - POW(t, 6.0); // Unused //tmp = (t2 * t2 * t2) - pow(t, 6.0); // Unused
/* Precalculate l */ /* Precalculate l */
l = lambda - lambda0; l = lambda - lambda0;
...@@ -216,17 +229,17 @@ void MapLatLonToXY (FLOAT phi, FLOAT lambda, FLOAT lambda0, FLOAT &x, FLOAT &y) ...@@ -216,17 +229,17 @@ void MapLatLonToXY (FLOAT phi, FLOAT lambda, FLOAT lambda0, FLOAT &x, FLOAT &y)
l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2); l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
/* Calculate easting (x) */ /* Calculate easting (x) */
x = N * COS(phi) * l x = N * cos(phi) * l
+ (N / 6.0 * POW(COS(phi), 3.0) * l3coef * POW(l, 3.0)) + (N / 6.0 * pow(cos(phi), 3.0) * l3coef * pow(l, 3.0))
+ (N / 120.0 * POW(COS(phi), 5.0) * l5coef * POW(l, 5.0)) + (N / 120.0 * pow(cos(phi), 5.0) * l5coef * pow(l, 5.0))
+ (N / 5040.0 * POW(COS(phi), 7.0) * l7coef * POW(l, 7.0)); + (N / 5040.0 * pow(cos(phi), 7.0) * l7coef * pow(l, 7.0));
/* Calculate northing (y) */ /* Calculate northing (y) */
y = ArcLengthOfMeridian (phi) y = ArcLengthOfMeridian (phi)
+ (t / 2.0 * N * POW(COS(phi), 2.0) * POW(l, 2.0)) + (t / 2.0 * N * pow(cos(phi), 2.0) * pow(l, 2.0))
+ (t / 24.0 * N * POW(COS(phi), 4.0) * l4coef * POW(l, 4.0)) + (t / 24.0 * N * pow(cos(phi), 4.0) * l4coef * pow(l, 4.0))
+ (t / 720.0 * N * POW(COS(phi), 6.0) * l6coef * POW(l, 6.0)) + (t / 720.0 * N * pow(cos(phi), 6.0) * l6coef * pow(l, 6.0))
+ (t / 40320.0 * N * POW(COS(phi), 8.0) * l8coef * POW(l, 8.0)); + (t / 40320.0 * N * pow(cos(phi), 8.0) * l8coef * pow(l, 8.0));
return; return;
} }
...@@ -260,31 +273,31 @@ void MapLatLonToXY (FLOAT phi, FLOAT lambda, FLOAT lambda0, FLOAT &x, FLOAT &y) ...@@ -260,31 +273,31 @@ void MapLatLonToXY (FLOAT phi, FLOAT lambda, FLOAT lambda0, FLOAT &x, FLOAT &y)
// //
// x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and // x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
// to optimize computations. // to optimize computations.
void MapXYToLatLon (FLOAT x, FLOAT y, FLOAT lambda0, FLOAT& phi, FLOAT& lambda) void MapXYToLatLon (double x, double y, double lambda0, double& phi, double& lambda)
{ {
FLOAT phif, Nf, Nfpow, nuf2, ep2, tf, tf2, tf4, cf; double phif, Nf, Nfpow, nuf2, ep2, tf, tf2, tf4, cf;
FLOAT x1frac, x2frac, x3frac, x4frac, x5frac, x6frac, x7frac, x8frac; double x1frac, x2frac, x3frac, x4frac, x5frac, x6frac, x7frac, x8frac;
FLOAT x2poly, x3poly, x4poly, x5poly, x6poly, x7poly, x8poly; double x2poly, x3poly, x4poly, x5poly, x6poly, x7poly, x8poly;
/* Get the value of phif, the footpoint latitude. */ /* Get the value of phif, the footpoint latitude. */
phif = FootpointLatitude (y); phif = FootpointLatitude (y);
/* Precalculate ep2 */ /* Precalculate ep2 */
ep2 = (POW(sm_a, 2.0) - POW(sm_b, 2.0)) ep2 = (pow(sm_a, 2.0) - pow(sm_b, 2.0))
/ POW(sm_b, 2.0); / pow(sm_b, 2.0);
/* Precalculate cos (phif) */ /* Precalculate cos (phif) */
cf = COS(phif); cf = cos(phif);
/* Precalculate nuf2 */ /* Precalculate nuf2 */
nuf2 = ep2 * POW(cf, 2.0); nuf2 = ep2 * pow(cf, 2.0);
/* Precalculate Nf and initialize Nfpow */ /* Precalculate Nf and initialize Nfpow */
Nf = POW(sm_a, 2.0) / (sm_b * SQRT(1 + nuf2)); Nf = pow(sm_a, 2.0) / (sm_b * sqrt(1 + nuf2));
Nfpow = Nf; Nfpow = Nf;
/* Precalculate tf */ /* Precalculate tf */
tf = TAN(phif); tf = tan(phif);
tf2 = tf * tf; tf2 = tf * tf;
tf4 = tf2 * tf2; tf4 = tf2 * tf2;
...@@ -333,15 +346,15 @@ void MapXYToLatLon (FLOAT x, FLOAT y, FLOAT lambda0, FLOAT& phi, FLOAT& lambda) ...@@ -333,15 +346,15 @@ void MapXYToLatLon (FLOAT x, FLOAT y, FLOAT lambda0, FLOAT& phi, FLOAT& lambda)
/* Calculate latitude */ /* Calculate latitude */
phi = phif + x2frac * x2poly * (x * x) phi = phif + x2frac * x2poly * (x * x)
+ x4frac * x4poly * POW(x, 4.0) + x4frac * x4poly * pow(x, 4.0)
+ x6frac * x6poly * POW(x, 6.0) + x6frac * x6poly * pow(x, 6.0)
+ x8frac * x8poly * POW(x, 8.0); + x8frac * x8poly * pow(x, 8.0);
/* Calculate longitude */ /* Calculate longitude */
lambda = lambda0 + x1frac * x lambda = lambda0 + x1frac * x
+ x3frac * x3poly * POW(x, 3.0) + x3frac * x3poly * pow(x, 3.0)
+ x5frac * x5poly * POW(x, 5.0) + x5frac * x5poly * pow(x, 5.0)
+ x7frac * x7poly * POW(x, 7.0); + x7frac * x7poly * pow(x, 7.0);
return; return;
} }
...@@ -366,9 +379,9 @@ void MapXYToLatLon (FLOAT x, FLOAT y, FLOAT lambda0, FLOAT& phi, FLOAT& lambda) ...@@ -366,9 +379,9 @@ void MapXYToLatLon (FLOAT x, FLOAT y, FLOAT lambda0, FLOAT& phi, FLOAT& lambda)
// //
// Returns: // Returns:
// The UTM zone used for calculating the values of x and y. // The UTM zone used for calculating the values of x and y.
int LatLonToUTMXY (FLOAT lat, FLOAT lon, int zone, FLOAT& x, FLOAT& y) { int LatLonToUTMXY (double lat, double lon, int zone, double& x, double& y) {
if ( (zone < 1) || (zone > 60) ) if ( (zone < 1) || (zone > 60) )
zone = FLOOR((lon + 180.0) / 6) + 1; zone = floor((lon + 180.0) / 6) + 1;
MapLatLonToXY (DegToRad(lat), DegToRad(lon), UTMCentralMeridian(zone), x, y); MapLatLonToXY (DegToRad(lat), DegToRad(lon), UTMCentralMeridian(zone), x, y);
...@@ -401,8 +414,8 @@ int LatLonToUTMXY (FLOAT lat, FLOAT lon, int zone, FLOAT& x, FLOAT& y) { ...@@ -401,8 +414,8 @@ int LatLonToUTMXY (FLOAT lat, FLOAT lon, int zone, FLOAT& x, FLOAT& y) {
// //
// Returns: // Returns:
// The function does not return a value. // The function does not return a value.
void UTMXYToLatLon (FLOAT x, FLOAT y, int zone, bool southhemi, FLOAT& lat, FLOAT& lon) { void UTMXYToLatLon (double x, double y, int zone, bool southhemi, double& lat, double& lon) {
FLOAT cmeridian; double cmeridian;
x -= 500000.0; x -= 500000.0;
x /= UTMScaleFactor; x /= UTMScaleFactor;
......
...@@ -3,9 +3,6 @@ ...@@ -3,9 +3,6 @@
// Original Javascript by Chuck Taylor // Original Javascript by Chuck Taylor
// Port to C++ by Alex Hajnal // Port to C++ by Alex Hajnal
// //
// *** THIS CODE USES 32-BIT FLOATS BY DEFAULT ***
// *** For 64-bit double-precision edit this file: undefine FLOAT_32 and define FLOAT_64 (see below)
//
// This is a simple port of the code on the Geographic/UTM Coordinate Converter (1) page from Javascript to C++. // 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). // 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). // Accuracy seems to be around 50cm (I suspect rounding errors are limiting precision).
...@@ -15,58 +12,18 @@ ...@@ -15,58 +12,18 @@
// //
// 1) http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html // 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 #ifndef UTM_H
#define UTM_H #define UTM_H
// Choose floating point precision:
// 32-bit (for Teensy 3.5/3.6 ARM boards, etc.)
#define FLOAT_64
// 64-bit (for desktop/server use)
//#define FLOAT_64
#ifdef FLOAT_64
#define FLOAT double
#define SIN sin
#define COS cos
#define TAN tan
#define POW pow
#define SQRT sqrt
#define FLOOR floor
#else
#ifdef FLOAT_32
#define FLOAT float
#define SIN sinf
#define COS cosf
#define TAN tanf
#define POW powf
#define SQRT sqrtf
#define FLOOR floorf
#endif
#endif
#include <math.h>
#define pi 3.14159265358979
/* Ellipsoid model constants (actual values here are for WGS84) */
#define sm_a 6378137.0
#define sm_b 6356752.314
#define sm_EccSquared 6.69437999013e-03
#define UTMScaleFactor 0.9996
// DegToRad // DegToRad
// Converts degrees to radians. // Converts degrees to radians.
FLOAT DegToRad(FLOAT deg); double DegToRad(double deg);
// RadToDeg // RadToDeg
// Converts radians to degrees. // Converts radians to degrees.
FLOAT RadToDeg(FLOAT rad); double RadToDeg(double rad);
// ArcLengthOfMeridian // ArcLengthOfMeridian
// Computes the ellipsoidal distance from the equator to a point at a // Computes the ellipsoidal distance from the equator to a point at a
...@@ -84,7 +41,7 @@ FLOAT RadToDeg(FLOAT rad); ...@@ -84,7 +41,7 @@ FLOAT RadToDeg(FLOAT rad);
// //
// Returns: // Returns:
// The ellipsoidal distance of the point from the equator, in meters. // The ellipsoidal distance of the point from the equator, in meters.
FLOAT ArcLengthOfMeridian (FLOAT phi); double ArcLengthOfMeridian (double phi);
// UTMCentralMeridian // UTMCentralMeridian
// Determines the central meridian for the given UTM zone. // Determines the central meridian for the given UTM zone.
...@@ -95,7 +52,7 @@ FLOAT ArcLengthOfMeridian (FLOAT phi); ...@@ -95,7 +52,7 @@ FLOAT ArcLengthOfMeridian (FLOAT phi);
// Returns: // Returns:
// The central meridian for the given UTM zone, in radians // The central meridian for the given UTM zone, in radians
// Range of the central meridian is the radian equivalent of [-177,+177]. // Range of the central meridian is the radian equivalent of [-177,+177].
FLOAT UTMCentralMeridian(int zone); double UTMCentralMeridian(int zone);
// FootpointLatitude // FootpointLatitude
// //
...@@ -110,7 +67,7 @@ FLOAT UTMCentralMeridian(int zone); ...@@ -110,7 +67,7 @@ FLOAT UTMCentralMeridian(int zone);
// //
// Returns: // Returns:
// The footpoint latitude, in radians. // The footpoint latitude, in radians.
FLOAT FootpointLatitude(FLOAT y); double FootpointLatitude(double y);
// MapLatLonToXY // MapLatLonToXY
// Converts a latitude/longitude pair to x and y coordinates in the // Converts a latitude/longitude pair to x and y coordinates in the
...@@ -131,7 +88,7 @@ FLOAT FootpointLatitude(FLOAT y); ...@@ -131,7 +88,7 @@ FLOAT FootpointLatitude(FLOAT y);
// //
// Returns: // Returns:
// The function does not return a value. // The function does not return a value.
void MapLatLonToXY (FLOAT phi, FLOAT lambda, FLOAT lambda0, FLOAT &x, FLOAT &y); void MapLatLonToXY (double phi, double lambda, double lambda0, double &x, double &y);
// MapXYToLatLon // MapXYToLatLon
// Converts x and y coordinates in the Transverse Mercator projection to // Converts x and y coordinates in the Transverse Mercator projection to
...@@ -160,7 +117,7 @@ void MapLatLonToXY (FLOAT phi, FLOAT lambda, FLOAT lambda0, FLOAT &x, FLOAT &y); ...@@ -160,7 +117,7 @@ void MapLatLonToXY (FLOAT phi, FLOAT lambda, FLOAT lambda0, FLOAT &x, FLOAT &y);
// //
// x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and // x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
// to optimize computations. // to optimize computations.
void MapXYToLatLon (FLOAT x, FLOAT y, FLOAT lambda0, FLOAT& phi, FLOAT& lambda); void MapXYToLatLon (double x, double y, double lambda0, double& phi, double& lambda);
// LatLonToUTMXY // LatLonToUTMXY
// Converts a latitude/longitude pair to x and y coordinates in the // Converts a latitude/longitude pair to x and y coordinates in the
...@@ -179,7 +136,7 @@ void MapXYToLatLon (FLOAT x, FLOAT y, FLOAT lambda0, FLOAT& phi, FLOAT& lambda); ...@@ -179,7 +136,7 @@ void MapXYToLatLon (FLOAT x, FLOAT y, FLOAT lambda0, FLOAT& phi, FLOAT& lambda);
// //
// Returns: // Returns:
// The UTM zone used for calculating the values of x and y. // The UTM zone used for calculating the values of x and y.
int LatLonToUTMXY (FLOAT lat, FLOAT lon, int zone, FLOAT& x, FLOAT& y); int LatLonToUTMXY (double lat, double lon, int zone, double& x, double& y);
// UTMXYToLatLon // UTMXYToLatLon
// //
...@@ -200,7 +157,7 @@ int LatLonToUTMXY (FLOAT lat, FLOAT lon, int zone, FLOAT& x, FLOAT& y); ...@@ -200,7 +157,7 @@ int LatLonToUTMXY (FLOAT lat, FLOAT lon, int zone, FLOAT& x, FLOAT& y);
// //
// Returns: // Returns:
// The function does not return a value. // The function does not return a value.
void UTMXYToLatLon (FLOAT x, FLOAT y, int zone, bool southhemi, FLOAT& lat, FLOAT& lon); void UTMXYToLatLon (double x, double y, int zone, bool southhemi, double& lat, double& lon);
#endif #endif
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