Use ellipsoids defined by the map file for the projection computation

gpxsee.pro

@@ -90,7 +90,8 @@ HEADERS += src/config.h \
     src/transversemercator.h \
     src/latlon.h \
     src/utm.h \
-    src/lambertconic.h
+    src/lambertconic.h \
+    src/ellipsoid.h
 SOURCES += src/main.cpp \
     src/gui.cpp \
     src/poi.cpp \
@@ -154,7 +155,8 @@ SOURCES += src/main.cpp \
     src/mercator.cpp \
     src/transversemercator.cpp \
     src/utm.cpp \
-    src/lambertconic.cpp
+    src/lambertconic.cpp \
+    src/ellipsoid.cpp
 RESOURCES += gpxsee.qrc
 TRANSLATIONS = lang/gpxsee_cs.ts \
     lang/gpxsee_sv.ts

src/ellipsoid.cpp

@@ -0,0 +1,92 @@
+#include <cmath>
+#include <QString>
+#include "wgs84.h"
+#include "ellipsoid.h"
+
+
+#define INTERNATIONAL_RADIUS     6378388.0
+#define INTERNATIONAL_FLATTENING (1.0/297.0)
+#define KRASSOVSKY_RADIUS        6378245.0
+#define KRASSOVSKY_FLATTENING    (1.0/298.3)
+#define BESSEL_RADIUS            6377397.155
+#define BESSEL_FLATTENING        (1.0/299.1528128)
+#define GRS80_RADIUS             6378137.0
+#define GRS80_FLATTENING         (1.0/298.257222101)
+#define WGS70_RADIUS             6378135.0
+#define WGS70_FLATTENING         (1.0/298.26)
+
+#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
+
+#define ARRAY_SIZE(array) (sizeof(array) / sizeof((array)[0]))
+
+typedef struct {
+	const char *name;
+	double radius;
+	double flattening;
+} Entry;
+
+static Entry list[] = {
+	{"S42", KRASSOVSKY_RADIUS, KRASSOVSKY_FLATTENING},
+	{"Pulkovo 1942", KRASSOVSKY_RADIUS, KRASSOVSKY_FLATTENING},
+	{"European 1950", INTERNATIONAL_RADIUS, INTERNATIONAL_FLATTENING},
+	{"European 1979", INTERNATIONAL_RADIUS, INTERNATIONAL_FLATTENING},
+	{"NZGD1949", INTERNATIONAL_RADIUS, INTERNATIONAL_FLATTENING},
+	{"NAD27", BESSEL_RADIUS, BESSEL_FLATTENING},
+	{"NAD83", GRS80_RADIUS, GRS80_FLATTENING},
+	{"WGS 72", WGS70_RADIUS, WGS70_FLATTENING}
+};
+
+
+Ellipsoid::Ellipsoid()
+{
+	_radius = WGS84_RADIUS;
+	_flattening = WGS84_FLATTENING;
+}
+
+Ellipsoid::Ellipsoid(const QString &datum)
+{
+	for (size_t i = 0; i < ARRAY_SIZE(list); i++) {
+		if (datum.startsWith(list[i].name)) {
+			_radius = list[i].radius;
+			_flattening = list[i].flattening;
+			return;
+		}
+	}
+
+	_radius = WGS84_RADIUS;
+	_flattening = WGS84_FLATTENING;
+}
+
+double Ellipsoid::q(double b) const
+{
+	double e = sqrt(_flattening * (2. - _flattening));
+	double esb = e * sin(b);
+	return log(tan(M_PI_4 + b / 2.) * pow((1. - esb) / (1. + esb), e / 2.));
+}
+
+double Ellipsoid::iq(double q) const
+{
+	double e = sqrt(_flattening * (2. - _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;
+}
+
+double Ellipsoid::nradius(double phi) const
+{
+	double sin_phi = sin(phi);
+	return (_radius / sqrt(1. - (_flattening * (2. - _flattening)) * sin_phi
+	  * sin_phi));
+}

src/ellipsoid.h

@@ -0,0 +1,24 @@
+#ifndef ELLIPSOID_H
+#define ELLIPSOID_H
+
+class QString;
+
+class Ellipsoid
+{
+public:
+	Ellipsoid();
+	Ellipsoid(const QString &datum);
+
+	double radius() const {return _radius;}
+	double flattening() const {return _flattening;}
+
+	double q(double b) const;
+	double iq(double q) const;
+	double nradius(double phi) const;
+
+private:
+	double _radius;
+	double _flattening;
+};
+
+#endif // ELLIPSOID_H

src/lambertconic.cpp

@@ -1,48 +1,11 @@
 #include <cmath>
 #include "rd.h"
-#include "wgs84.h"
 #include "lambertconic.h"
 
 
-#ifndef M_PI_2
+LambertConic::LambertConic(const Ellipsoid &ellipsoid, double standardParallel1,
-	#define M_PI_2 1.57079632679489661923
+  double standardParallel2, double centralParallel, double centralMeridian,
-#endif // M_PI_2
+  double scale, double falseEasting, double falseNorthing) : _e(ellipsoid)
-#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;
@@ -51,12 +14,12 @@ LambertConic::LambertConic(double standardParallel1, double standardParallel2,
 	double sp1 = deg2rad(standardParallel1);
 	double sp2 = deg2rad(standardParallel2);
 
-	double N1 = nradius(sp1);
+	double N1 = _e.nradius(sp1);
-	double N2 = nradius(sp2);
+	double N2 = _e.nradius(sp2);
 
-	_q0 = q(deg2rad(centralParallel));
+	_q0 = ellipsoid.q(deg2rad(centralParallel));
-	double q1 = q(sp1);
+	double q1 = _e.q(sp1);
-	double q2 = q(sp2);
+	double q2 = _e.q(sp2);
 
 	_n = log((N1 * cos(sp1)) / (N2 * cos(sp2))) / (q2 - q1);
 	double R1 = N1 * cos(sp1) / _n;
@@ -66,7 +29,7 @@ LambertConic::LambertConic(double standardParallel1, double standardParallel2,
 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()))));
