diff --git a/src/datum.cpp b/src/datum.cpp index c1e934a1..7c9dfd90 100644 --- a/src/datum.cpp +++ b/src/datum.cpp @@ -1,5 +1,7 @@ +#include #include #include "wgs84.h" +#include "rd.h" #include "datum.h" @@ -81,3 +83,39 @@ bool Datum::loadList(const QString &path) return true; } + +// Abridged Molodensky transformation +Coordinates Datum::toWGS84(const Coordinates &c) const +{ + if (_ellipsoid.radius() == WGS84_RADIUS + && _ellipsoid.flattening() == WGS84_FLATTENING + && _dx == 0 && _dy == 0 && _dz == 0) + return c; + + double rlat = deg2rad(c.lat()); + double rlon = deg2rad(c.lon()); + + double slat = sin(rlat); + double clat = cos(rlat); + double slon = sin(rlon); + double clon = cos(rlon); + double ssqlat = slat * slat; + + double from_f = ellipsoid().flattening(); + double df = WGS84_FLATTENING - from_f; + double from_a = ellipsoid().radius(); + double da = WGS84_RADIUS - from_a; + double from_esq = ellipsoid().flattening() + * (2.0 - ellipsoid().flattening()); + double adb = 1.0 / (1.0 - from_f); + double rn = from_a / sqrt(1 - from_esq * ssqlat); + double rm = from_a * (1 - from_esq) / pow((1 - from_esq * ssqlat), 1.5); + double from_h = 0.0; + + double dlat = (-dx() * slat * clon - dy() * slat * slon + dz() * clat + da + * rn * from_esq * slat * clat / from_a + df * (rm * adb + rn / adb) * slat + * clat) / (rm + from_h); + double dlon = (-dx() * slon + dy() * clon) / ((rn + from_h) * clat); + + return Coordinates(c.lon() + rad2deg(dlon), c.lat() + rad2deg(dlat)); +} diff --git a/src/datum.h b/src/datum.h index 128b003d..a49b5674 100644 --- a/src/datum.h +++ b/src/datum.h @@ -3,6 +3,7 @@ #include #include "ellipsoid.h" +#include "coordinates.h" class Datum { @@ -17,8 +18,8 @@ public: double dz() const {return _dz;} bool isNull() const {return _ellipsoid.isNull();} - bool isWGS84() const - {return _ellipsoid.isWGS84() && _dx == 0 && _dy == 0 && _dz == 0;} + + Coordinates toWGS84(const Coordinates &c) const; static bool loadList(const QString &path); static const QString &errorString() {return _errorString;} diff --git a/src/ellipsoid.h b/src/ellipsoid.h index c128ede6..dbdb680e 100644 --- a/src/ellipsoid.h +++ b/src/ellipsoid.h @@ -3,7 +3,6 @@ #include #include -#include "wgs84.h" class Ellipsoid { @@ -16,8 +15,6 @@ public: double flattening() const {return _flattening;} bool isNull() const {return _radius < 0 || _flattening < 0;} - bool isWGS84() const - {return _radius == WGS84_RADIUS && _flattening == WGS84_FLATTENING;} static bool loadList(const QString &path); static const QString &errorString() {return _errorString;} diff --git a/src/offlinemap.cpp b/src/offlinemap.cpp index b6e4016d..974954a8 100644 --- a/src/offlinemap.cpp +++ b/src/offlinemap.cpp @@ -23,39 +23,6 @@ #include "offlinemap.h" -// Abridged Molodensky transformation -static Coordinates toWGS84(Coordinates c, const Datum &datum) -{ - double rlat = deg2rad(c.lat()); - double rlon = deg2rad(c.lon()); - - double slat = sin(rlat); - double clat = cos(rlat); - double slon = sin(rlon); - double clon = cos(rlon); - double ssqlat = slat * slat; - - double from_f = datum.ellipsoid().flattening(); - double df = WGS84_FLATTENING - from_f; - double from_a = datum.ellipsoid().radius(); - 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); - double rn = from_a / sqrt(1 - from_esq * ssqlat); - double rm = from_a * (1 - from_esq) / pow((1 - from_esq * ssqlat), 1.5); - double from_h = 0.0; - - double dlat = (-datum.dx() * slat * clon - datum.dy() * slat * slon - + datum.dz() * clat + da * rn * from_esq * slat * clat / from_a + df - * (rm * adb + rn / adb) * slat * clat) / (rm + from_h); - double dlon = (-datum.dx() * slon + datum.dy() * clon) / ((rn + from_h) - * clat); - - return Coordinates(c.lon() + rad2deg(dlon), c.lat() + rad2deg(dlat)); -} - - int OfflineMap::parse(QIODevice &device, QList &points, QString &projection, ProjectionSetup &setup, QString &datum) { @@ -263,15 +230,10 @@ bool OfflineMap::createProjection(const QString &datum, } for (int i = 0; i < points.size(); i++) { - if (points.at(i).ll.isNull()) { - if (d.isWGS84()) - points[i].ll = _projection->xy2ll(points.at(i).pp); - else - points[i].ll = toWGS84(_projection->xy2ll(points.at(i).pp), d); - } else { - if (!d.isWGS84()) - points[i].ll = toWGS84(points[i].ll, d); - } + if (points.at(i).ll.isNull()) + points[i].ll = d.toWGS84(_projection->xy2ll(points.at(i).pp)); + else + points[i].ll = d.toWGS84(points[i].ll); } return true;