Added support for Mercator projection (the real one)

gpxsee.pro

@@ -127,7 +127,8 @@ HEADERS += src/config.h \
     src/data/nmeaparser.h \
     src/data/oziparsers.h \
     src/map/rectd.h \
-    src/map/geocentric.h
+    src/map/geocentric.h \
+    src/map/mercator.h
 SOURCES += src/main.cpp \
     src/common/coordinates.cpp \
     src/common/rectc.cpp \
@@ -222,7 +223,8 @@ SOURCES += src/main.cpp \
     src/data/igcparser.cpp \
     src/data/nmeaparser.cpp \
     src/data/oziparsers.cpp \
-    src/map/geocentric.cpp
+    src/map/geocentric.cpp \
+    src/map/mercator.cpp
 RESOURCES += gpxsee.qrc
 TRANSLATIONS = lang/gpxsee_cs.ts \
     lang/gpxsee_sv.ts \

src/map/geotiff.cpp

@@ -9,49 +9,51 @@
 #define TIFF_SHORT  3
 #define TIFF_LONG   4
 
-#define ModelPixelScaleTag           33550
+#define ModelPixelScaleTag            33550
-#define ModelTiepointTag             33922
+#define ModelTiepointTag              33922
-#define ModelTransformationTag       34264
+#define ModelTransformationTag        34264
-#define GeoKeyDirectoryTag           34735
+#define GeoKeyDirectoryTag            34735
-#define GeoDoubleParamsTag           34736
+#define GeoDoubleParamsTag            34736
 
-#define GTModelTypeGeoKey            1024
+#define GTModelTypeGeoKey             1024
-#define GTRasterTypeGeoKey           1025
+#define GTRasterTypeGeoKey            1025
-#define GeographicTypeGeoKey         2048
+#define GeographicTypeGeoKey          2048
-#define GeogGeodeticDatumGeoKey      2050
+#define GeogGeodeticDatumGeoKey       2050
-#define GeogPrimeMeridianGeoKey      2051
+#define GeogPrimeMeridianGeoKey       2051
-#define GeogAngularUnitsGeoKey       2054
+#define GeogAngularUnitsGeoKey        2054
-#define GeogEllipsoidGeoKey          2056
+#define GeogEllipsoidGeoKey           2056
-#define GeogAzimuthUnitsGeoKey       2060
+#define GeogAzimuthUnitsGeoKey        2060
-#define ProjectedCSTypeGeoKey        3072
+#define ProjectedCSTypeGeoKey         3072
-#define ProjectionGeoKey             3074
+#define ProjectionGeoKey              3074
-#define ProjCoordTransGeoKey         3075
+#define ProjCoordTransGeoKey          3075
-#define ProjLinearUnitsGeoKey        3076
+#define ProjLinearUnitsGeoKey         3076
-#define ProjStdParallel1GeoKey       3078
+#define ProjStdParallel1GeoKey        3078
-#define ProjStdParallel2GeoKey       3079
+#define ProjStdParallel2GeoKey        3079
-#define ProjNatOriginLongGeoKey      3080
+#define ProjNatOriginLongGeoKey       3080
-#define ProjNatOriginLatGeoKey       3081
+#define ProjNatOriginLatGeoKey        3081
-#define ProjFalseEastingGeoKey       3082
+#define ProjFalseEastingGeoKey        3082
-#define ProjFalseNorthingGeoKey      3083
+#define ProjFalseNorthingGeoKey       3083
-#define ProjFalseOriginLatGeoKey     3085
+#define ProjFalseOriginLatGeoKey      3085
-#define ProjCenterLongGeoKey         3088
+#define ProjCenterLongGeoKey          3088
-#define ProjCenterLatGeoKey          3089
+#define ProjCenterLatGeoKey           3089
-#define ProjScaleAtNatOriginGeoKey   3092
+#define ProjCenterEastingGeoKey       3090
-#define ProjScaleAtCenterGeoKey      3093
+#define ProjFalseOriginNorthingGeoKey 3091
-#define ProjAzimuthAngleGeoKey       3094
+#define ProjScaleAtNatOriginGeoKey    3092
-#define ProjRectifiedGridAngleGeoKey 3096
+#define ProjScaleAtCenterGeoKey       3093
+#define ProjAzimuthAngleGeoKey        3094
+#define ProjRectifiedGridAngleGeoKey  3096
 
