Added support for Lambert Conformal Conic projection

gpxsee.pro

@@ -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

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

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;