Added support for polyconic projections

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 \

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,,,,,,

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;

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<quint16, Value> &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);

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);

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<const Polyconic*>(&ct);
+	return (other != 0 && _a == other->_a && _b == other->_b
+	  && _latitudeOrigin == other->_latitudeOrigin
+	  && _longitudeOrigin == other->_longitudeOrigin
+	  && _falseNorthing == other->_falseNorthing
+	  && _falseEasting == other->_falseEasting);
+}

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

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(),