-#define ModelTypeProjected           1
+#define ModelTypeProjected            1
-#define ModelTypeGeographic          2
+#define ModelTypeGeographic           2
-#define ModelTypeGeocentric          3
+#define ModelTypeGeocentric           3
 
-#define CT_TransverseMercator        1
+#define CT_TransverseMercator         1
-#define CT_ObliqueMercator           3
+#define CT_ObliqueMercator            3
-#define CT_Mercator                  7
+#define CT_Mercator                   7
-#define CT_LambertConfConic_2SP      8
+#define CT_LambertConfConic_2SP       8
-#define CT_LambertConfConic_1SP      9
+#define CT_LambertConfConic_1SP       9
-#define CT_LambertAzimEqualArea      10
+#define CT_LambertAzimEqualArea       10
-#define CT_AlbersEqualArea           11
+#define CT_AlbersEqualArea            11
 
 
 #define IS_SET(map, key) \
@@ -253,6 +255,8 @@ bool GeoTIFF::readKeys(TIFFFile &file, Ctx &ctx, QMap<quint16, Value> &kv) const
 			case ProjAzimuthAngleGeoKey:
 			case ProjRectifiedGridAngleGeoKey:
 			case ProjFalseOriginLatGeoKey:
+			case ProjCenterEastingGeoKey:
+			case ProjFalseOriginNorthingGeoKey:
 				if (!readGeoValue(file, ctx.values, entry.ValueOffset,
 				  value.DOUBLE))
 					return false;
@@ -327,7 +331,7 @@ Projection::Method GeoTIFF::method(QMap<quint16, Value> &kv)
 		case CT_ObliqueMercator:
 			return Projection::Method(9815);
 		case CT_Mercator:
-			return Projection::Method(1024);
+			return Projection::Method(9804);
 		case CT_LambertConfConic_2SP:
 			return Projection::Method(9802);
 		case CT_LambertConfConic_1SP:
@@ -421,10 +425,14 @@ bool GeoTIFF::projectedModel(QMap<quint16, Value> &kv)
 			sp2 = NAN;
 		if (kv.contains(ProjFalseNorthingGeoKey))
 			fn = lu.toMeters(kv.value(ProjFalseNorthingGeoKey).DOUBLE);
+		else if (kv.contains(ProjFalseOriginNorthingGeoKey))
+			fn = lu.toMeters(kv.value(ProjFalseOriginNorthingGeoKey).DOUBLE);
 		else
 			fn = NAN;
 		if (kv.contains(ProjFalseEastingGeoKey))
 			fe = lu.toMeters(kv.value(ProjFalseEastingGeoKey).DOUBLE);
+		else if (kv.contains(ProjCenterEastingGeoKey))
+			fe = lu.toMeters(kv.value(ProjCenterEastingGeoKey).DOUBLE);
 		else
 			fe = NAN;
 

src/map/mercator.cpp

