From b3e8081942cd6d80bafacbde68e44731609d9b3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20T=C5=AFma?= Date: Wed, 29 Mar 2017 00:17:47 +0200 Subject: [PATCH] Added support for Transversal Mercator projection --- gpxsee.pro | 9 +++- src/coordinates.cpp | 10 ---- src/coordinates.h | 7 +-- src/emptymap.cpp | 7 +-- src/mercator.cpp | 13 +++++ src/mercator.h | 13 +++++ src/offlinemap.cpp | 69 +++++++++++++++++++++++---- src/offlinemap.h | 7 ++- src/onlinemap.cpp | 7 +-- src/projection.h | 15 ++++++ src/transversemercator.cpp | 97 ++++++++++++++++++++++++++++++++++++++ src/transversemercator.h | 22 +++++++++ src/wgs84.h | 3 +- 13 files changed, 244 insertions(+), 35 deletions(-) create mode 100644 src/mercator.cpp create mode 100644 src/mercator.h create mode 100644 src/projection.h create mode 100644 src/transversemercator.cpp create mode 100644 src/transversemercator.h diff --git a/gpxsee.pro b/gpxsee.pro index 40c72107..7bbda49e 100644 --- a/gpxsee.pro +++ b/gpxsee.pro @@ -84,7 +84,10 @@ HEADERS += src/config.h \ src/mapdir.h \ src/matrix.h \ src/tar.h \ - src/atlas.h + src/atlas.h \ + src/projection.h \ + src/mercator.h \ + src/transversemercator.h SOURCES += src/main.cpp \ src/gui.cpp \ src/poi.cpp \ @@ -144,7 +147,9 @@ SOURCES += src/main.cpp \ src/mapdir.cpp \ src/matrix.cpp \ src/tar.cpp \ - src/atlas.cpp + src/atlas.cpp \ + src/mercator.cpp \ + src/transversemercator.cpp RESOURCES += gpxsee.qrc TRANSLATIONS = lang/gpxsee_cs.ts \ lang/gpxsee_sv.ts diff --git a/src/coordinates.cpp b/src/coordinates.cpp index 02b78396..436afea4 100644 --- a/src/coordinates.cpp +++ b/src/coordinates.cpp @@ -17,16 +17,6 @@ qreal Coordinates::distanceTo(const Coordinates &c) const return (WGS84_RADIUS * (2.0 * atan2(sqrt(a), sqrt(1.0 - a)))); } -QPointF Coordinates::toMercator() const -{ - return QPointF(_lon, rad2deg(log(tan(M_PI/4.0 + deg2rad(_lat)/2.0)))); -} - -Coordinates Coordinates::fromMercator(const QPointF &m) -{ - return Coordinates(m.x(), rad2deg(2 * atan(exp(deg2rad(m.y()))) - M_PI/2)); -} - QDebug operator<<(QDebug dbg, const Coordinates &coordinates) { dbg.nospace() << "Coordinates(" << coordinates.lon() << ", " diff --git a/src/coordinates.h b/src/coordinates.h index 1251fe90..3a7db43e 100644 --- a/src/coordinates.h +++ b/src/coordinates.h @@ -11,7 +11,9 @@ public: Coordinates() {_lon = NAN; _lat = NAN;} Coordinates(const Coordinates &c) {_lon = c._lon; _lat = c._lat;} Coordinates(qreal lon, qreal lat) {_lon = lon; _lat = lat;} + Coordinates(const QPointF &p) {_lon = p.x(), _lat = p.y();} + QPointF toPointF() const {return QPointF(_lon, _lat);} qreal &rlon() {return _lon;} qreal &rlat() {return _lat;} @@ -29,11 +31,6 @@ public: qreal distanceTo(const Coordinates &c) const; QPair boundingRect(qreal distance) const; - QPointF toMercator() const; - static Coordinates fromMercator(const QPointF &m); - - QPointF toPointF() const {return QPointF(_lon, _lat);} - private: qreal _lat, _lon; }; diff --git a/src/emptymap.cpp b/src/emptymap.cpp index 011c0c76..71355584 100644 --- a/src/emptymap.cpp +++ b/src/emptymap.cpp @@ -4,6 +4,7 @@ #include "rd.h" #include "wgs84.h" #include "coordinates.h" +#include "mercator.h" #include "emptymap.h" @@ -27,7 +28,7 @@ qreal EmptyMap::zoomFit(const QSize &size, const QRectF &br) else { Coordinates topLeft(br.topLeft()); Coordinates bottomRight(br.bottomRight()); - QRectF tbr(topLeft.toMercator(), bottomRight.toMercator()); + QRectF tbr(Mercator().ll2xy(topLeft), Mercator().ll2xy(bottomRight)); QPointF sc(tbr.width() / size.width(), tbr.height() / size.height()); @@ -65,12 +66,12 @@ void EmptyMap::draw(QPainter *painter, const QRectF &rect) QPointF EmptyMap::ll2xy(const Coordinates &c) const { - QPointF m = c.toMercator(); + QPointF m = Mercator().ll2xy(c); return QPointF(m.x() / _scale, m.y() / -_scale); } Coordinates EmptyMap::xy2ll(const QPointF &p) const { QPointF m(p.x() * _scale, -p.y() * _scale); - return Coordinates::fromMercator(m); + return Mercator().xy2ll(m); } diff --git a/src/mercator.cpp b/src/mercator.cpp new file mode 100644 index 00000000..1d0e3dfa --- /dev/null +++ b/src/mercator.cpp @@ -0,0 +1,13 @@ +#include +#include "rd.h" +#include "mercator.h" + +QPointF Mercator::ll2xy(const Coordinates &c) const +{ + return QPointF(c.lon(), rad2deg(log(tan(M_PI/4.0 + deg2rad(c.lat())/2.0)))); +} + +Coordinates Mercator::xy2ll(const QPointF &p) const +{ + return Coordinates(p.x(), rad2deg(2 * atan(exp(deg2rad(p.y()))) - M_PI/2)); +} diff --git a/src/mercator.h b/src/mercator.h new file mode 100644 index 00000000..0ec7ca96 --- /dev/null +++ b/src/mercator.h @@ -0,0 +1,13 @@ +#ifndef MERCATOR_H +#define MERCATOR_H + +#include "projection.h" + +class Mercator : public Projection +{ +public: + virtual QPointF ll2xy(const Coordinates &c) const; + virtual Coordinates xy2ll(const QPointF &p) const; +}; + +#endif // MERCATOR_H diff --git a/src/offlinemap.cpp b/src/offlinemap.cpp index b208733a..b14f3fa2 100644 --- a/src/offlinemap.cpp +++ b/src/offlinemap.cpp @@ -12,10 +12,13 @@ #include "wgs84.h" #include "coordinates.h" #include "matrix.h" +#include "mercator.h" +#include "transversemercator.h" #include "offlinemap.h" -int OfflineMap::parseMapFile(QIODevice &device, QList &points) +int OfflineMap::parseMapFile(QIODevice &device, QList &points, + QString &projection, double params[8]) { bool res; int ln = 1; @@ -71,6 +74,15 @@ int OfflineMap::parseMapFile(QIODevice &device, QList &points) if (!res) return ln; _size = QSize(w, h); + } else if (key == "Map Projection") { + projection = list.at(1); + } else if (key == "Projection Setup") { + if (list.count() < 9) + return ln; + for (int i = 1; i < 9; i++) { + params[i-1] = 0; + params[i-1] = list.at(i).trimmed().toFloat(&res); + } } } @@ -80,6 +92,22 @@ int OfflineMap::parseMapFile(QIODevice &device, QList &points) return 0; } +bool OfflineMap::createProjection(const QString &projection, double params[8]) +{ + if (projection == "Mercator") { + _projection = new Mercator(); + return true; + } else if (projection == "Transverse Mercator") { + _projection = new TransverseMercator(params[1], params[2], params[3], + params[4]); + return true; + } + + qWarning("%s: %s: unsupported map projection", qPrintable(_name), + qPrintable(projection)); + return false; +} + bool OfflineMap::computeTransformation(const QList &points) { if (points.count() < 2) { @@ -94,7 +122,7 @@ bool OfflineMap::computeTransformation(const QList &points) for (size_t k = 0; k < c.h(); k++) { for (int i = 0; i < points.size(); i++) { double f[3], t[2]; - QPointF p = points.at(i).second.toMercator(); + QPointF p = _projection->ll2xy(points.at(i).second); f[0] = p.x(); f[1] = p.y(); @@ -110,7 +138,7 @@ bool OfflineMap::computeTransformation(const QList &points) Q.zeroize(); for (int qi = 0; qi < points.size(); qi++) { double v[3]; - QPointF p = points.at(qi).second.toMercator(); + QPointF p = _projection->ll2xy(points.at(qi).second); v[0] = p.x(); v[1] = p.y(); @@ -242,8 +270,13 @@ OfflineMap::OfflineMap(const QString &path, QObject *parent) : Map(parent) { int errorLine = -2; QList points; + QString proj; + double params[8]; + _valid = false; + _img = 0; + _projection = 0; QFileInfo fi(path); _name = fi.fileName(); @@ -263,7 +296,7 @@ OfflineMap::OfflineMap(const QString &path, QObject *parent) : Map(parent) if (tarFiles.at(j).endsWith(".map")) { QByteArray ba = _tar.file(tarFiles.at(j)); QBuffer buffer(&ba); - errorLine = parseMapFile(buffer, points); + errorLine = parseMapFile(buffer, points, proj, params); _imgPath = QString(); break; } @@ -271,13 +304,15 @@ OfflineMap::OfflineMap(const QString &path, QObject *parent) : Map(parent) break; } else if (fileName.endsWith(".map")) { QFile mapFile(mapFiles.at(i).absoluteFilePath()); - errorLine = parseMapFile(mapFile, points); + errorLine = parseMapFile(mapFile, points, proj, params); break; } } if (!mapLoaded(errorLine)) return; + if (!createProjection(proj, params)) + return; if (!computeTransformation(points)) return; computeResolution(points); @@ -301,7 +336,6 @@ OfflineMap::OfflineMap(const QString &path, QObject *parent) : Map(parent) } } - _img = 0; _valid = true; } @@ -310,8 +344,13 @@ OfflineMap::OfflineMap(Tar &tar, const QString &path, QObject *parent) { int errorLine = -2; QList points; + QString proj; + double params[8]; + _valid = false; + _img = 0; + _projection = 0; QFileInfo fi(path); _name = fi.fileName(); @@ -323,12 +362,15 @@ OfflineMap::OfflineMap(Tar &tar, const QString &path, QObject *parent) if (tarFiles.at(j).startsWith(prefix)) { QByteArray ba = tar.file(tarFiles.at(j)); QBuffer buffer(&ba); - errorLine = parseMapFile(buffer, points); + errorLine = parseMapFile(buffer, points, proj, params); break; } } if (!mapLoaded(errorLine)) return; + + if (!createProjection(proj, params)) + return; if (!totalSizeSet()) return; if (!computeTransformation(points)) @@ -344,10 +386,17 @@ OfflineMap::OfflineMap(Tar &tar, const QString &path, QObject *parent) } _imgPath = QString(); - _img = 0; _valid = true; } +OfflineMap::~OfflineMap() +{ + if (_img) + delete _img; + if (_projection) + delete _projection; +} + void OfflineMap::load() { if (!_tarPath.isNull() && !_tileSize.isValid()) { @@ -454,10 +503,10 @@ void OfflineMap::draw(QPainter *painter, const QRectF &rect) QPointF OfflineMap::ll2xy(const Coordinates &c) const { - return _transform.map(c.toMercator()); + return _transform.map(_projection->ll2xy(c)); } Coordinates OfflineMap::xy2ll(const QPointF &p) const { - return Coordinates::fromMercator(_transform.inverted().map(p)); + return _projection->xy2ll(_transform.inverted().map(p)); } diff --git a/src/offlinemap.h b/src/offlinemap.h index c3524932..bbcd0919 100644 --- a/src/offlinemap.h +++ b/src/offlinemap.h @@ -8,6 +8,7 @@ class QIODevice; class QImage; +class Projection; class OfflineMap : public Map { @@ -16,6 +17,7 @@ class OfflineMap : public Map public: OfflineMap(const QString &path, QObject *parent = 0); OfflineMap(Tar &tar, const QString &path, QObject *parent = 0); + ~OfflineMap(); const QString &name() const {return _name;} @@ -40,9 +42,11 @@ public: private: typedef QPair ReferencePoint; - int parseMapFile(QIODevice &device, QList &points); + int parseMapFile(QIODevice &device, QList &points, + QString &projection, double params[8]); bool mapLoaded(int res); bool totalSizeSet(); + bool createProjection(const QString &projection, double params[8]); bool computeTransformation(const QList &points); bool computeResolution(QList &points); bool getTileInfo(const QStringList &tiles, const QString &path = QString()); @@ -50,6 +54,7 @@ private: QString _name; QSize _size; + Projection *_projection; QTransform _transform; qreal _resolution; diff --git a/src/onlinemap.cpp b/src/onlinemap.cpp index 2fc10d45..315eca1a 100644 --- a/src/onlinemap.cpp +++ b/src/onlinemap.cpp @@ -7,6 +7,7 @@ #include "wgs84.h" #include "misc.h" #include "coordinates.h" +#include "mercator.h" #include "onlinemap.h" @@ -173,7 +174,7 @@ qreal OnlineMap::zoomFit(const QSize &size, const QRectF &br) else { Coordinates topLeft(br.topLeft()); Coordinates bottomRight(br.bottomRight()); - QRectF tbr(topLeft.toMercator(), bottomRight.toMercator()); + QRectF tbr(Mercator().ll2xy(topLeft), Mercator().ll2xy(bottomRight)); QPointF sc(tbr.width() / size.width(), tbr.height() / size.height()); @@ -234,12 +235,12 @@ void OnlineMap::draw(QPainter *painter, const QRectF &rect) QPointF OnlineMap::ll2xy(const Coordinates &c) const { - QPointF m = c.toMercator(); + QPointF m = Mercator().ll2xy(c); return QPointF(m.x() / _scale, m.y() / -_scale); } Coordinates OnlineMap::xy2ll(const QPointF &p) const { QPointF m(p.x() * _scale, -p.y() * _scale); - return Coordinates::fromMercator(m); + return Mercator().xy2ll(m); } diff --git a/src/projection.h b/src/projection.h new file mode 100644 index 00000000..e487468a --- /dev/null +++ b/src/projection.h @@ -0,0 +1,15 @@ +#ifndef PROJECTION_H +#define PROJECTION_H + +#include +#include "coordinates.h" + +class Projection { +public: + virtual ~Projection() {} + + virtual QPointF ll2xy(const Coordinates &c) const = 0; + virtual Coordinates xy2ll(const QPointF &p) const = 0; +}; + +#endif // PROJECTION_H diff --git a/src/transversemercator.cpp b/src/transversemercator.cpp new file mode 100644 index 00000000..e4235d92 --- /dev/null +++ b/src/transversemercator.cpp @@ -0,0 +1,97 @@ +#include +#include "rd.h" +#include "wgs84.h" +#include "transversemercator.h" + +TransverseMercator::TransverseMercator(double centralMeridian, double scale, + double falseEasting, double falseNorthing) +{ + _centralMeridian = centralMeridian; + _scale = scale; + _falseEasting = falseEasting; + _falseNorthing = falseNorthing; +} + +QPointF TransverseMercator::ll2xy(const Coordinates &c) const +{ + QPointF p; + + const double e2 = WGS84_FLATTENING * (2 - WGS84_FLATTENING); + const double n = WGS84_FLATTENING / (2 - WGS84_FLATTENING); + const double rectifyingRadius = WGS84_RADIUS / (1 + n) + * (1 + 0.25*pow(n, 2) + 0.015625*pow(n, 4)); + + double A = e2; + double B = (5 * pow(e2, 2) - pow(e2, 3)) / 6.0; + double C = (104 * pow(e2, 3) - 45 * pow(e2, 4)) / 120.0; + double D = (1237 * pow(e2, 4)) / 1260.0; + + double phi = deg2rad(c.lat()); + double lambda = deg2rad(c.lon()); + double lambda0 = deg2rad(_centralMeridian); + + double deltaLambda = lambda - lambda0; + + double phiStar = phi - sin(phi) * cos(phi) * (A + B*pow(sin(phi), 2) + + C*pow(sin(phi), 4) + D*pow(sin(phi), 6)); + + double xiPrim = atan(tan(phiStar) / cos(deltaLambda)); + double etaPrim = atanh(cos(phiStar) * sin(deltaLambda)); + + double beta1 = 1/2.0 * n - 2/3.0 * pow(n, 2) + 5/16.0 * pow(n, 3) + + 41/180.0 * pow(n, 4); + double beta2 = 13/48.0 * pow(n, 2) - 3/5.0 * pow(n, 3) + 557/1440.0 + * pow(n, 4); + double beta3 = 61/240.0 * pow(n, 3) - 103/140.0 * pow(n, 4); + double beta4 = 49561/161280.0 * pow(n, 4); + + p.ry() = _falseNorthing + _scale * rectifyingRadius * (xiPrim + beta1 + * sin(2*xiPrim) * cosh(2*etaPrim) + beta2 * sin(4*xiPrim) + * cosh(4*etaPrim) + beta3 * sin(6*xiPrim) * cosh(6*etaPrim) + beta4 + * sin(8*xiPrim) * cosh(8*etaPrim)); + p.rx() = _falseEasting + _scale * rectifyingRadius * (etaPrim + beta1 + * cos(2*xiPrim) * sinh(2*etaPrim) + beta2 * cos(4*xiPrim) + * sinh(4*etaPrim) + beta3 * cos(6*xiPrim) * sinh(6*etaPrim) + beta4 + * cos(8*xiPrim) * sinh(8*etaPrim)); + + return p; +} + +Coordinates TransverseMercator::xy2ll(const QPointF &p) const +{ + const double e2 = WGS84_FLATTENING * (2 - WGS84_FLATTENING); + const double n = WGS84_FLATTENING / (2 - WGS84_FLATTENING); + const double rectifyingRadius = WGS84_RADIUS / (1 + n) + * (1 + 0.25*pow(n, 2) + 0.015625*pow(n, 4)); + + double xi = (p.y() - _falseNorthing) / (_scale * rectifyingRadius); + double eta = (p.x() - _falseEasting) / (_scale * rectifyingRadius); + + double delta1 = 1/2.0 * n - 2/3.0 * pow(n, 2) + 37/96.0 * pow(n, 3) + - 1/360.0 * pow(n, 4); + double delta2 = 1/48.0 * pow(n, 2) + 1/15.0 * pow(n, 3) - 437/1440.0 + * pow(n, 4); + double delta3 = 17/480.0 * pow(n, 3) - 37/840.0 * pow(n, 4); + double delta4 = 4397/161280.0 * pow(n, 4); + + double xiPrim = xi - delta1 * sin(2*xi) * cosh(2*eta) - delta2 * sin(4*xi) + * cosh(4*eta) - delta3 * sin(6*xi) * cosh(6*eta) - delta4 * sin(8*xi) + * cosh(8*eta); + double etaPrim = eta - delta1 * cos(2*xi) * sinh(2*eta) - delta2 * cos(4*xi) + * sinh(4*eta) - delta3 * cos(6*xi) * sinh(6*eta) - delta4 * cos(8*xi) + * sinh(8*eta); + + double phiStar = asin(sin(xiPrim) / cosh(etaPrim)); + double deltaLambda = atan(sinh(etaPrim) / cos(xiPrim)); + + double AStar = e2 + pow(e2, 2) + pow(e2, 3) + pow(e2, 4); + double BStar = (7 * pow(e2, 2) + 17 * pow(e2, 3) + 30 * pow(e2, 4)) / -6; + double CStar = (224 * pow(e2, 3) + 889 * pow(e2, 4)) / 120; + double DStar = (4279 * pow(e2, 4)) / -1260; + + double phi = phiStar + sin(phiStar) * cos(phiStar) * (AStar + BStar + * pow(sin(phiStar), 2) + CStar * pow(sin(phiStar), 4) + DStar + * pow(sin(phiStar), 6)); + + return Coordinates(_centralMeridian + rad2deg(deltaLambda), rad2deg(phi)); +} diff --git a/src/transversemercator.h b/src/transversemercator.h new file mode 100644 index 00000000..00fdf74d --- /dev/null +++ b/src/transversemercator.h @@ -0,0 +1,22 @@ +#ifndef TRANSVERSEMERCATOR_H +#define TRANSVERSEMERCATOR_H + +#include "projection.h" + +class TransverseMercator : public Projection +{ +public: + TransverseMercator(double centralMeridian, double scale, + double falseEasting, double falseNorthing); + + virtual QPointF ll2xy(const Coordinates &c) const; + virtual Coordinates xy2ll(const QPointF &p) const; + +private: + double _centralMeridian; + double _scale; + double _falseEasting; + double _falseNorthing; +}; + +#endif // TRANSVERSEMERCATOR_H diff --git a/src/wgs84.h b/src/wgs84.h index 766f3440..858fbeb1 100644 --- a/src/wgs84.h +++ b/src/wgs84.h @@ -1,6 +1,7 @@ #ifndef WGS84_H #define WGS84_H -#define WGS84_RADIUS 6378137.0 +#define WGS84_RADIUS 6378137.0 +#define WGS84_FLATTENING (1.0/298.257223563) #endif // WGS84_H