Added support for Lambert Azimuthal Equal Area projections

gpxsee.pro

@@ -103,7 +103,8 @@ HEADERS += src/config.h \
     src/heartrategraphitem.h \
     src/temperaturegraphitem.h \
     src/cadencegraphitem.h \
-    src/powergraphitem.h
+    src/powergraphitem.h \
+    src/azimuthalequalarea.h
 SOURCES += src/main.cpp \
     src/gui.cpp \
     src/poi.cpp \
@@ -179,7 +180,8 @@ SOURCES += src/main.cpp \
     src/heartrategraphitem.cpp \
     src/temperaturegraphitem.cpp \
     src/cadencegraphitem.cpp \
-    src/powergraphitem.cpp
+    src/powergraphitem.cpp \
+    src/azimuthalequalarea.cpp
 RESOURCES += gpxsee.qrc
 TRANSLATIONS = lang/gpxsee_cs.ts \
     lang/gpxsee_sv.ts \

src/azimuthalequalarea.cpp

@@ -0,0 +1,194 @@
+/*
+ * Based on libgeotrans with the following Source Code Disclaimer:
+
+1. The GEOTRANS source code ("the software") is provided free of charge by
+the National Imagery and Mapping Agency (NIMA) of the United States
+Department of Defense. Although NIMA makes no copyright claim under Title 17
+U.S.C., NIMA claims copyrights in the source code under other legal regimes.
+NIMA hereby grants to each user of the software a license to use and
+distribute the software, and develop derivative works.
+
+2. Warranty Disclaimer: The software was developed to meet only the internal
+requirements of the U.S. National Imagery and Mapping Agency. The software
+is provided "as is," and no warranty, express or implied, including but not
+limited to the implied warranties of merchantability and fitness for
+particular purpose or arising by statute or otherwise in law or from a
+course of dealing or usage in trade, is made by NIMA as to the accuracy and
+functioning of the software.
+
+3. NIMA and its personnel are not required to provide technical support or
+general assistance with respect to the software.
+
+4. Neither NIMA nor its personnel will be liable for any claims, losses, or
+damages arising from or connected with the use of the software. The user
+agrees to hold harmless the United States National Imagery and Mapping
+Agency. The user's sole and exclusive remedy is to stop using the software.
+
+5. NIMA requests that products developed using the software credit the
+source of the software with the following statement, "The product was
+developed using GEOTRANS, a product of the National Imagery and Mapping
+Agency and U.S. Army Engineering Research and Development Center."
+
+6. For any products developed using the software, NIMA requires a disclaimer
+that use of the software does not indicate endorsement or approval of the
+product by the Secretary of Defense or the National Imagery and Mapping
+Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
+the National Imagery and Mapping Agency, the initials "NIMA", the seal of
+the National Imagery and Mapping Agency, or any colorable imitation thereof
+shall not be used to imply approval, endorsement, or authorization of a
+product without prior written permission from United States Secretary of
+Defense.
+
+*/
+
+#include "rd.h"
+#include "azimuthalequalarea.h"
+
+
+#ifndef M_PI_2
+	#define M_PI_2 1.57079632679489661923
+#endif // M_PI_2
+
+AzimuthalEqualArea::AzimuthalEqualArea(const Ellipsoid &ellipsoid,
+  double latitudeOrigin, double longitudeOrigin, double falseEasting,
+  double falseNorthing)
+{
+	double es2, es4, es6;
+
+	_e = ellipsoid;
+	es2 = 2 * _e.flattening() - _e.flattening() * _e.flattening();
+
+	es4 = es2 * es2;
+	es6 = es4 * es2;
+	_ra = _e.radius() * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6
+	  / 3024.0);
+	_latOrigin = deg2rad(latitudeOrigin);
+	_sinLatOrigin = sin(_latOrigin);
+	_cosLatOrigin = cos(_latOrigin);
+	_absLatOrigin = fabs(_latOrigin);
+	_lonOrigin = deg2rad(longitudeOrigin);
+	if (_lonOrigin > M_PI)
+		_lonOrigin -= 2.0 * M_PI;
+
+	_falseNorthing = falseNorthing;
+	_falseEasting = falseEasting;
+}
+
+QPointF AzimuthalEqualArea::ll2xy(const Coordinates &c) const
+{
+	double dlam;
+	double k_prime;
+	double cd;
+	double rlat = deg2rad(c.lat());
+	double slat = sin(rlat);
+	double clat = cos(rlat);
+	double cos_c;
+	double sin_dlam, cos_dlam;
+	double Ra_kprime;
+	double Ra_PI_OVER_2_Lat;
+	QPointF p;
+
+
+	dlam = deg2rad(c.lon()) - _lonOrigin;
+	if (dlam > M_PI)
+		dlam -= 2.0 * M_PI;
+	if (dlam < -M_PI)
+		dlam += 2.0 * M_PI;
+
+	sin_dlam = sin(dlam);
+	cos_dlam = cos(dlam);
+	if (fabs(_absLatOrigin - M_PI_2) < 1.0e-10) {
+		if (_latOrigin >= 0.0) {
+			Ra_PI_OVER_2_Lat = _ra * (M_PI_2 - rlat);
+			p.rx() = Ra_PI_OVER_2_Lat * sin_dlam + _falseEasting;
+			p.ry() = -1.0 * (Ra_PI_OVER_2_Lat * cos_dlam) + _falseNorthing;
+		} else {
+			Ra_PI_OVER_2_Lat = _ra * (M_PI_2 + rlat);
+			p.rx() = Ra_PI_OVER_2_Lat * sin_dlam + _falseEasting;
+			p.ry() = Ra_PI_OVER_2_Lat * cos_dlam + _falseNorthing;
+		}
+	} else if (_absLatOrigin <= 1.0e-10) {
+		cos_c = clat * cos_dlam;
+		if (fabs(fabs(cos_c) - 1.0) < 1.0e-14) {
+			if (cos_c >= 0.0) {
+				p.rx() = _falseEasting;
+				p.ry() = _falseNorthing;
+			} else
+				return QPointF(NAN, NAN);
+		} else {
+			cd = acos(cos_c);
+			k_prime = cd / sin(cd);
+			Ra_kprime = _ra * k_prime;
+			p.rx() = Ra_kprime * clat * sin_dlam + _falseEasting;
+			p.ry() = Ra_kprime * slat + _falseNorthing;
+		}
+	} else {
+		cos_c = (_sinLatOrigin * slat) + (_cosLatOrigin * clat * cos_dlam);
+		if (fabs(fabs(cos_c) - 1.0) < 1.0e-14) {
+			if (cos_c >= 0.0) {
+				p.rx() = _falseEasting;
+				p.ry() = _falseNorthing;
+			} else
+				return QPointF(NAN, NAN);
+		} else {
+			cd = acos(cos_c);
+			k_prime = cd / sin(cd);
+			Ra_kprime = _ra * k_prime;
+			p.rx() = Ra_kprime * clat * sin_dlam + _falseEasting;
+			p.ry() = Ra_kprime * (_cosLatOrigin * slat - _sinLatOrigin * clat
+			  * cos_dlam) + _falseNorthing;
+		}
+	}
+
+	return p;
+}
+
+Coordinates AzimuthalEqualArea::xy2ll(const QPointF &p) const
+{
+	double dx, dy;
+	double rho;
+	double cd;
+	double sin_c, cos_c, dy_sinc;
+	double lat, lon;
+
+
+	dy = p.y() - _falseNorthing;
+	dx = p.x() - _falseEasting;
+	rho = sqrt(dx * dx + dy * dy);
+	if (fabs(rho) <= 1.0e-10) {
+		lat = _latOrigin;
+		lon = _lonOrigin;
+	} else {
+		cd = rho / _ra;
+		sin_c = sin(cd);
+		cos_c = cos(cd);
+		dy_sinc = dy * sin_c;
+		lat = asin((cos_c * _sinLatOrigin) + ((dy_sinc * _cosLatOrigin) / rho));
+		if (fabs(_absLatOrigin - M_PI_2) < 1.0e-10) {
+			if (_latOrigin >= 0.0)
+				lon = _lonOrigin + atan2(dx, -dy);
+			else
+				lon = _lonOrigin + atan2(dx, dy);
+			}
+		else
+			lon = _lonOrigin + atan2((dx * sin_c), ((rho * _cosLatOrigin
+			  * cos_c) - (dy_sinc * _sinLatOrigin)));
+	}
+
+	if (lat > M_PI_2)
+		lat = M_PI_2;
+	else if (lat < -M_PI_2)
+		lat = -M_PI_2;
+
+	if (lon > M_PI)
+		lon -= 2.0 * M_PI;
+	if (lon < -M_PI)
+		lon += 2.0 * M_PI;
+
+	if (lon > M_PI)
+		lon = M_PI;
+	else if (lon < -M_PI)
+		lon = -M_PI;
+
+	return Coordinates(rad2deg(lon), rad2deg(lat));
+}