+	double R = _R0 * exp(_n * (_q0 - _e.q(deg2rad(c.lat()))));
 
 	return QPointF(_fe + R * sin(dl), _fn + _R0 - R * cos(dl));
 }
@@ -79,5 +42,5 @@ Coordinates LambertConic::xy2ll(const QPointF &p) const
 	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)));
+	return Coordinates(rad2deg(deg2rad(_cm) + dl / _n), rad2deg(_e.iq(q)));
 }

src/lambertconic.h

@@ -1,19 +1,22 @@
 #ifndef LAMBERTCONIC_H
 #define LAMBERTCONIC_H
 
+#include "ellipsoid.h"
 #include "projection.h"
 
 class LambertConic : public Projection
 {
 public:
-	LambertConic(double standardParallel1, double standardParallel2,
+	LambertConic(const Ellipsoid &ellipsoid, double standardParallel1,
-	  double centralParallel, double centralMeridian, double scale,
+	  double standardParallel2, double centralParallel, double centralMeridian,
-	  double falseEasting, double falseNorthing);
+	  double scale, double falseEasting, double falseNorthing);
 
 	virtual QPointF ll2xy(const Coordinates &c) const;
 	virtual Coordinates xy2ll(const QPointF &p) const;
 
 private:
+	Ellipsoid _e;
+
 	double _cm;
 	double _fe;
 	double _fn;

src/offlinemap.cpp

@@ -12,6 +12,7 @@
 #include "wgs84.h"
 #include "coordinates.h"
 #include "matrix.h"
+#include "ellipsoid.h"
 #include "latlon.h"
 #include "mercator.h"
 #include "transversemercator.h"
@@ -21,7 +22,7 @@
 
 
 int OfflineMap::parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
-  QString &projection, ProjectionSetup &setup)
+  QString &projection, ProjectionSetup &setup, QString &datum)
 {
 	bool res;
 	int ln = 1;
@@ -38,6 +39,8 @@ int OfflineMap::parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
 				return ln;
 		} else if (ln == 3)
 			_imgPath = line.trimmed();
