Added support for Transversal Mercator projection

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

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() << ", "

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<Coordinates, Coordinates> boundingRect(qreal distance) const;
 
-	QPointF toMercator() const;
-	static Coordinates fromMercator(const QPointF &m);
-
-	QPointF toPointF() const {return QPointF(_lon, _lat);}
-
 private:
 	qreal _lat, _lon;
 };

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);
 }

src/mercator.cpp

@@ -0,0 +1,13 @@
+#include <cmath>
+#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));
+}

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

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<ReferencePoint> &points)
+int OfflineMap::parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
+  QString &projection, double params[8])
 {
 	bool res;
 	int ln = 1;
@@ -71,6 +74,15 @@ int OfflineMap::parseMapFile(QIODevice &device, QList<ReferencePoint> &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<ReferencePoint> &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<ReferencePoint> &points)
 {
 	if (points.count() < 2) {
@@ -94,7 +122,7 @@ bool OfflineMap::computeTransformation(const QList<ReferencePoint> &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<ReferencePoint> &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<ReferencePoint> 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<ReferencePoint> 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));
 }

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<QPoint, Coordinates> ReferencePoint;
 
-	int parseMapFile(QIODevice &device, QList<ReferencePoint> &points);
+	int parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
+	  QString &projection, double params[8]);
 	bool mapLoaded(int res);
 	bool totalSizeSet();
+	bool createProjection(const QString &projection, double params[8]);
 	bool computeTransformation(const QList<ReferencePoint> &points);
 	bool computeResolution(QList<ReferencePoint> &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;
 

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);
 }

src/projection.h

@@ -0,0 +1,15 @@
+#ifndef PROJECTION_H
+#define PROJECTION_H
+
+#include <QPointF>
+#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

src/transversemercator.cpp

@@ -0,0 +1,97 @@
+#include <cmath>
+#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));
+}

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

src/wgs84.h

@@ -2,5 +2,6 @@
 #define WGS84_H
 
 #define WGS84_RADIUS     6378137.0
+#define WGS84_FLATTENING (1.0/298.257223563)
 
 #endif // WGS84_H