src/azimuthalequalarea.h

@@ -0,0 +1,30 @@
+#ifndef AZIMUTHALEQUALAREA_H
+#define AZIMUTHALEQUALAREA_H
+
+#include "ellipsoid.h"
+#include "projection.h"
+
+class AzimuthalEqualArea : public Projection
+{
+public:
+	AzimuthalEqualArea(const Ellipsoid &ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double falseEasting, double falseNorthing);
+
+	virtual QPointF ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const QPointF &p) const;
+
+private:
+	Ellipsoid _e;
+
+	double _ra;
+	double _sinLatOrigin;
+	double _cosLatOrigin;
+	double _absLatOrigin;
+
+	double _latOrigin;
+	double _lonOrigin;
+	double _falseNorthing;
+	double _falseEasting;
+};
+
+#endif // AZIMUTHALEQUALAREA_H

src/offlinemap.cpp

@@ -19,6 +19,7 @@
 #include "utm.h"
 #include "lambertconic.h"
 #include "albersequal.h"
+#include "azimuthalequalarea.h"
 #include "ozf.h"
 #include "rectc.h"
 #include "offlinemap.h"
@@ -208,6 +209,9 @@ bool OfflineMap::createProjection(const QString &datum,
 		_projection = new AlbersEqual(d.ellipsoid(), setup.standardParallel1,
 		  setup.standardParallel2, setup.latitudeOrigin, setup.longitudeOrigin,
 		  setup.falseEasting, setup.falseNorthing);
+	else if (projection == "(A)Lambert Azimuthual Equal Area")
+		_projection = new AzimuthalEqualArea(d.ellipsoid(), setup.latitudeOrigin,
+		  setup.longitudeOrigin, setup.falseEasting, setup.falseNorthing);
 	else if (projection == "(UTM) Universal Transverse Mercator") {
 		if (setup.zone)
 			_projection = new UTM(d.ellipsoid(), setup.zone);