From be1c7fa4c20b4bff11cac77ee4826bf54c988b87 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20T=C5=AFma?= Date: Sun, 13 May 2018 00:40:03 +0200 Subject: [PATCH] Added support for Mercator projection (the real one) --- gpxsee.pro | 6 +- src/map/geotiff.cpp | 90 ++++++++++++++++-------------- src/map/mercator.cpp | 122 +++++++++++++++++++++++++++++++++++++++++ src/map/mercator.h | 29 ++++++++++ src/map/projection.cpp | 9 ++- 5 files changed, 211 insertions(+), 45 deletions(-) create mode 100644 src/map/mercator.cpp create mode 100644 src/map/mercator.h diff --git a/gpxsee.pro b/gpxsee.pro index 168c3c79..10012377 100644 --- a/gpxsee.pro +++ b/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 \ diff --git a/src/map/geotiff.cpp b/src/map/geotiff.cpp index 4e3ab352..f33de18f 100644 --- a/src/map/geotiff.cpp +++ b/src/map/geotiff.cpp @@ -9,49 +9,51 @@ #define TIFF_SHORT 3 #define TIFF_LONG 4 -#define ModelPixelScaleTag 33550 -#define ModelTiepointTag 33922 -#define ModelTransformationTag 34264 -#define GeoKeyDirectoryTag 34735 -#define GeoDoubleParamsTag 34736 +#define ModelPixelScaleTag 33550 +#define ModelTiepointTag 33922 +#define ModelTransformationTag 34264 +#define GeoKeyDirectoryTag 34735 +#define GeoDoubleParamsTag 34736 -#define GTModelTypeGeoKey 1024 -#define GTRasterTypeGeoKey 1025 -#define GeographicTypeGeoKey 2048 -#define GeogGeodeticDatumGeoKey 2050 -#define GeogPrimeMeridianGeoKey 2051 -#define GeogAngularUnitsGeoKey 2054 -#define GeogEllipsoidGeoKey 2056 -#define GeogAzimuthUnitsGeoKey 2060 -#define ProjectedCSTypeGeoKey 3072 -#define ProjectionGeoKey 3074 -#define ProjCoordTransGeoKey 3075 -#define ProjLinearUnitsGeoKey 3076 -#define ProjStdParallel1GeoKey 3078 -#define ProjStdParallel2GeoKey 3079 -#define ProjNatOriginLongGeoKey 3080 -#define ProjNatOriginLatGeoKey 3081 -#define ProjFalseEastingGeoKey 3082 -#define ProjFalseNorthingGeoKey 3083 -#define ProjFalseOriginLatGeoKey 3085 -#define ProjCenterLongGeoKey 3088 -#define ProjCenterLatGeoKey 3089 -#define ProjScaleAtNatOriginGeoKey 3092 -#define ProjScaleAtCenterGeoKey 3093 -#define ProjAzimuthAngleGeoKey 3094 -#define ProjRectifiedGridAngleGeoKey 3096 +#define GTModelTypeGeoKey 1024 +#define GTRasterTypeGeoKey 1025 +#define GeographicTypeGeoKey 2048 +#define GeogGeodeticDatumGeoKey 2050 +#define GeogPrimeMeridianGeoKey 2051 +#define GeogAngularUnitsGeoKey 2054 +#define GeogEllipsoidGeoKey 2056 +#define GeogAzimuthUnitsGeoKey 2060 +#define ProjectedCSTypeGeoKey 3072 +#define ProjectionGeoKey 3074 +#define ProjCoordTransGeoKey 3075 +#define ProjLinearUnitsGeoKey 3076 +#define ProjStdParallel1GeoKey 3078 +#define ProjStdParallel2GeoKey 3079 +#define ProjNatOriginLongGeoKey 3080 +#define ProjNatOriginLatGeoKey 3081 +#define ProjFalseEastingGeoKey 3082 +#define ProjFalseNorthingGeoKey 3083 +#define ProjFalseOriginLatGeoKey 3085 +#define ProjCenterLongGeoKey 3088 +#define ProjCenterLatGeoKey 3089 +#define ProjCenterEastingGeoKey 3090 +#define ProjFalseOriginNorthingGeoKey 3091 +#define ProjScaleAtNatOriginGeoKey 3092 +#define ProjScaleAtCenterGeoKey 3093 +#define ProjAzimuthAngleGeoKey 3094 +#define ProjRectifiedGridAngleGeoKey 3096 -#define ModelTypeProjected 1 -#define ModelTypeGeographic 2 -#define ModelTypeGeocentric 3 +#define ModelTypeProjected 1 +#define ModelTypeGeographic 2 +#define ModelTypeGeocentric 3 -#define CT_TransverseMercator 1 -#define CT_ObliqueMercator 3 -#define CT_Mercator 7 -#define CT_LambertConfConic_2SP 8 -#define CT_LambertConfConic_1SP 9 -#define CT_LambertAzimEqualArea 10 -#define CT_AlbersEqualArea 11 +#define CT_TransverseMercator 1 +#define CT_ObliqueMercator 3 +#define CT_Mercator 7 +#define CT_LambertConfConic_2SP 8 +#define CT_LambertConfConic_1SP 9 +#define CT_LambertAzimEqualArea 10 +#define CT_AlbersEqualArea 11 #define IS_SET(map, key) \ @@ -253,6 +255,8 @@ bool GeoTIFF::readKeys(TIFFFile &file, Ctx &ctx, QMap &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 &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 &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; diff --git a/src/map/mercator.cpp b/src/map/mercator.cpp new file mode 100644 index 00000000..f26e5603 --- /dev/null +++ b/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)); +} diff --git a/src/map/mercator.h b/src/map/mercator.h new file mode 100644 index 00000000..d55db268 --- /dev/null +++ b/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 diff --git a/src/map/projection.cpp b/src/map/projection.cpp index c60b58f9..a40487d4 100644 --- a/src/map/projection.cpp +++ b/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(),