+		else if (ln == 5)
+			datum = line.split(',').at(0).trimmed();
 		else {
 			QList<QByteArray> list = line.split(',');
 			QString key(list.at(0).trimmed());
@@ -135,8 +138,9 @@ int OfflineMap::parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
 	return 0;
 }
 
-bool OfflineMap::createProjection(const QString &projection,
+bool OfflineMap::createProjection(const QString &datum,
-  const ProjectionSetup &setup, QList<ReferencePoint> &points)
+  const QString &projection, const ProjectionSetup &setup,
+  QList<ReferencePoint> &points)
 {
 	if (points.count() < 2) {
 		qWarning("%s: insufficient number of reference points",
@@ -147,12 +151,13 @@ bool OfflineMap::createProjection(const QString &projection,
 	if (projection == "Mercator")
 		_projection = new Mercator();
 	else if (projection == "Transverse Mercator")
-		_projection = new TransverseMercator(setup.centralMeridian, setup.scale,
+		_projection = new TransverseMercator(Ellipsoid(datum),
-		  setup.falseEasting, setup.falseNorthing);
+		  setup.centralMeridian, setup.scale, setup.falseEasting,
+		  setup.falseNorthing);
 	else if (projection == "Latitude/Longitude")
 		_projection = new LatLon();
 	else if (projection == "Lambert Conformal Conic")
-		_projection = new LambertConic(setup.standardParallel1,
+		_projection = new LambertConic(Ellipsoid(datum), setup.standardParallel1,
 		  setup.standardParallel2, setup.centralParallel, setup.centralMeridian,
 		  setup.scale, setup.falseEasting, setup.falseNorthing);
 	else if (projection == "(UTM) Universal Transverse Mercator") {
@@ -335,7 +340,7 @@ OfflineMap::OfflineMap(const QString &path, QObject *parent) : Map(parent)
 {
 	int errorLine = -2;
 	QList<ReferencePoint> points;
-	QString proj;
+	QString proj, datum;
 	ProjectionSetup setup;
 
 
@@ -362,7 +367,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, proj, setup);
+					errorLine = parseMapFile(buffer, points, proj, setup, datum);
 					_imgPath = QString();
 					break;
 				}
@@ -370,14 +375,14 @@ 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, proj, setup);
+			errorLine = parseMapFile(mapFile, points, proj, setup, datum);
 			break;
 		}
 	}
 	if (!mapLoaded(errorLine))
 		return;
 
-	if (!createProjection(proj, setup, points))
+	if (!createProjection(datum, proj, setup, points))
 		return;
 	if (!computeTransformation(points))
 		return;
@@ -410,7 +415,7 @@ OfflineMap::OfflineMap(Tar &tar, const QString &path, QObject *parent)
 {
 	int errorLine = -2;
 	QList<ReferencePoint> points;
-	QString proj;
+	QString proj, datum;
 	ProjectionSetup setup;
 
 
@@ -428,14 +433,14 @@ 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, proj, setup);
+			errorLine = parseMapFile(buffer, points, proj, setup, datum);
 			break;
 		}
 	}
 	if (!mapLoaded(errorLine))
 		return;
 
-	if (!createProjection(proj, setup, points))
+	if (!createProjection(datum, proj, setup, points))
 		return;
 	if (!totalSizeSet())
 		return;

src/offlinemap.h

@@ -62,10 +62,10 @@ private:
 	} ProjectionSetup;
 
 	int parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
