From abd1dc2450c9b3998fe3507ab557453db7cf2287 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20T=C5=AFma?= Date: Fri, 4 Dec 2020 00:25:57 +0100 Subject: [PATCH] Added support for polyconic projections --- gpxsee.pro | 2 + pkg/csv/pcs.csv | 4 + src/map/bsbmap.cpp | 3 + src/map/geotiff.cpp | 3 + src/map/mapfile.cpp | 2 + src/map/polyconic.cpp | 211 +++++++++++++++++++++++++++++++++++++++++ src/map/polyconic.h | 37 ++++++++ src/map/projection.cpp | 7 ++ 8 files changed, 269 insertions(+) create mode 100644 src/map/polyconic.cpp create mode 100644 src/map/polyconic.h diff --git a/gpxsee.pro b/gpxsee.pro index 66409362..0d668643 100644 --- a/gpxsee.pro +++ b/gpxsee.pro @@ -106,6 +106,7 @@ HEADERS += src/common/config.h \ src/map/IMG/textpathitem.h \ src/map/IMG/textpointitem.h \ src/map/bsbmap.h \ + src/map/polyconic.h \ src/map/projection.h \ src/map/ellipsoid.h \ src/map/datum.h \ @@ -281,6 +282,7 @@ SOURCES += src/main.cpp \ src/map/downloader.cpp \ src/map/emptymap.cpp \ src/map/ozimap.cpp \ + src/map/polyconic.cpp \ src/map/tar.cpp \ src/map/atlas.cpp \ src/map/ozf.cpp \ diff --git a/pkg/csv/pcs.csv b/pkg/csv/pcs.csv index 7d195714..2bf845e1 100644 --- a/pkg/csv/pcs.csv +++ b/pkg/csv/pcs.csv @@ -1887,6 +1887,7 @@ CI1971 / Chatham Islands Map Grid,5518,4672,5517,9001,9807,4500,8801,-44,9110,88 CI1979 / Chatham Islands Map Grid,5519,4673,5517,9001,9807,4500,8801,-44,9110,8802,-176.3,9110,8805,1,9201,8806,350000,9001,8807,650000,9001,,,,,, DHDN / 3-degree Gauss-Kruger zone 1,5520,4314,16261,9001,9807,4530,8801,0,9102,8802,3,9102,8805,1,9201,8806,1500000,9001,8807,0,9001,,,,,, WGS 84 / Gabon TM 2011,5523,4326,5522,9001,9807,4499,8801,0,9102,8802,11.3,9110,8805,0.9996,9201,8806,1500000,9001,8807,5500000,9001,,,,,, +SAD69(96) / Brazil Polyconic,5530,5527,19941,9001,9818,4499,8801,0,9102,8802,-54,9102,8806,5000000,9001,8807,10000000,9001,,,,,,,,, SAD69(96) / UTM zone 21S,5531,5527,16121,9001,9807,4400,8801,0,9102,8802,-57,9102,8805,0.9996,9201,8806,500000,9001,8807,10000000,9001,,,,,, SAD69(96) / UTM zone 22S,5532,4618,16122,9001,9807,4400,8801,0,9102,8802,-51,9102,8805,0.9996,9201,8806,500000,9001,8807,10000000,9001,,,,,, SAD69(96) / UTM zone 23S,5533,5527,16123,9001,9807,4400,8801,0,9102,8802,-45,9102,8805,0.9996,9201,8806,500000,9001,8807,10000000,9001,,,,,, @@ -1970,6 +1971,7 @@ SAD69(96) / UTM zone 18S,5875,5527,16118,9001,9807,4400,8801,0,9102,8802,-75,910 SAD69(96) / UTM zone 19S,5876,5527,16119,9001,9807,4400,8801,0,9102,8802,-69,9102,8805,0.9996,9201,8806,500000,9001,8807,10000000,9001,,,,,, SAD69(96) / UTM zone 20S,5877,5527,16120,9001,9807,4400,8801,0,9102,8802,-63,9102,8805,0.9996,9201,8806,500000,9001,8807,10000000,9001,,,,,, Cadastre 1997 / UTM zone 38S,5879,4475,16138,9001,9807,4400,8801,0,9102,8802,45,9102,8805,0.9996,9201,8806,500000,9001,8807,10000000,9001,,,,,, +SIRGAS 2000 / Brazil Polyconic,5880,4674,19941,9001,9818,4499,8801,0,9102,8802,-54,9102,8806,5000000,9001,8807,10000000,9001,,,,,,,,, WGS 84 / EPSG Arctic Regional zone A1,5921,4326,5906,9001,9802,4400,8821,81.19020136,9110,8822,-111,9102,8823,85,9102,8824,77,9102,8826,0,9001,8827,0,9001,,, WGS 84 / EPSG Arctic Regional zone A2,5922,4326,5907,9001,9802,4400,8821,81.19020136,9110,8822,-39,9102,8823,85,9102,8824,77,9102,8826,0,9001,8827,0,9001,,, WGS 84 / EPSG Arctic Regional zone A3,5923,4326,5908,9001,9802,4400,8821,81.19020136,9110,8822,33,9102,8823,85,9102,8824,77,9102,8826,0,9001,8827,0,9001,,, @@ -3046,6 +3048,8 @@ Pulkovo 1942 / Gauss-Kruger 32N,28492,4284,16332,9001,9807,4530,8801,0,9102,8802 Qatar 1974 / Qatar National Grid,28600,4285,19919,9001,9807,4400,8801,24.27,9110,8802,51.13,9110,8805,0.99999,9201,8806,200000,9001,8807,300000,9001,,,,,, Amersfoort / RD Old,28991,4289,19913,9001,9809,4499,8801,52.0922178,9110,8802,5.23155,9110,8805,0.9999079,9201,8806,0,9001,8807,0,9001,,,,,, Amersfoort / RD New,28992,4289,19914,9001,9809,4499,8801,52.0922178,9110,8802,5.23155,9110,8805,0.9999079,9201,8806,155000,9001,8807,463000,9001,,,,,, +SAD69 / Brazil Polyconic,29100,4291,19941,9001,9818,4499,8801,0,9102,8802,-54,9102,8806,5000000,9001,8807,10000000,9001,,,,,,,,, +SAD69 / Brazil Polyconic,29101,4618,19941,9001,9818,4499,8801,0,9102,8802,-54,9102,8806,5000000,9001,8807,10000000,9001,,,,,,,,, SAD69 / UTM zone 18N,29118,4291,16018,9001,9807,4400,8801,0,9102,8802,-75,9102,8805,0.9996,9201,8806,500000,9001,8807,0,9001,,,,,, SAD69 / UTM zone 19N,29119,4291,16019,9001,9807,4400,8801,0,9102,8802,-69,9102,8805,0.9996,9201,8806,500000,9001,8807,0,9001,,,,,, SAD69 / UTM zone 20N,29120,4291,16020,9001,9807,4400,8801,0,9102,8802,-63,9102,8805,0.9996,9201,8806,500000,9001,8807,0,9001,,,,,, diff --git a/src/map/bsbmap.cpp b/src/map/bsbmap.cpp index c168657b..99af49fd 100644 --- a/src/map/bsbmap.cpp +++ b/src/map/bsbmap.cpp @@ -315,6 +315,9 @@ bool BSBMap::createProjection(const QString &datum, const QString &proj, } else if (!proj.compare("LAMBERT CONFORMAL CONIC", Qt::CaseInsensitive)) { Projection::Setup setup(0, params[0], NAN, 0, 0, params[2], params[3]); pcs = PCS(gcs, 9802, setup, 9001); + } else if (!proj.compare("POLYCONIC", Qt::CaseInsensitive)) { + Projection::Setup setup(0, params[0], NAN, 0, 0, NAN, NAN); + pcs = PCS(gcs, 9818, setup, 9001); } else { _errorString = proj + ": Unknown/missing projection"; return false; diff --git a/src/map/geotiff.cpp b/src/map/geotiff.cpp index 0ee12c10..7da554f1 100644 --- a/src/map/geotiff.cpp +++ b/src/map/geotiff.cpp @@ -55,6 +55,7 @@ #define CT_AlbersEqualArea 11 #define CT_PolarStereographic 15 #define CT_ObliqueStereographic 16 +#define CT_Polyconic 22 #define IS_SET(map, key) \ @@ -344,6 +345,8 @@ Projection::Method GeoTIFF::method(QMap &kv) return Projection::Method(9829); case CT_ObliqueStereographic: return Projection::Method(9809); + case CT_Polyconic: + return Projection::Method(9818); default: _errorString = QString("%1: unknown coordinate transformation method") .arg(kv.value(ProjCoordTransGeoKey).SHORT); diff --git a/src/map/mapfile.cpp b/src/map/mapfile.cpp index 39c0877c..95b3fc63 100644 --- a/src/map/mapfile.cpp +++ b/src/map/mapfile.cpp @@ -190,6 +190,8 @@ bool MapFile::createProjection(const GCS *gcs, const QString &name, pcs = PCS(gcs, 9822, setup, 9001); else if (name == "(A)Lambert Azimuthual Equal Area") pcs = PCS(gcs, 9820, setup, 9001); + else if (name == "Polyconic (American)") + pcs = PCS(gcs, 9818, setup, 9001); else if (name == "(NZTM2) New Zealand TM 2000") pcs = PCS(gcs, 9807, Projection::Setup(0, 173.0, 0.9996, 1600000, 10000000, NAN, NAN), 9001); diff --git a/src/map/polyconic.cpp b/src/map/polyconic.cpp new file mode 100644 index 00000000..09354e2c --- /dev/null +++ b/src/map/polyconic.cpp @@ -0,0 +1,211 @@ +/* + * Based on libgeotrans with the following Source Code Disclaimer: + +1. The GEOTRANS source code ("the software") is provided free of charge by +the National Imagery and Mapping Agency (NIMA) of the United States +Department of Defense. Although NIMA makes no copyright claim under Title 17 +U.S.C., NIMA claims copyrights in the source code under other legal regimes. +NIMA hereby grants to each user of the software a license to use and +distribute the software, and develop derivative works. + +2. Warranty Disclaimer: The software was developed to meet only the internal +requirements of the U.S. National Imagery and Mapping Agency. The software +is provided "as is," and no warranty, express or implied, including but not +limited to the implied warranties of merchantability and fitness for +particular purpose or arising by statute or otherwise in law or from a +course of dealing or usage in trade, is made by NIMA as to the accuracy and +functioning of the software. + +3. NIMA and its personnel are not required to provide technical support or +general assistance with respect to the software. + +4. Neither NIMA nor its personnel will be liable for any claims, losses, or +damages arising from or connected with the use of the software. The user +agrees to hold harmless the United States National Imagery and Mapping +Agency. The user's sole and exclusive remedy is to stop using the software. + +5. NIMA requests that products developed using the software credit the +source of the software with the following statement, "The product was +developed using GEOTRANS, a product of the National Imagery and Mapping +Agency and U.S. Army Engineering Research and Development Center." + +6. For any products developed using the software, NIMA requires a disclaimer +that use of the software does not indicate endorsement or approval of the +product by the Secretary of Defense or the National Imagery and Mapping +Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of +the National Imagery and Mapping Agency, the initials "NIMA", the seal of +the National Imagery and Mapping Agency, or any colorable imitation thereof +shall not be used to imply approval, endorsement, or authorization of a +product without prior written permission from United States Secretary of +Defense. + +*/ + +#include "ellipsoid.h" +#include "polyconic.h" + +#define POLY_COEFF_TIMES_SIN(coeff, x, latit) \ + (coeff * (sin (x * latit))) +#define POLY_M(c0lat, c1s2lat, c2s4lat, c3s6lat) \ + (_a * (c0lat - c1s2lat + c2s4lat - c3s6lat)) +#define FLOAT_EQ(x, v, epsilon) \ + (((v - epsilon) < x) && (x < (v + epsilon))) + +Polyconic::Polyconic(const Ellipsoid *ellipsoid, double latitudeOrigin, + double longitudeOrigin, double falseEasting, double falseNorthing) +{ + double j, three_es4; + double lat, sin2lat, sin4lat, sin6lat; + double a2; + double b2; + + _longitudeOrigin = deg2rad(longitudeOrigin); + _latitudeOrigin = deg2rad(latitudeOrigin); + _a = ellipsoid->radius(); + _b = ellipsoid->b(); + if (_longitudeOrigin > M_PI) + _longitudeOrigin -= 2 * M_PI; + _falseNorthing = falseNorthing; + _falseEasting = falseEasting; + a2 = _a * _a; + b2 = _b * _b; + _es2 = (a2 - b2) / a2; + _es4 = _es2 * _es2; + _es6 = _es4 * _es2; + + j = 45.0 * _es6 / 1024.0; + three_es4 = 3.0 * _es4; + _c0 = 1.0 - _es2 / 4.0 - three_es4 / 64.0 - 5.0 * _es6 / 256.0; + _c1 = 3.0 * _es2 / 8.0 + three_es4 / 32.0 + j; + _c2 = 15.0 * _es4 / 256.0 + j; + _c3 = 35.0 * _es6 / 3072.0; + + lat = _c0 * _latitudeOrigin; + sin2lat = POLY_COEFF_TIMES_SIN(_c1, 2.0, _latitudeOrigin); + sin4lat = POLY_COEFF_TIMES_SIN(_c2, 4.0, _latitudeOrigin); + sin6lat = POLY_COEFF_TIMES_SIN(_c3, 6.0, _latitudeOrigin); + _M0 = POLY_M(lat, sin2lat, sin4lat, sin6lat); +} + +PointD Polyconic::ll2xy(const Coordinates &c) const +{ + double Longitude = deg2rad(c.lon()); + double Latitude = deg2rad(c.lat()); + double slat = sin(Latitude); + double lat, sin2lat, sin4lat, sin6lat; + double dlam; + double NN; + double NN_OVER_tlat; + double MM; + double EE; + + dlam = Longitude - _longitudeOrigin; + if (dlam > M_PI) + dlam -= 2 * M_PI; + if (dlam < -M_PI) + dlam += 2 * M_PI; + + if (Latitude == 0.0) { + return PointD(_a * dlam + _falseEasting, + -_M0 + _falseNorthing); + } else { + NN = _a / sqrt(1.0 - _es2 * (slat * slat)); + NN_OVER_tlat = NN / tan(Latitude); + lat = _c0 * Latitude; + sin2lat = POLY_COEFF_TIMES_SIN(_c1, 2.0, Latitude); + sin4lat = POLY_COEFF_TIMES_SIN(_c2, 4.0, Latitude); + sin6lat = POLY_COEFF_TIMES_SIN(_c3, 6.0, Latitude); + MM = POLY_M(lat, sin2lat, sin4lat, sin6lat); + EE = dlam * slat; + return PointD(NN_OVER_tlat * sin(EE) + _falseEasting, + MM - _M0 + NN_OVER_tlat * (1.0 - cos(EE)) + _falseNorthing); + } +} + +Coordinates Polyconic::xy2ll(const PointD &p) const +{ + double dx; + double dy; + double dx_OVER_Poly_a; + double AA; + double BB; + double CC = 0.0; + double PHIn, Delta_PHI = 1.0; + double sin_PHIn; + double PHI, sin2PHI, sin4PHI, sin6PHI; + double Mn, Mn_prime, Ma; + double AA_Ma; + double Ma2_PLUS_BB; + double AA_MINUS_Ma; + double tolerance = 1.0e-12; + double Latitude; + double Longitude; + + + dy = p.y() - _falseNorthing; + dx = p.x() - _falseEasting; + dx_OVER_Poly_a = dx / _a; + + if (FLOAT_EQ(dy,-_M0,1)) { + Latitude = 0.0; + Longitude = dx_OVER_Poly_a + _longitudeOrigin; + } else { + AA = (_M0 + dy) / _a; + BB = dx_OVER_Poly_a * dx_OVER_Poly_a + (AA * AA); + PHIn = AA; + + while (fabs(Delta_PHI) > tolerance) { + sin_PHIn = sin(PHIn); + CC = sqrt(1.0 - _es2 * sin_PHIn * sin_PHIn) * tan(PHIn); + PHI = _c0 * PHIn; + sin2PHI = POLY_COEFF_TIMES_SIN(_c1, 2.0, PHIn); + sin4PHI = POLY_COEFF_TIMES_SIN(_c2, 4.0, PHIn); + sin6PHI = POLY_COEFF_TIMES_SIN(_c3, 6.0, PHIn); + Mn = POLY_M(PHI, sin2PHI, sin4PHI, sin6PHI); + Mn_prime = _c0 - 2.0 * _c1 * cos(2.0 * PHIn) + 4.0 * _c2 + * cos(4.0 * PHIn) - 6.0 * _c3 * cos(6.0 * PHIn); + Ma = Mn / _a; + AA_Ma = AA * Ma; + Ma2_PLUS_BB = Ma * Ma + BB; + AA_MINUS_Ma = AA - Ma; + Delta_PHI = (AA_Ma * CC + AA_MINUS_Ma - 0.5 * (Ma2_PLUS_BB) * CC) / + (_es2 * sin2PHI * (Ma2_PLUS_BB - 2.0 * AA_Ma) / 4.0 * CC + + (AA_MINUS_Ma) * (CC * Mn_prime - 2.0 / sin2PHI) - Mn_prime); + PHIn -= Delta_PHI; + } + Latitude = PHIn; + + if (Latitude > M_PI_2) + Latitude = M_PI_2; + else if (Latitude < -M_PI_2) + Latitude = -M_PI_2; + + if (FLOAT_EQ(fabs(Latitude), M_PI_2, 0.00001) || (Latitude == 0)) + Longitude = _longitudeOrigin; + else + Longitude = (asin(dx_OVER_Poly_a * CC)) / sin(Latitude) + + _longitudeOrigin; + } + + if (Longitude > M_PI) + Longitude -= 2 * M_PI; + if (Longitude < -M_PI) + Longitude += 2 *M_PI; + + if (Longitude > M_PI) + Longitude = M_PI; + else if (Longitude < -M_PI) + Longitude = -M_PI; + + return Coordinates(rad2deg(Longitude), rad2deg(Latitude)); +} + +bool Polyconic::operator==(const CT &ct) const +{ + const Polyconic *other = dynamic_cast(&ct); + return (other != 0 && _a == other->_a && _b == other->_b + && _latitudeOrigin == other->_latitudeOrigin + && _longitudeOrigin == other->_longitudeOrigin + && _falseNorthing == other->_falseNorthing + && _falseEasting == other->_falseEasting); +} diff --git a/src/map/polyconic.h b/src/map/polyconic.h new file mode 100644 index 00000000..36ae2e79 --- /dev/null +++ b/src/map/polyconic.h @@ -0,0 +1,37 @@ +#ifndef POLYCONIC_H +#define POLYCONIC_H + +#include "ct.h" + +class Ellipsoid; + +class Polyconic : public CT +{ +public: + Polyconic(const Ellipsoid *ellipsoid, double latitudeOrigin, + double longitudeOrigin, double falseEasting, double falseNorthing); + + virtual CT *clone() const {return new Polyconic(*this);} + virtual bool operator==(const CT &ct) const; + + virtual PointD ll2xy(const Coordinates &c) const; + virtual Coordinates xy2ll(const PointD &p) const; + +private: + double _a; + double _b; + double _es2; + double _es4; + double _es6; + double _M0; + double _c0; + double _c1; + double _c2; + double _c3; + double _longitudeOrigin; + double _latitudeOrigin; + double _falseEasting; + double _falseNorthing; +}; + +#endif // POLYCONIC_H diff --git a/src/map/projection.cpp b/src/map/projection.cpp index c6542ac3..e4f51cf5 100644 --- a/src/map/projection.cpp +++ b/src/map/projection.cpp @@ -8,6 +8,7 @@ #include "krovak.h" #include "polarstereographic.h" #include "obliquestereographic.h" +#include "polyconic.h" #include "latlon.h" #include "gcs.h" #include "pcs.h" @@ -25,6 +26,7 @@ Projection::Method::Method(int id) case 9807: case 9809: case 9815: + case 9818: case 9819: case 9820: case 9822: @@ -85,6 +87,11 @@ Projection::Projection(const PCS *pcs) : _gcs(0), _ct(0), _geographic(false) setup.longitudeOrigin(), setup.scale(), setup.falseEasting(), setup.falseNorthing()); break; + case 9818: + _ct = new Polyconic(ellipsoid, setup.latitudeOrigin(), + setup.longitudeOrigin(), setup.falseEasting(), + setup.falseNorthing()); + break; case 9819: _ct = new Krovak(ellipsoid, setup.standardParallel1(), setup.standardParallel2(), setup.scale(), setup.latitudeOrigin(),