@@ -0,0 +1,122 @@
+/*
+ * 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 "mercator.h"
+
+Mercator::Mercator(const Ellipsoid *ellipsoid, double latitudeOrigin,
+  double longitudeOrigin, double falseEasting, double falseNorthing)
+{
+	double es2;
+	double es3;
+	double es4;
+	double sin_olat;
+
+	_a = ellipsoid->radius();
+	_latitudeOrigin = deg2rad(latitudeOrigin);
+	_longitudeOrigin = deg2rad(longitudeOrigin);
+	if (_longitudeOrigin > M_PI)
+		_longitudeOrigin -= 2*M_PI;
+	_falseNorthing = falseNorthing;
+	_falseEasting = falseEasting;
+	_es = 2 * ellipsoid->flattening() - ellipsoid->flattening()
+	  * ellipsoid->flattening();
+	_e = sqrt(_es);
+	sin_olat = sin(_latitudeOrigin);
+	_scaleFactor = 1.0 / (sqrt(1.e0 - _es * sin_olat * sin_olat)
+	  / cos(_latitudeOrigin));
+	es2 = _es * _es;
+	es3 = es2 * _es;
+	es4 = es3 * _es;
+	_ab = _es / 2.e0 + 5.e0 * es2 / 24.e0 + es3 / 12.e0 + 13.e0 * es4 / 360.e0;
+	_bb = 7.e0 * es2 / 48.e0 + 29.e0 * es3 / 240.e0 + 811.e0 * es4 / 11520.e0;
+	_cb = 7.e0 * es3 / 120.e0 + 81.e0 * es4 / 1120.e0;
+	_db = 4279.e0 * es4 / 161280.e0;
+}
+
+PointD Mercator::ll2xy(const Coordinates &c) const
+{
+	double lon = deg2rad(c.lon());
+	double lat = deg2rad(c.lat());
+	double ctanz2;
+	double e_x_sinlat;
+	double delta_lon;
+	double tan_temp;
+	double pow_temp;
+
+	if (lon > M_PI)
+		lon -= 2*M_PI;
+	e_x_sinlat = _e * sin(lat);
+	tan_temp = tan(M_PI / 4.e0 + lat / 2.e0);
+	pow_temp = pow((1.e0 - e_x_sinlat) / (1.e0 + e_x_sinlat), _e / 2.e0);
+	ctanz2 = tan_temp * pow_temp;
+	delta_lon = lon - _longitudeOrigin;
+	if (delta_lon > M_PI)
+	  delta_lon -= 2 * M_PI;
+	if (delta_lon < -M_PI)
+	  delta_lon += 2 * M_PI;
+
+	return PointD(_scaleFactor * _a * delta_lon + _falseEasting,
+	  _scaleFactor * _a * log(ctanz2) + _falseNorthing);
+}
+
+Coordinates Mercator::xy2ll(const PointD &p) const
+{
+	double dx;
+	double dy;
+	double xphi;
+	double lat, lon;
+
+	dy = p.y() - _falseNorthing;
+	dx = p.x() - _falseEasting;
+	lon = _longitudeOrigin + dx / (_scaleFactor * _a);
+	xphi = M_PI / 2.e0 - 2.e0 * atan(1.e0 / exp(dy / (_scaleFactor * _a)));
+	lat = xphi + _ab * sin(2.e0 * xphi) + _bb * sin(4.e0 * xphi)
+	  + _cb * sin(6.e0 * xphi) + _db * sin(8.e0 * xphi);
+	if (lon > M_PI)
+		lon -= 2 * M_PI;
+	if (lon < -M_PI)
+		lon += 2 * M_PI;
+
+	return Coordinates(rad2deg(lon), rad2deg(lat));
+}

src/map/mercator.h

@@ -0,0 +1,29 @@
+#ifndef MERCATOR_H
+#define MERCATOR_H
+
+#include "ct.h"
+
+class Ellipsoid;
+
+class Mercator : public CT
+{
+public:
+	Mercator(const Ellipsoid *ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double falseEasting, double falseNorthing);
+
+	virtual CT *clone() const {return new Mercator(*this);}
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _a, _es, _e;
+	double _latitudeOrigin;
+	double _longitudeOrigin;
+	double _falseNorthing;
+	double _falseEasting;
+	double _scaleFactor;
+	double _ab, _bb, _cb, _db;
+};
+
+#endif // MERCATOR_H

src/map/projection.cpp

@@ -1,4 +1,5 @@
 #include "datum.h"
+#include "mercator.h"
 #include "webmercator.h"
 #include "transversemercator.h"
 #include "lambertconic.h"
@@ -16,11 +17,11 @@ Projection::Method::Method(int id)
 		case 1024:
 		case 9801:
 		case 9802:
+		case 9804:
 		case 9807:
 		case 9815:
 		case 9820:
 		case 9822:
-		case 9841:
 			_id = id;
 			break;
 		default:
@@ -36,7 +37,6 @@ Projection::Projection(const PCS *pcs) : _gcs(pcs->gcs()), _units(pcs->units()),
 
 	switch (pcs->method().id()) {
 		case 1024:
-		case 9841:
 			_ct = new WebMercator();
 			break;
 		case 9801:
@@ -51,6 +51,11 @@ Projection::Projection(const PCS *pcs) : _gcs(pcs->gcs()), _units(pcs->units()),
 			  setup.longitudeOrigin(), setup.falseEasting(),
 			  setup.falseNorthing());
 			break;
+		case 9804:
+			_ct = new Mercator(ellipsoid, setup.latitudeOrigin(),
+			  setup.longitudeOrigin(), setup.falseEasting(),
+			  setup.falseNorthing());
+			break;
 		case 9807:
 			_ct = new TransverseMercator(ellipsoid, setup.latitudeOrigin(),
 			  setup.longitudeOrigin(), setup.scale(), setup.falseEasting(),