diff --git a/src/ellipsoid.cpp b/src/ellipsoid.cpp index 5d61fb9a..b28c299a 100644 --- a/src/ellipsoid.cpp +++ b/src/ellipsoid.cpp @@ -12,85 +12,31 @@ #define CLARKE1866_FLATTENING (1.0/294.9786982) #define GRS80_RADIUS 6378137.0 #define GRS80_FLATTENING (1.0/298.257222101) -#define WGS70_RADIUS 6378135.0 -#define WGS70_FLATTENING (1.0/298.26) -#define SAD69_RADIUS 6378160.0 -#define SAD69_FLATTENING (1.0/298.25) +#define WGS72_RADIUS 6378135.0 +#define WGS72_FLATTENING (1.0/298.26) +#define GRS67_RADIUS 6378160.0 +#define GRS67_FLATTENING (1.0/298.25) -#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", CLARKE1866_RADIUS, CLARKE1866_FLATTENING}, - {"NAD83", GRS80_RADIUS, GRS80_FLATTENING}, - {"WGS 72", WGS70_RADIUS, WGS70_FLATTENING}, - {"South American 1969", SAD69_RADIUS, SAD69_FLATTENING} +// Must be in Ellipsoid::Name order! +struct {double radius; double flattening;} static ellipsoids[] = { + {CLARKE1866_RADIUS, CLARKE1866_FLATTENING}, + {GRS80_RADIUS, GRS80_FLATTENING}, + {INTERNATIONAL_RADIUS, INTERNATIONAL_FLATTENING}, + {KRASSOVSKY_RADIUS, KRASSOVSKY_FLATTENING}, + {WGS84_RADIUS, WGS84_FLATTENING}, + {WGS72_RADIUS, WGS72_FLATTENING}, + {GRS67_RADIUS, GRS67_FLATTENING} }; - Ellipsoid::Ellipsoid() { _radius = WGS84_RADIUS; _flattening = WGS84_FLATTENING; } -Ellipsoid::Ellipsoid(const QString &datum) +Ellipsoid::Ellipsoid(Name name) { - 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)); + _radius = ellipsoids[name].radius; + _flattening = ellipsoids[name].flattening; } diff --git a/src/ellipsoid.h b/src/ellipsoid.h index 5d2a2445..75a55474 100644 --- a/src/ellipsoid.h +++ b/src/ellipsoid.h @@ -1,21 +1,25 @@ #ifndef ELLIPSOID_H #define ELLIPSOID_H -class QString; - class Ellipsoid { public: + enum Name { + Clarke1866, + GRS80, + International1924, + Krassowsky1940, + WGS84, + WGS72, + GRS67 + }; + Ellipsoid(); - Ellipsoid(const QString &datum); + Ellipsoid(Name name); 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; diff --git a/src/lambertconic.cpp b/src/lambertconic.cpp index ada62662..618d8f1f 100644 --- a/src/lambertconic.cpp +++ b/src/lambertconic.cpp @@ -2,24 +2,59 @@ #include "rd.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(const Ellipsoid &el, double b) +{ + double e = sqrt(el.flattening() * (2. - el.flattening())); + double esb = e * sin(b); + return log(tan(M_PI_4 + b / 2.) * pow((1. - esb) / (1. + esb), e / 2.)); +} + +static double iq(const Ellipsoid &el, double q) +{ + double e = sqrt(el.flattening() * (2. - el.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(const Ellipsoid &el, double phi) +{ + double sin_phi = sin(phi); + return (el.radius() / sqrt(1. - (el.flattening() * (2. - el.flattening())) + * sin_phi * sin_phi)); +} LambertConic::LambertConic(const Ellipsoid &ellipsoid, double standardParallel1, - double standardParallel2, double centralParallel, double centralMeridian, + double standardParallel2, double latitudeOrigin, double longitudeOrigin, double scale, double falseEasting, double falseNorthing) : _e(ellipsoid) { - _cm = centralMeridian; + _cm = longitudeOrigin; _fe = falseEasting; _fn = falseNorthing; double sp1 = deg2rad(standardParallel1); double sp2 = deg2rad(standardParallel2); - double N1 = _e.nradius(sp1); - double N2 = _e.nradius(sp2); + double N1 = nradius(_e, sp1); + double N2 = nradius(_e, sp2); - _q0 = ellipsoid.q(deg2rad(centralParallel)); - double q1 = _e.q(sp1); - double q2 = _e.q(sp2); + _q0 = q(_e, deg2rad(latitudeOrigin)); + double q1 = q(_e, sp1); + double q2 = q(_e, sp2); _n = log((N1 * cos(sp1)) / (N2 * cos(sp2))) / (q2 - q1); double R1 = N1 * cos(sp1) / _n; @@ -29,7 +64,7 @@ LambertConic::LambertConic(const Ellipsoid &ellipsoid, double standardParallel1, QPointF LambertConic::ll2xy(const Coordinates &c) const { double dl = _n * (deg2rad(c.lon()) - deg2rad(_cm)); - double R = _R0 * exp(_n * (_q0 - _e.q(deg2rad(c.lat())))); + double R = _R0 * exp(_n * (_q0 - q(_e, deg2rad(c.lat())))); return QPointF(_fe + R * sin(dl), _fn + _R0 - R * cos(dl)); } @@ -42,5 +77,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(_e.iq(q))); + return Coordinates(rad2deg(deg2rad(_cm) + dl / _n), rad2deg(iq(_e, q))); } diff --git a/src/lambertconic.h b/src/lambertconic.h index a065e302..c331d8e4 100644 --- a/src/lambertconic.h +++ b/src/lambertconic.h @@ -8,7 +8,7 @@ class LambertConic : public Projection { public: LambertConic(const Ellipsoid &ellipsoid, double standardParallel1, - double standardParallel2, double centralParallel, double centralMeridian, + double standardParallel2, double latitudeOrigin, double longitudeOrigin, double scale, double falseEasting, double falseNorthing); virtual QPointF ll2xy(const Coordinates &c) const; diff --git a/src/misc.h b/src/misc.h index c4c16ac9..307ee373 100644 --- a/src/misc.h +++ b/src/misc.h @@ -3,6 +3,9 @@ #include +#define ARRAY_SIZE(array) \ + (sizeof(array) / sizeof((array)[0])) + double niceNum(double x, int round); int str2int(const char *str, int len); QRectF scaled(const QRectF &rect, qreal factor); diff --git a/src/offlinemap.cpp b/src/offlinemap.cpp index 18ea9eed..e7d623d3 100644 --- a/src/offlinemap.cpp +++ b/src/offlinemap.cpp @@ -21,6 +21,18 @@ #include "offlinemap.h" +struct {const char *name; Ellipsoid::Name ellipsoid;} static datums[] = { + {"S42", Ellipsoid::Krassowsky1940}, + {"Pulkovo 1942", Ellipsoid::Krassowsky1940}, + {"European 1950", Ellipsoid::International1924}, + {"European 1979", Ellipsoid::International1924}, + {"NZGD1949", Ellipsoid::International1924}, + {"NAD27", Ellipsoid::Clarke1866}, + {"NAD83", Ellipsoid::GRS80}, + {"WGS 72", Ellipsoid::WGS72}, + {"South American 1969", Ellipsoid::GRS67} +}; + int OfflineMap::parseMapFile(QIODevice &device, QList &points, QString &projection, ProjectionSetup &setup, QString &datum) { @@ -108,12 +120,12 @@ int OfflineMap::parseMapFile(QIODevice &device, QList &points, } else if (key == "Projection Setup") { if (list.count() < 8) return ln; - setup.centralParallel = list.at(1).trimmed().toFloat(&res); + setup.latitudeOrigin = list.at(1).trimmed().toFloat(&res); if (!res) - setup.centralParallel = 0; - setup.centralMeridian = list.at(2).trimmed().toFloat(&res); + setup.latitudeOrigin = 0; + setup.longitudeOrigin = list.at(2).trimmed().toFloat(&res); if (!res) - setup.centralMeridian = 0; + setup.longitudeOrigin = 0; setup.scale = list.at(3).trimmed().toFloat(&res); if (!res) setup.scale = 1.0; @@ -142,29 +154,39 @@ bool OfflineMap::createProjection(const QString &datum, const QString &projection, const ProjectionSetup &setup, QList &points) { + Ellipsoid::Name ellipsoid = Ellipsoid::WGS84; + if (points.count() < 2) { qWarning("%s: insufficient number of reference points", qPrintable(_name)); return false; } + for (size_t i = 0; i < ARRAY_SIZE(datums); i++) { + if (datum.startsWith(datums[i].name)) { + ellipsoid = datums[i].ellipsoid; + break; + } + } + if (projection == "Mercator") _projection = new Mercator(); else if (projection == "Transverse Mercator") - _projection = new TransverseMercator(Ellipsoid(datum), - setup.centralMeridian, setup.scale, setup.falseEasting, + _projection = new TransverseMercator(Ellipsoid(ellipsoid), + setup.longitudeOrigin, setup.scale, setup.falseEasting, setup.falseNorthing); else if (projection == "Latitude/Longitude") _projection = new LatLon(); else if (projection == "Lambert Conformal Conic") - _projection = new LambertConic(Ellipsoid(datum), setup.standardParallel1, - setup.standardParallel2, setup.centralParallel, setup.centralMeridian, - setup.scale, setup.falseEasting, setup.falseNorthing); + _projection = new LambertConic(Ellipsoid(ellipsoid), + setup.standardParallel1, setup.standardParallel2, + setup.latitudeOrigin, setup.longitudeOrigin, setup.scale, + setup.falseEasting, setup.falseNorthing); else if (projection == "(UTM) Universal Transverse Mercator") { if (setup.zone) - _projection = new UTM(setup.zone); + _projection = new UTM(Ellipsoid(ellipsoid), setup.zone); else if (!points.first().ll.isNull()) - _projection = new UTM(points.first().ll); + _projection = new UTM(Ellipsoid(ellipsoid), points.first().ll); else { qWarning("%s: Can not determine UTM zone", qPrintable(_name)); return false; @@ -496,36 +518,6 @@ void OfflineMap::unload() } } -QRectF OfflineMap::bounds() const -{ - return QRectF(QPointF(0, 0), _size); -} - -qreal OfflineMap::zoomFit(const QSize &size, const QRectF &br) -{ - Q_UNUSED(size); - Q_UNUSED(br); - - return 1.0; -} - -qreal OfflineMap::resolution(const QPointF &p) const -{ - Q_UNUSED(p); - - return _resolution; -} - -qreal OfflineMap::zoomIn() -{ - return 1.0; -} - -qreal OfflineMap::zoomOut() -{ - return 1.0; -} - void OfflineMap::draw(QPainter *painter, const QRectF &rect) { if (_tileSize.isValid()) { diff --git a/src/offlinemap.h b/src/offlinemap.h index e0726f16..1f2d8022 100644 --- a/src/offlinemap.h +++ b/src/offlinemap.h @@ -22,13 +22,13 @@ public: const QString &name() const {return _name;} - QRectF bounds() const; - qreal resolution(const QPointF &p) const; + QRectF bounds() const {return QRectF(QPointF(0, 0), _size);} + qreal resolution(const QPointF &) const {return _resolution;} qreal zoom() const {return 1.0;} - qreal zoomFit(const QSize &size, const QRectF &br); - qreal zoomIn(); - qreal zoomOut(); + qreal zoomFit(const QSize &, const QRectF &) {return 1.0;} + qreal zoomIn() {return 1.0;} + qreal zoomOut() {return 1.0;} QPointF ll2xy(const Coordinates &c) const {return _transform.map(_projection->ll2xy(c));} @@ -57,8 +57,8 @@ private: } ReferencePoint; typedef struct { - double centralParallel; - double centralMeridian; + double latitudeOrigin; + double longitudeOrigin; double scale; double falseEasting; double falseNorthing; diff --git a/src/utm.cpp b/src/utm.cpp index 9b2a6294..aef6cfe7 100644 --- a/src/utm.cpp +++ b/src/utm.cpp @@ -1,12 +1,14 @@ #include "ellipsoid.h" #include "utm.h" -UTM::UTM(int zone) : _tm(Ellipsoid(), (qAbs(zone) - 1)*6 - 180 + 3, 0.9996, - 500000, zone < 0 ? 10000000 : 0) +UTM::UTM(const Ellipsoid &ellipsoid, 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) +UTM::UTM(const Ellipsoid &ellipsoid, const Coordinates &c) + : _tm(ellipsoid, 0, 0, 0, 0) { int zone = int((c.lon() + 180)/6) + 1; @@ -24,6 +26,6 @@ UTM::UTM(const Coordinates &c) : _tm(Ellipsoid(), 0, 0, 0, 0) } double cm = (zone - 1)*6 - 180 + 3; - _tm = TransverseMercator(Ellipsoid(), cm, 0.9996, 500000, + _tm = TransverseMercator(ellipsoid, cm, 0.9996, 500000, (c.lat() < 0) ? 10000000 : 0); } diff --git a/src/utm.h b/src/utm.h index f368eb0e..3b989952 100644 --- a/src/utm.h +++ b/src/utm.h @@ -7,8 +7,8 @@ class UTM : public Projection { public: - UTM(int zone); - UTM(const Coordinates &c); + UTM(const Ellipsoid &ellipsoid, int zone); + UTM(const Ellipsoid &ellipsoid, const Coordinates &c); virtual QPointF ll2xy(const Coordinates &c) const {return _tm.ll2xy(c);}