version++

README.md

@@ -3,7 +3,8 @@ GPXSee is a Qt-based GPS log file viewer and analyzer that supports GPX, TCX,
 KML, FIT, IGC and NMEA files.
 
 ## Features
-* User-definable map sources.
+* User-definable online maps.
+* Offline maps (image-based OziExplorer maps and tiled TrekBuddy maps/atlases).
 * Elevation, speed, heart rate, cadence, power and temperature graphs.
 * Support for multiple tracks in one view.
 * Support for POI files.

gpxsee.pro

@@ -1,5 +1,5 @@
 TARGET = GPXSee
-VERSION = 4.0
+VERSION = 4.1
 QT += core \
     gui \
     network
@@ -89,7 +89,8 @@ HEADERS += src/config.h \
     src/mercator.h \
     src/transversemercator.h \
     src/latlon.h \
-    src/utm.h
+    src/utm.h \
+    src/lambertconic.h
 SOURCES += src/main.cpp \
     src/gui.cpp \
     src/poi.cpp \
@@ -152,7 +153,8 @@ SOURCES += src/main.cpp \
     src/atlas.cpp \
     src/mercator.cpp \
     src/transversemercator.cpp \
-    src/utm.cpp
+    src/utm.cpp \
+    src/lambertconic.cpp
 RESOURCES += gpxsee.qrc
 TRANSLATIONS = lang/gpxsee_cs.ts \
     lang/gpxsee_sv.ts

pkg/gpxsee.nsi

@@ -5,7 +5,7 @@
 ; The name of the installer
 Name "GPXSee"
 ; Program version
-!define VERSION "4.0"
+!define VERSION "4.1"
 
 ; The file to write
 OutFile "GPXSee-${VERSION}.exe"

pkg/gpxsee64.nsi

@@ -5,7 +5,7 @@
 ; The name of the installer
 Name "GPXSee"
 ; Program version
-!define VERSION "4.0"
+!define VERSION "4.1"
 
 ; The file to write
 OutFile "GPXSee-${VERSION}_x64.exe"

src/atlas.cpp

@@ -107,6 +107,7 @@ void Atlas::computeBounds()
 Atlas::Atlas(const QString &path, QObject *parent) : Map(parent)
 {
 	_valid = false;
+	_zoom = 0;
 
 	QFileInfo fi(path);
 	_name = fi.fileName();
@@ -139,7 +140,6 @@ Atlas::Atlas(const QString &path, QObject *parent) : Map(parent)
 	computeZooms();
 	computeBounds();
 
-	_zoom = 0;
 	_valid = true;
 }
 

src/lambertconic.cpp

@@ -0,0 +1,83 @@
+#include <cmath>
+#include "rd.h"
+#include "wgs84.h"
+#include "lambertconic.h"
+
+
+#ifndef M_PI_2
+	#define M_PI_2 1.57079632679489661923
+#endif // M_PI_2
+#ifndef M_PI_4
+	#define M_PI_4 0.78539816339744830962
+#endif // M_PI_4
+
+static double q(double b)
+{
+	double e = sqrt(WGS84_FLATTENING * (2. - WGS84_FLATTENING));
+	double esb = e * sin(b);
+	return log(tan(M_PI_4 + b / 2.) * pow((1. - esb) / (1. + esb), e / 2.));
+}
+
+static double inv_q(double q)
+{
+	double e = sqrt(WGS84_FLATTENING * (2. - WGS84_FLATTENING));
+	double b0 = 0.;
+	double b = 2. * atan(exp(q)) - M_PI_2;
+
+	do {
+		b0 = b;
+		double esb = e * sin(b);
+		b = 2. * atan(exp(q) * pow((1. - esb) / (1. + esb), -e / 2.)) - M_PI_2;
+	} while (fabs(b - b0) > 1e-10);
+
+	return b;
+}
+
+static double nradius(double phi)
+{
+	double sin_phi = sin(phi);
+	return (WGS84_RADIUS / sqrt(1. - (WGS84_FLATTENING
+	  * (2. - WGS84_FLATTENING)) * sin_phi * sin_phi));
+}
+
+LambertConic::LambertConic(double standardParallel1, double standardParallel2,
+  double centralParallel, double centralMeridian, double scale,
+  double falseEasting, double falseNorthing)
+{
+	_cm = centralMeridian;
+	_fe = falseEasting;
+	_fn = falseNorthing;
+
+	double sp1 = deg2rad(standardParallel1);
+	double sp2 = deg2rad(standardParallel2);
+
+	double N1 = nradius(sp1);
+	double N2 = nradius(sp2);
+
+	_q0 = q(deg2rad(centralParallel));
+	double q1 = q(sp1);
+	double q2 = q(sp2);
+
+	_n = log((N1 * cos(sp1)) / (N2 * cos(sp2))) / (q2 - q1);
+	double R1 = N1 * cos(sp1) / _n;
+	_R0 = scale * R1 * exp(_n * (q1 - _q0));
+}
+
+QPointF LambertConic::ll2xy(const Coordinates &c) const
+{
+	double dl = _n * (deg2rad(c.lon()) - deg2rad(_cm));
+	double R = _R0 * exp(_n * (_q0 - q(deg2rad(c.lat()))));
+
+	return QPointF(_fe + R * sin(dl), _fn + _R0 - R * cos(dl));
+}
+
+Coordinates LambertConic::xy2ll(const QPointF &p) const
+{
+	double dl = atan((p.x() - _fe) / (_fn + _R0 - p.y()));
+	double dx = p.x() - _fe;
+	double dy = p.y() - _fn - _R0;
+	double R = sqrt(dx * dx + dy * dy);
+	double q = _q0 - log(R / _R0) / _n;
+
+	return Coordinates(rad2deg(deg2rad(_cm) + dl / _n), rad2deg(inv_q(q)));
+}

