diff --git a/gpxsee.pro b/gpxsee.pro index 6b91b145..904062a5 100644 --- a/gpxsee.pro +++ b/gpxsee.pro @@ -90,7 +90,8 @@ HEADERS += src/config.h \ src/transversemercator.h \ src/latlon.h \ src/utm.h \ - src/lambertconic.h + src/lambertconic.h \ + src/ellipsoid.h SOURCES += src/main.cpp \ src/gui.cpp \ src/poi.cpp \ @@ -154,7 +155,8 @@ SOURCES += src/main.cpp \ src/mercator.cpp \ src/transversemercator.cpp \ src/utm.cpp \ - src/lambertconic.cpp + src/lambertconic.cpp \ + src/ellipsoid.cpp RESOURCES += gpxsee.qrc TRANSLATIONS = lang/gpxsee_cs.ts \ lang/gpxsee_sv.ts diff --git a/src/ellipsoid.cpp b/src/ellipsoid.cpp new file mode 100644 index 00000000..38ec58e3 --- /dev/null +++ b/src/ellipsoid.cpp @@ -0,0 +1,92 @@ +#include +#include +#include "wgs84.h" +#include "ellipsoid.h" + + +#define INTERNATIONAL_RADIUS 6378388.0 +#define INTERNATIONAL_FLATTENING (1.0/297.0) +#define KRASSOVSKY_RADIUS 6378245.0 +#define KRASSOVSKY_FLATTENING (1.0/298.3) +#define BESSEL_RADIUS 6377397.155 +#define BESSEL_FLATTENING (1.0/299.1528128) +#define GRS80_RADIUS 6378137.0 +#define GRS80_FLATTENING (1.0/298.257222101) +#define WGS70_RADIUS 6378135.0 +#define WGS70_FLATTENING (1.0/298.26) + +#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 + +#define ARRAY_SIZE(array) (sizeof(array) / sizeof((array)[0])) + +typedef struct { + const char *name; + double radius; + double flattening; +} Entry; + +static Entry list[] = { + {"S42", KRASSOVSKY_RADIUS, KRASSOVSKY_FLATTENING}, + {"Pulkovo 1942", KRASSOVSKY_RADIUS, KRASSOVSKY_FLATTENING}, + {"European 1950", INTERNATIONAL_RADIUS, INTERNATIONAL_FLATTENING}, + {"European 1979", INTERNATIONAL_RADIUS, INTERNATIONAL_FLATTENING}, + {"NZGD1949", INTERNATIONAL_RADIUS, INTERNATIONAL_FLATTENING}, + {"NAD27", BESSEL_RADIUS, BESSEL_FLATTENING}, + {"NAD83", GRS80_RADIUS, GRS80_FLATTENING}, + {"WGS 72", WGS70_RADIUS, WGS70_FLATTENING} +}; + + +Ellipsoid::Ellipsoid() +{ + _radius = WGS84_RADIUS; + _flattening = WGS84_FLATTENING; +} + +Ellipsoid::Ellipsoid(const QString &datum) +{ + for (size_t i = 0; i < ARRAY_SIZE(list); i++) { + if (datum.startsWith(list[i].name)) { + _radius = list[i].radius; + _flattening = list[i].flattening; + return; + } + } + + _radius = WGS84_RADIUS; + _flattening = WGS84_FLATTENING; +} + +double Ellipsoid::q(double b) const +{ + double e = sqrt(_flattening * (2. - _flattening)); + double esb = e * sin(b); + return log(tan(M_PI_4 + b / 2.) * pow((1. - esb) / (1. + esb), e / 2.)); +} + +double Ellipsoid::iq(double q) const +{ + double e = sqrt(_flattening * (2. - _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; +} + +double Ellipsoid::nradius(double phi) const +{ + double sin_phi = sin(phi); + return (_radius / sqrt(1. - (_flattening * (2. - _flattening)) * sin_phi + * sin_phi)); +} diff --git a/src/ellipsoid.h b/src/ellipsoid.h new file mode 100644 index 00000000..5d2a2445 --- /dev/null +++ b/src/ellipsoid.h @@ -0,0 +1,24 @@ +#ifndef ELLIPSOID_H +#define ELLIPSOID_H + +class QString; + +class Ellipsoid +{ +public: + Ellipsoid(); + Ellipsoid(const QString &datum); + + double radius() const {return _radius;} + double flattening() const {return _flattening;} + + double q(double b) const; + double iq(double q) const; + double nradius(double phi) const; + +private: + double _radius; + double _flattening; +}; + +#endif // ELLIPSOID_H diff --git a/src/lambertconic.cpp b/src/lambertconic.cpp index 3bda15b1..ada62662 100644 --- a/src/lambertconic.cpp +++ b/src/lambertconic.cpp @@ -1,48 +1,11 @@ #include #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) +LambertConic::LambertConic(const Ellipsoid &ellipsoid, double standardParallel1, + double standardParallel2, double centralParallel, double centralMeridian, + double scale, double falseEasting, double falseNorthing) : _e(ellipsoid) { _cm = centralMeridian; _fe = falseEasting; @@ -51,12 +14,12 @@ LambertConic::LambertConic(double standardParallel1, double standardParallel2, double sp1 = deg2rad(standardParallel1); double sp2 = deg2rad(standardParallel2); - double N1 = nradius(sp1); - double N2 = nradius(sp2); + double N1 = _e.nradius(sp1); + double N2 = _e.nradius(sp2); - _q0 = q(deg2rad(centralParallel)); - double q1 = q(sp1); - double q2 = q(sp2); + _q0 = ellipsoid.q(deg2rad(centralParallel)); + double q1 = _e.q(sp1); + double q2 = _e.q(sp2); _n = log((N1 * cos(sp1)) / (N2 * cos(sp2))) / (q2 - q1); double R1 = N1 * cos(sp1) / _n; @@ -66,7 +29,7 @@ LambertConic::LambertConic(double standardParallel1, double standardParallel2, 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())))); + double R = _R0 * exp(_n * (_q0 - _e.q(deg2rad(c.lat())))); return QPointF(_fe + R * sin(dl), _fn + _R0 - R * cos(dl)); } @@ -79,5 +42,5 @@ Coordinates LambertConic::xy2ll(const QPointF &p) const 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))); + return Coordinates(rad2deg(deg2rad(_cm) + dl / _n), rad2deg(_e.iq(q))); } diff --git a/src/lambertconic.h b/src/lambertconic.h index 78ee1699..a065e302 100644 --- a/src/lambertconic.h +++ b/src/lambertconic.h @@ -1,19 +1,22 @@ #ifndef LAMBERTCONIC_H #define LAMBERTCONIC_H +#include "ellipsoid.h" #include "projection.h" class LambertConic : public Projection { public: - LambertConic(double standardParallel1, double standardParallel2, - double centralParallel, double centralMeridian, double scale, - double falseEasting, double falseNorthing); + LambertConic(const Ellipsoid &ellipsoid, 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: + Ellipsoid _e; + double _cm; double _fe; double _fn; diff --git a/src/offlinemap.cpp b/src/offlinemap.cpp index f4aec59b..01f90a94 100644 --- a/src/offlinemap.cpp +++ b/src/offlinemap.cpp @@ -12,6 +12,7 @@ #include "wgs84.h" #include "coordinates.h" #include "matrix.h" +#include "ellipsoid.h" #include "latlon.h" #include "mercator.h" #include "transversemercator.h" @@ -21,7 +22,7 @@ int OfflineMap::parseMapFile(QIODevice &device, QList &points, - QString &projection, ProjectionSetup &setup) + QString &projection, ProjectionSetup &setup, QString &datum) { bool res; int ln = 1; @@ -38,6 +39,8 @@ int OfflineMap::parseMapFile(QIODevice &device, QList &points, return ln; } else if (ln == 3) _imgPath = line.trimmed(); + else if (ln == 5) + datum = line.split(',').at(0).trimmed(); else { QList list = line.split(','); QString key(list.at(0).trimmed()); @@ -135,8 +138,9 @@ int OfflineMap::parseMapFile(QIODevice &device, QList &points, return 0; } -bool OfflineMap::createProjection(const QString &projection, - const ProjectionSetup &setup, QList &points) +bool OfflineMap::createProjection(const QString &datum, + const QString &projection, const ProjectionSetup &setup, + QList &points) { if (points.count() < 2) { qWarning("%s: insufficient number of reference points", @@ -147,12 +151,13 @@ bool OfflineMap::createProjection(const QString &projection, if (projection == "Mercator") _projection = new Mercator(); else if (projection == "Transverse Mercator") - _projection = new TransverseMercator(setup.centralMeridian, setup.scale, - setup.falseEasting, setup.falseNorthing); + _projection = new TransverseMercator(Ellipsoid(datum), + setup.centralMeridian, setup.scale, setup.falseEasting, + setup.falseNorthing); else if (projection == "Latitude/Longitude") _projection = new LatLon(); else if (projection == "Lambert Conformal Conic") - _projection = new LambertConic(setup.standardParallel1, + _projection = new LambertConic(Ellipsoid(datum), setup.standardParallel1, setup.standardParallel2, setup.centralParallel, setup.centralMeridian, setup.scale, setup.falseEasting, setup.falseNorthing); else if (projection == "(UTM) Universal Transverse Mercator") { @@ -335,7 +340,7 @@ OfflineMap::OfflineMap(const QString &path, QObject *parent) : Map(parent) { int errorLine = -2; QList points; - QString proj; + QString proj, datum; ProjectionSetup setup; @@ -362,7 +367,7 @@ OfflineMap::OfflineMap(const QString &path, QObject *parent) : Map(parent) if (tarFiles.at(j).endsWith(".map")) { QByteArray ba = _tar.file(tarFiles.at(j)); QBuffer buffer(&ba); - errorLine = parseMapFile(buffer, points, proj, setup); + errorLine = parseMapFile(buffer, points, proj, setup, datum); _imgPath = QString(); break; } @@ -370,14 +375,14 @@ OfflineMap::OfflineMap(const QString &path, QObject *parent) : Map(parent) break; } else if (fileName.endsWith(".map")) { QFile mapFile(mapFiles.at(i).absoluteFilePath()); - errorLine = parseMapFile(mapFile, points, proj, setup); + errorLine = parseMapFile(mapFile, points, proj, setup, datum); break; } } if (!mapLoaded(errorLine)) return; - if (!createProjection(proj, setup, points)) + if (!createProjection(datum, proj, setup, points)) return; if (!computeTransformation(points)) return; @@ -410,7 +415,7 @@ OfflineMap::OfflineMap(Tar &tar, const QString &path, QObject *parent) { int errorLine = -2; QList points; - QString proj; + QString proj, datum; ProjectionSetup setup; @@ -428,14 +433,14 @@ OfflineMap::OfflineMap(Tar &tar, const QString &path, QObject *parent) if (tarFiles.at(j).startsWith(prefix)) { QByteArray ba = tar.file(tarFiles.at(j)); QBuffer buffer(&ba); - errorLine = parseMapFile(buffer, points, proj, setup); + errorLine = parseMapFile(buffer, points, proj, setup, datum); break; } } if (!mapLoaded(errorLine)) return; - if (!createProjection(proj, setup, points)) + if (!createProjection(datum, proj, setup, points)) return; if (!totalSizeSet()) return; diff --git a/src/offlinemap.h b/src/offlinemap.h index 619e5674..138bfd26 100644 --- a/src/offlinemap.h +++ b/src/offlinemap.h @@ -62,10 +62,10 @@ private: } ProjectionSetup; int parseMapFile(QIODevice &device, QList &points, - QString &projection, ProjectionSetup &setup); + QString &projection, ProjectionSetup &setup, QString &datum); bool mapLoaded(int res); bool totalSizeSet(); - bool createProjection(const QString &projection, + bool createProjection(const QString &datum, const QString &projection, const ProjectionSetup &setup, QList &points); bool computeTransformation(const QList &points); bool computeResolution(QList &points); diff --git a/src/transversemercator.cpp b/src/transversemercator.cpp index 7e1b64dd..14f1b4d5 100644 --- a/src/transversemercator.cpp +++ b/src/transversemercator.cpp @@ -1,20 +1,22 @@ #include #include "rd.h" -#include "wgs84.h" +#include "ellipsoid.h" #include "transversemercator.h" -TransverseMercator::TransverseMercator(double centralMeridian, double scale, - double falseEasting, double falseNorthing) +TransverseMercator::TransverseMercator(const Ellipsoid &ellipsoid, + double centralMeridian, double scale, double falseEasting, + double falseNorthing) { _centralMeridian = centralMeridian; _scale = scale; _falseEasting = falseEasting; _falseNorthing = falseNorthing; - const double e2 = WGS84_FLATTENING * (2 - WGS84_FLATTENING); - const double n = WGS84_FLATTENING / (2 - WGS84_FLATTENING); - _rectifyingRadius = WGS84_RADIUS / (1 + n) + + const double e2 = ellipsoid.flattening() * (2 - ellipsoid.flattening()); + const double n = ellipsoid.flattening() / (2 - ellipsoid.flattening()); + _rectifyingRadius = ellipsoid.radius() / (1 + n) * (1 + 0.25*pow(n, 2) + 0.015625*pow(n, 4)); _A = e2; diff --git a/src/transversemercator.h b/src/transversemercator.h index d06d7b78..6ed1ee1b 100644 --- a/src/transversemercator.h +++ b/src/transversemercator.h @@ -3,11 +3,13 @@ #include "projection.h" +class Ellipsoid; + class TransverseMercator : public Projection { public: - TransverseMercator(double centralMeridian, double scale, - double falseEasting, double falseNorthing); + TransverseMercator(const Ellipsoid &ellipsoid, double centralMeridian, + double scale, double falseEasting, double falseNorthing); virtual QPointF ll2xy(const Coordinates &c) const; virtual Coordinates xy2ll(const QPointF &p) const; diff --git a/src/utm.cpp b/src/utm.cpp index 80af54e7..9b2a6294 100644 --- a/src/utm.cpp +++ b/src/utm.cpp @@ -1,6 +1,12 @@ +#include "ellipsoid.h" #include "utm.h" -UTM::UTM(const Coordinates &c) : _tm(0, 1.0, 0, 0) +UTM::UTM(int zone) : _tm(Ellipsoid(), (qAbs(zone) - 1)*6 - 180 + 3, 0.9996, + 500000, zone < 0 ? 10000000 : 0) +{ +} + +UTM::UTM(const Coordinates &c) : _tm(Ellipsoid(), 0, 0, 0, 0) { int zone = int((c.lon() + 180)/6) + 1; @@ -18,5 +24,6 @@ UTM::UTM(const Coordinates &c) : _tm(0, 1.0, 0, 0) } double cm = (zone - 1)*6 - 180 + 3; - _tm = TransverseMercator(cm, 0.9996, 500000, (c.lat() < 0) ? 10000000 : 0); + _tm = TransverseMercator(Ellipsoid(), cm, 0.9996, 500000, + (c.lat() < 0) ? 10000000 : 0); } diff --git a/src/utm.h b/src/utm.h index 5a7a3eb5..f368eb0e 100644 --- a/src/utm.h +++ b/src/utm.h @@ -7,8 +7,7 @@ class UTM : public Projection { public: - UTM(int zone) : _tm((qAbs(zone) - 1)*6 - 180 + 3, 0.9996, 500000, - zone < 0 ? 10000000 : 0) {} + UTM(int zone); UTM(const Coordinates &c); virtual QPointF ll2xy(const Coordinates &c) const