-	  QString &projection, ProjectionSetup &setup);
+	  QString &projection, ProjectionSetup &setup, QString &datum);
 	bool mapLoaded(int res);
 	bool totalSizeSet();
-	bool createProjection(const QString &projection,
+	bool createProjection(const QString &datum, const QString &projection,
 	  const ProjectionSetup &setup, QList<ReferencePoint> &points);
 	bool computeTransformation(const QList<ReferencePoint> &points);
 	bool computeResolution(QList<ReferencePoint> &points);

src/transversemercator.cpp

@@ -1,20 +1,22 @@
 #include <cmath>
 #include "rd.h"
-#include "wgs84.h"
+#include "ellipsoid.h"
 #include "transversemercator.h"
 
 
-TransverseMercator::TransverseMercator(double centralMeridian, double scale,
+TransverseMercator::TransverseMercator(const Ellipsoid &ellipsoid,
-  double falseEasting, double falseNorthing)
+  double centralMeridian, double scale, double falseEasting,
+  double falseNorthing)
 {
 	_centralMeridian = centralMeridian;
 	_scale = scale;
 	_falseEasting = falseEasting;
 	_falseNorthing = falseNorthing;
 
-	const double e2 = WGS84_FLATTENING * (2 - WGS84_FLATTENING);
+
-	const double n = WGS84_FLATTENING / (2 - WGS84_FLATTENING);
+	const double e2 = ellipsoid.flattening() * (2 - ellipsoid.flattening());
-	_rectifyingRadius = WGS84_RADIUS / (1 + n)
+	const double n = ellipsoid.flattening() / (2 - ellipsoid.flattening());
+	_rectifyingRadius = ellipsoid.radius() / (1 + n)
 	  * (1 + 0.25*pow(n, 2) + 0.015625*pow(n, 4));
 
 	_A = e2;

src/transversemercator.h

@@ -3,11 +3,13 @@
 
 #include "projection.h"
 
+class Ellipsoid;
+
 class TransverseMercator : public Projection
 {
 public:
-	TransverseMercator(double centralMeridian, double scale,
+	TransverseMercator(const Ellipsoid &ellipsoid, double centralMeridian,
-	  double falseEasting, double falseNorthing);
+	  double scale, double falseEasting, double falseNorthing);
 
 	virtual QPointF ll2xy(const Coordinates &c) const;
 	virtual Coordinates xy2ll(const QPointF &p) const;

src/utm.cpp

@@ -1,6 +1,12 @@
+#include "ellipsoid.h"
 #include "utm.h"
 
-UTM::UTM(const Coordinates &c) : _tm(0, 1.0, 0, 0)
+UTM::UTM(int zone) : _tm(Ellipsoid(), (qAbs(zone) - 1)*6 - 180 + 3, 0.9996,
+  500000, zone < 0 ? 10000000 : 0)
+{
+}
+
+UTM::UTM(const Coordinates &c) : _tm(Ellipsoid(), 0, 0, 0, 0)
 {
 	int zone = int((c.lon() + 180)/6) + 1;
 
@@ -18,5 +24,6 @@ UTM::UTM(const Coordinates &c) : _tm(0, 1.0, 0, 0)
 	}
 	double cm = (zone - 1)*6 - 180 + 3;
 
-	_tm = TransverseMercator(cm, 0.9996, 500000, (c.lat() < 0) ? 10000000 : 0);
+	_tm = TransverseMercator(Ellipsoid(), cm, 0.9996, 500000,
+	  (c.lat() < 0) ? 10000000 : 0);
 }

src/utm.h

@@ -7,8 +7,7 @@
 class UTM : public Projection
 {
 public:
-	UTM(int zone) : _tm((qAbs(zone) - 1)*6 - 180 + 3, 0.9996, 500000,
+	UTM(int zone);
-	  zone < 0 ? 10000000 : 0) {}
 	UTM(const Coordinates &c);
 
 	virtual QPointF ll2xy(const Coordinates &c) const