diff --git a/gpxsee.pro b/gpxsee.pro index 40d5ee21..f873a1ea 100644 --- a/gpxsee.pro +++ b/gpxsee.pro @@ -92,7 +92,8 @@ HEADERS += src/config.h \ src/ellipsoid.h \ src/ozf.h \ src/datum.h \ - src/maplist.h + src/maplist.h \ + src/albersequal.h SOURCES += src/main.cpp \ src/gui.cpp \ src/poi.cpp \ @@ -158,7 +159,8 @@ SOURCES += src/main.cpp \ src/ellipsoid.cpp \ src/ozf.cpp \ src/datum.cpp \ - src/maplist.cpp + src/maplist.cpp \ + src/albersequal.cpp RESOURCES += gpxsee.qrc TRANSLATIONS = lang/gpxsee_cs.ts \ lang/gpxsee_sv.ts \ diff --git a/src/albersequal.cpp b/src/albersequal.cpp new file mode 100644 index 00000000..4aefc4a5 --- /dev/null +++ b/src/albersequal.cpp @@ -0,0 +1,190 @@ +#include "ellipsoid.h" +#include "rd.h" +#include "albersequal.h" + + +#ifndef M_PI_2 + #define M_PI_2 1.57079632679489661923 +#endif // M_PI_2 + +#define ONE_MINUS_SQR(x) (1.0 - (x) * (x)) +#define ALBERS_Q(slat, one_minus_sqr_es_sin, es_sin) \ + (_one_minus_es2 * ((slat) / (one_minus_sqr_es_sin) - \ + (1 / (_two_es)) * log((1 - (es_sin)) / (1 + (es_sin))))) +#define ALBERS_M(clat, one_minus_sqr_es_sin) \ + ((clat) / sqrt(one_minus_sqr_es_sin)) + + +AlbersEqual::AlbersEqual(const Ellipsoid &ellipsoid, double standardParallel1, + double standardParallel2, double latitudeOrigin, double longitudeOrigin, + double falseEasting, double falseNorthing) +{ + double sin_lat, sin_lat1, sin_lat2, cos_lat1, cos_lat2; + double m1, m2, sqr_m1, sqr_m2; + double q0, q1, q2; + double es_sin, es_sin1, es_sin2; + double one_minus_sqr_es_sin1, one_minus_sqr_es_sin2; + double nq0; + double sp1, sp2; + + + _e = ellipsoid; + _latitudeOrigin = deg2rad(latitudeOrigin); + _longitudeOrigin = deg2rad(longitudeOrigin); + _falseEasting = falseEasting; + _falseNorthing = falseNorthing; + + sp1 = deg2rad(standardParallel1); + sp2 = deg2rad(standardParallel2); + + if (_latitudeOrigin > M_PI) + _latitudeOrigin -= M_PI * 2.0; + + _es2 = 2 * _e.flattening() - _e.flattening() * _e.flattening(); + _es = sqrt(_es2); + _one_minus_es2 = 1 - _es2; + _two_es = 2 * _es; + + sin_lat = sin(_latitudeOrigin); + es_sin = _es * sin_lat; + q0 = ALBERS_Q(sin_lat, ONE_MINUS_SQR(es_sin), es_sin); + + sin_lat1 = sin(sp1); + cos_lat1 = cos(sp1); + es_sin1 = _es * sin_lat1; + one_minus_sqr_es_sin1 = ONE_MINUS_SQR(es_sin1); + m1 = ALBERS_M(cos_lat1, one_minus_sqr_es_sin1); + q1 = ALBERS_Q(sin_lat1, one_minus_sqr_es_sin1, es_sin1); + + sqr_m1 = m1 * m1; + if (fabs(sp1 - sp2) > 1.0e-10) { + sin_lat2 = sin(sp2); + cos_lat2 = cos(sp2); + es_sin2 = _es * sin_lat2; + one_minus_sqr_es_sin2 = ONE_MINUS_SQR(es_sin2); + m2 = ALBERS_M(cos_lat2, one_minus_sqr_es_sin2); + q2 = ALBERS_Q(sin_lat2, one_minus_sqr_es_sin2, es_sin2); + sqr_m2 = m2 * m2; + _n = (sqr_m1 - sqr_m2) / (q2 - q1); + } else + _n = sin_lat1; + + _C = sqr_m1 + _n * q1; + _a_over_n = _e.radius() / _n; + nq0 = _n * q0; + _rho0 = (_C < nq0) ? 0 : _a_over_n * sqrt(_C - nq0); +} + +QPointF AlbersEqual::ll2xy(const Coordinates &c) const +{ + double dlam; + double sin_lat; + double es_sin; + double q; + double rho; + double theta; + double nq; + + + dlam = deg2rad(c.lon()) - _longitudeOrigin; + if (dlam > M_PI) + dlam -= 2.0 * M_PI; + if (dlam < -M_PI) + dlam += 2.0 * M_PI; + + sin_lat = sin(deg2rad(c.lat())); + es_sin = _es * sin_lat; + q = ALBERS_Q(sin_lat, ONE_MINUS_SQR(es_sin), es_sin); + nq = _n * q; + rho = (_C < nq) ? 0 : _a_over_n * sqrt(_C - nq); + theta = _n * dlam; + + return QPointF(rho * sin(theta) + _falseEasting, + _rho0 - rho * cos(theta) + _falseNorthing); +} + +Coordinates AlbersEqual::xy2ll(const QPointF &p) const +{ + double dy, dx; + double rho0_minus_dy; + double q, qc, q_over_2; + double rho, rho_n; + double phi, delta_phi = 1.0; + double sin_phi; + double es_sin, one_minus_sqr_es_sin; + double theta = 0.0; + int count = 30; + double tolerance = 4.85e-10; + double lat, lon; + + + dy = p.y() - _falseNorthing; + dx = p.x() - _falseEasting; + + rho0_minus_dy = _rho0 - dy; + rho = sqrt(dx * dx + rho0_minus_dy * rho0_minus_dy); + + if (_n < 0) { + rho *= -1.0; + dy *= -1.0; + dx *= -1.0; + rho0_minus_dy *= -1.0; + } + + if (rho != 0.0) + theta = atan2(dx, rho0_minus_dy); + rho_n = rho * _n; + q = (_C - (rho_n * rho_n) / (_e.radius() * _e.radius())) / _n; + qc = 1 - ((_one_minus_es2) / (_two_es)) * log((1.0 - _es) / (1.0 + _es)); + if (fabs(fabs(qc) - fabs(q)) > 1.0e-6) { + q_over_2 = q / 2.0; + if (q_over_2 > 1.0) + lat = M_PI_2; + else if (q_over_2 < -1.0) + lat = -M_PI_2; + else { + phi = asin(q_over_2); + if (_es < 1.0e-10) + lat = phi; + else { + while ((fabs(delta_phi) > tolerance) && count) { + sin_phi = sin(phi); + es_sin = _es * sin_phi; + one_minus_sqr_es_sin = ONE_MINUS_SQR(es_sin); + delta_phi = (one_minus_sqr_es_sin * one_minus_sqr_es_sin) + / (2.0 * cos(phi)) * (q / (_one_minus_es2) - sin_phi + / one_minus_sqr_es_sin + (log((1.0 - es_sin) + / (1.0 + es_sin)) / (_two_es))); + phi += delta_phi; + count --; + } + + lat = phi; + } + + if (lat > M_PI_2) + lat = M_PI_2; + else if (lat < -M_PI_2) + lat = -M_PI_2; + } + } else { + if (q >= 0.0) + lat = M_PI_2; + else + lat = -M_PI_2; + } + + lon = _longitudeOrigin + theta / _n; + + if (lon > M_PI) + lon -= M_PI * 2; + if (lon < -M_PI) + lon += M_PI * 2; + + 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/albersequal.h b/src/albersequal.h new file mode 100644 index 00000000..c510a652 --- /dev/null +++ b/src/albersequal.h @@ -0,0 +1,36 @@ +#ifndef ALBERSEQUAL_H +#define ALBERSEQUAL_H + +#include "projection.h" + +class Ellipsoid; + +class AlbersEqual : public Projection +{ +public: + AlbersEqual(const Ellipsoid &ellipsoid, double standardParallel1, + double standardParallel2, 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 _latitudeOrigin; + double _longitudeOrigin; + double _falseEasting; + double _falseNorthing; + + double _rho0; + double _C; + double _n; + double _es; + double _es2; + double _a_over_n; + double _one_minus_es2; + double _two_es; +}; + +#endif // ALBERSEQUAL_H diff --git a/src/offlinemap.cpp b/src/offlinemap.cpp index af970652..5469007f 100644 --- a/src/offlinemap.cpp +++ b/src/offlinemap.cpp @@ -18,6 +18,7 @@ #include "transversemercator.h" #include "utm.h" #include "lambertconic.h" +#include "albersequal.h" #include "ozf.h" #include "offlinemap.h" @@ -25,8 +26,6 @@ // Abridged Molodensky transformation static Coordinates toWGS84(Coordinates c, const Datum &datum) { - Ellipsoid WGS84(WGS84_RADIUS, WGS84_FLATTENING); - double dX = datum.dx(); double dY = datum.dy(); double dZ = datum.dz(); @@ -38,9 +37,9 @@ static Coordinates toWGS84(Coordinates c, const Datum &datum) double ssqlat = slat * slat; double from_f = datum.ellipsoid().flattening(); - double df = WGS84.flattening() - from_f; + double df = WGS84_FLATTENING - from_f; double from_a = datum.ellipsoid().radius(); - double da = WGS84.radius() - from_a; + double da = WGS84_RADIUS - from_a; double from_esq = datum.ellipsoid().flattening() * (2.0 - datum.ellipsoid().flattening()); double adb = 1.0 / (1.0 - from_f); @@ -219,6 +218,10 @@ bool OfflineMap::createProjection(const QString &datum, setup.standardParallel1, setup.standardParallel2, setup.latitudeOrigin, setup.longitudeOrigin, setup.scale, setup.falseEasting, setup.falseNorthing); + else if (projection == "Albers Equal Area") + _projection = new AlbersEqual(d.ellipsoid(), setup.standardParallel1, + setup.standardParallel2, 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);