src/lambertconic.h

@@ -0,0 +1,26 @@
+#ifndef LAMBERTCONIC_H
+#define LAMBERTCONIC_H
+
+#include "projection.h"
+
+class LambertConic : public Projection
+{
+public:
+	LambertConic(double standardParallel1, double standardParallel2,
+	  double centralParallel, double centralMeridian, double scale,
+	  double falseEasting, double falseNorthing);
+
+	virtual QPointF ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const QPointF &p) const;
+
+private:
+	double _cm;
+	double _fe;
+	double _fn;
+
+	double _q0;
+	double _R0;
+	double _n;
+};
+
+#endif // LAMBERTCONIC_H

src/offlinemap.cpp

@@ -16,6 +16,7 @@
 #include "mercator.h"
 #include "transversemercator.h"
 #include "utm.h"
+#include "lambertconic.h"
 #include "offlinemap.h"
 
 
@@ -104,18 +105,27 @@ int OfflineMap::parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
 			} else if (key == "Projection Setup") {
 				if (list.count() < 8)
 					return ln;
+				setup.centralParallel = list.at(1).trimmed().toFloat(&res);
+				if (!res)
+					setup.centralParallel = 0;
 				setup.centralMeridian = list.at(2).trimmed().toFloat(&res);
 				if (!res)
 					setup.centralMeridian = 0;
 				setup.scale = list.at(3).trimmed().toFloat(&res);
 				if (!res)
-					setup.scale = 0;
+					setup.scale = 1.0;
 				setup.falseEasting = list.at(4).trimmed().toFloat(&res);
 				if (!res)
 					setup.falseEasting = 0;
 				setup.falseNorthing = list.at(5).trimmed().toFloat(&res);
 				if (!res)
 					setup.falseNorthing = 0;
+				setup.standardParallel1 = list.at(6).trimmed().toFloat(&res);
+				if (!res)
+					setup.standardParallel1 = 0;
+				setup.standardParallel2 = list.at(7).trimmed().toFloat(&res);
+				if (!res)
+					setup.standardParallel2 = 0;
 			}
 		}
 
@@ -141,6 +151,10 @@ bool OfflineMap::createProjection(const QString &projection,
 		  setup.falseEasting, setup.falseNorthing);
 	else if (projection == "Latitude/Longitude")
 		_projection = new LatLon();
+	else if (projection == "Lambert Conformal Conic")
+		_projection = new LambertConic(setup.standardParallel1,
+		  setup.standardParallel2, setup.centralParallel, setup.centralMeridian,
+		  setup.scale, setup.falseEasting, setup.falseNorthing);
 	else if (projection == "(UTM) Universal Transverse Mercator") {
 		if (setup.zone)
 			_projection = new UTM(setup.zone);
@@ -328,6 +342,7 @@ OfflineMap::OfflineMap(const QString &path, QObject *parent) : Map(parent)
 	_valid = false;
 	_img = 0;
 	_projection = 0;
+	_resolution = 0;
 
 	QFileInfo fi(path);
 	_name = fi.fileName();

src/offlinemap.h

@@ -51,10 +51,13 @@ private:
 	} ReferencePoint;
 
 	typedef struct {
+		double centralParallel;
 		double centralMeridian;
 		double scale;
 		double falseEasting;
 		double falseNorthing;
+		double standardParallel1;
+		double standardParallel2;
 		int zone;
 	} ProjectionSetup;
 

