diff --git a/gpxsee.pro b/gpxsee.pro index 50a26fc1..472a1aa2 100644 --- a/gpxsee.pro +++ b/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 \ diff --git a/src/azimuthalequalarea.cpp b/src/azimuthalequalarea.cpp new file mode 100644 index 00000000..7ba60fe3 --- /dev/null +++ b/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)); +} diff --git a/src/azimuthalequalarea.h b/src/azimuthalequalarea.h new file mode 100644 index 00000000..77048dc6 --- /dev/null +++ b/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 diff --git a/src/offlinemap.cpp b/src/offlinemap.cpp index d77da929..e356dc6f 100644 --- a/src/offlinemap.cpp +++ b/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);