diff --git a/gpxsee.pro b/gpxsee.pro index d0227c07..7981f0f6 100644 --- a/gpxsee.pro +++ b/gpxsee.pro @@ -89,7 +89,8 @@ HEADERS += src/config.h \ src/mercator.h \ src/transversemercator.h \ src/latlon.h \ - src/utm.h + src/utm.h \ + src/lambertconic.h SOURCES += src/main.cpp \ src/gui.cpp \ src/poi.cpp \ @@ -152,7 +153,8 @@ SOURCES += src/main.cpp \ src/atlas.cpp \ src/mercator.cpp \ src/transversemercator.cpp \ - src/utm.cpp + src/utm.cpp \ + src/lambertconic.cpp RESOURCES += gpxsee.qrc TRANSLATIONS = lang/gpxsee_cs.ts \ lang/gpxsee_sv.ts diff --git a/src/lambertconic.cpp b/src/lambertconic.cpp new file mode 100644 index 00000000..3bda15b1 --- /dev/null +++ b/src/lambertconic.cpp @@ -0,0 +1,83 @@ +#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) +{ + _cm = centralMeridian; + _fe = falseEasting; + _fn = falseNorthing; + + double sp1 = deg2rad(standardParallel1); + double sp2 = deg2rad(standardParallel2); + + double N1 = nradius(sp1); + double N2 = nradius(sp2); + + _q0 = q(deg2rad(centralParallel)); + double q1 = q(sp1); + double q2 = q(sp2); + + _n = log((N1 * cos(sp1)) / (N2 * cos(sp2))) / (q2 - q1); + double R1 = N1 * cos(sp1) / _n; + _R0 = scale * R1 * exp(_n * (q1 - _q0)); +} + +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())))); + + return QPointF(_fe + R * sin(dl), _fn + _R0 - R * cos(dl)); +} + +Coordinates LambertConic::xy2ll(const QPointF &p) const +{ + double dl = atan((p.x() - _fe) / (_fn + _R0 - p.y())); + double dx = p.x() - _fe; + double dy = p.y() - _fn - _R0; + 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))); +} diff --git a/src/lambertconic.h b/src/lambertconic.h new file mode 100644 index 00000000..78ee1699 --- /dev/null +++ b/src/lambertconic.h @@ -0,0 +1,26 @@ +#ifndef LAMBERTCONIC_H +#define LAMBERTCONIC_H + +#include "projection.h" + +class LambertConic : public Projection +{ +public: + LambertConic(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: + double _cm; + double _fe; + double _fn; + + double _q0; + double _R0; + double _n; +}; + +#endif // LAMBERTCONIC_H diff --git a/src/offlinemap.cpp b/src/offlinemap.cpp index 5d31ef3c..f4aec59b 100644 --- a/src/offlinemap.cpp +++ b/src/offlinemap.cpp @@ -16,6 +16,7 @@ #include "mercator.h" #include "transversemercator.h" #include "utm.h" +#include "lambertconic.h" #include "offlinemap.h" @@ -104,18 +105,27 @@ 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); + if (!res) + setup.centralParallel = 0; setup.centralMeridian = list.at(2).trimmed().toFloat(&res); if (!res) setup.centralMeridian = 0; setup.scale = list.at(3).trimmed().toFloat(&res); if (!res) - setup.scale = 0; + setup.scale = 1.0; setup.falseEasting = list.at(4).trimmed().toFloat(&res); if (!res) setup.falseEasting = 0; setup.falseNorthing = list.at(5).trimmed().toFloat(&res); if (!res) setup.falseNorthing = 0; + setup.standardParallel1 = list.at(6).trimmed().toFloat(&res); + if (!res) + setup.standardParallel1 = 0; + setup.standardParallel2 = list.at(7).trimmed().toFloat(&res); + if (!res) + setup.standardParallel2 = 0; } } @@ -141,6 +151,10 @@ bool OfflineMap::createProjection(const QString &projection, setup.falseEasting, setup.falseNorthing); else if (projection == "Latitude/Longitude") _projection = new LatLon(); + else if (projection == "Lambert Conformal Conic") + _projection = new LambertConic(setup.standardParallel1, + setup.standardParallel2, setup.centralParallel, setup.centralMeridian, + setup.scale, setup.falseEasting, setup.falseNorthing); else if (projection == "(UTM) Universal Transverse Mercator") { if (setup.zone) _projection = new UTM(setup.zone); diff --git a/src/offlinemap.h b/src/offlinemap.h index 395032bf..619e5674 100644 --- a/src/offlinemap.h +++ b/src/offlinemap.h @@ -51,10 +51,13 @@ private: } ReferencePoint; typedef struct { + double centralParallel; double centralMeridian; double scale; double falseEasting; double falseNorthing; + double standardParallel1; + double standardParallel2; int zone; } ProjectionSetup;