src/onlinemap.cpp

@@ -45,7 +45,6 @@ OnlineMap::OnlineMap(const QString &name, const QString &url, QObject *parent)
 {
 	_name = name;
 	_url = url;
-	downloader = downloader;
 	_block = false;
 	_scale = ((360.0/(qreal)(1<<ZOOM_MAX))/(qreal)TILE_SIZE);
 

src/transversemercator.cpp

@@ -4,14 +4,6 @@
 #include "transversemercator.h"
 
 
-TransverseMercator::TransverseMercator()
-{
-	_centralMeridian = 0;
-	_scale = 1.0;
-	_falseEasting = 0;
-	_falseNorthing = 0;
-}
-
 TransverseMercator::TransverseMercator(double centralMeridian, double scale,
   double falseEasting, double falseNorthing)
 {
@@ -19,48 +11,58 @@ TransverseMercator::TransverseMercator(double centralMeridian, double scale,
 	_scale = scale;
 	_falseEasting = falseEasting;
 	_falseNorthing = falseNorthing;
+
+	const double e2 = WGS84_FLATTENING * (2 - WGS84_FLATTENING);
+	const double n = WGS84_FLATTENING / (2 - WGS84_FLATTENING);
+	_rectifyingRadius = WGS84_RADIUS / (1 + n)
+	  * (1 + 0.25*pow(n, 2) + 0.015625*pow(n, 4));
+
+	_A = e2;
+	_B = (5 * pow(e2, 2) - pow(e2, 3)) / 6.0;
+	_C = (104 * pow(e2, 3) - 45 * pow(e2, 4)) / 120.0;
+	_D = (1237 * pow(e2, 4)) / 1260.0;
+
+	_beta1 = 1/2.0 * n - 2/3.0 * pow(n, 2) + 5/16.0 * pow(n, 3) + 41/180.0
+	  * pow(n, 4);
+	_beta2 = 13/48.0 * pow(n, 2) - 3/5.0 * pow(n, 3) + 557/1440.0 * pow(n, 4);
+	_beta3 = 61/240.0 * pow(n, 3) - 103/140.0 * pow(n, 4);
+	_beta4 = 49561/161280.0 * pow(n, 4);
+
+	_delta1 = 1/2.0 * n - 2/3.0 * pow(n, 2) + 37/96.0 * pow(n, 3) - 1/360.0
+	  * pow(n, 4);
+	_delta2 = 1/48.0 * pow(n, 2) + 1/15.0 * pow(n, 3) - 437/1440.0 * pow(n, 4);
+	_delta3 = 17/480.0 * pow(n, 3) - 37/840.0 * pow(n, 4);
+	_delta4 = 4397/161280.0 * pow(n, 4);
+
+	_AStar = e2 + pow(e2, 2) + pow(e2, 3) + pow(e2, 4);
+	_BStar = (7 * pow(e2, 2) + 17 * pow(e2, 3) + 30 * pow(e2, 4)) / -6;
+	_CStar = (224 * pow(e2, 3) + 889 * pow(e2, 4)) / 120;
+	_DStar = (4279 * pow(e2, 4)) / -1260;
 }
 
 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 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
+	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
+	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;
@@ -68,38 +70,21 @@ QPointF TransverseMercator::ll2xy(const Coordinates &c) const
 
 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 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)
+	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 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
+	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

@@ -6,7 +6,6 @@
 class TransverseMercator : public Projection
 {
 public:
-	TransverseMercator();
 	TransverseMercator(double centralMeridian, double scale,
 	  double falseEasting, double falseNorthing);
 
@@ -18,6 +17,12 @@ private:
 	double _scale;
 	double _falseEasting;
 	double _falseNorthing;
+
+	double _rectifyingRadius;
+	double _A, _B, _C, _D;
+	double _beta1, _beta2, _beta3, _beta4;
+	double _delta1, _delta2, _delta3, _delta4;
+	double _AStar, _BStar, _CStar, _DStar;
 };
 
 #endif // TRANSVERSEMERCATOR_H

src/utm.cpp

@@ -1,6 +1,6 @@
 #include "utm.h"
 
-UTM::UTM(const Coordinates &c)
+UTM::UTM(const Coordinates &c) : _tm(0, 1.0, 0, 0)
 {
 	int zone = int((c.lon() + 180)/6) + 1;