Added support for polar stereographic projection

closes #181

gpxsee.pro

@@ -148,7 +148,8 @@ HEADERS += src/common/config.h \
     src/data/oziparsers.h \
     src/data/locparser.h \
     src/data/slfparser.h \
-    src/data/dem.h
+    src/data/dem.h \
+    src/map/polarstereographic.h
 SOURCES += src/main.cpp \
     src/common/coordinates.cpp \
     src/common/rectc.cpp \
@@ -257,7 +258,9 @@ SOURCES += src/main.cpp \
     src/data/oziparsers.cpp \
     src/data/locparser.cpp \
     src/data/slfparser.cpp \
-    src/data/dem.cpp
+    src/data/dem.cpp \
+    src/map/polarstereographic.cpp \
+    src/map/rectd.cpp
 RESOURCES += gpxsee.qrc
 TRANSLATIONS = lang/gpxsee_en.ts \
     lang/gpxsee_cs.ts \

pkg/csv/pcs.csv

@@ -859,6 +859,7 @@ RT90 2.5 gon V,3021,4124,19929,9001,9807,4530,8801,0,9110,8802,15.48298,9110,880
 RT90 0 gon,3022,4124,17336,9001,9807,4530,8801,0,9110,8802,18.03298,9110,8805,1,9201,8806,1500000,9001,8807,0,9001,,,,,,
 RT90 2.5 gon O,3023,4124,17337,9001,9807,4530,8801,0,9110,8802,20.18298,9110,8805,1,9201,8806,1500000,9001,8807,0,9001,,,,,,
 RT90 5 gon O,3024,4124,17338,9001,9807,4530,8801,0,9110,8802,22.33298,9110,8805,1,9201,8806,1500000,9001,8807,0,9001,,,,,,
+WGS 84 / Antarctic Polar Stereographic,3031,4326,19992,9001,9829,4490,8806,0,9001,8807,0,9001,8832,-71,9102,8833,0,9102,,,,,,,,,
 WGS 84 / Australian Antarctic Lambert,3033,4326,19994,9001,9802,4400,8821,-50,9110,8822,70,9110,8823,-68.3,9110,8824,-74.3,9110,8826,6000000,9001,8827,6000000,9001,,,
 ETRS89 / LCC Europe,3034,4258,19985,9001,9802,4500,8821,52,9102,8822,10,9102,8823,35,9102,8824,65,9102,8826,4000000,9001,8827,2800000,9001,,,
 ETRS89 / LAEA Europe,3035,4258,19986,9001,9820,4532,8801,52,9102,8802,10,9102,8806,4321000,9001,8807,3210000,9001,,,,,,,,,

src/map/coordinatesystem.cpp

@@ -18,6 +18,7 @@ CoordinateSystem::CoordinateSystem(int code)
 		case 4467:
 		case 4469:
 		case 4470:
+		case 4490:
 		case 4495:
 		case 4496:
 		case 4497:

src/map/geotiff.cpp

@@ -42,6 +42,7 @@
 #define ProjScaleAtNatOriginGeoKey     3092
 #define ProjScaleAtCenterGeoKey        3093
 #define ProjAzimuthAngleGeoKey         3094
+#define ProjStraightVertPoleLongGeoKey 3095
 #define ProjRectifiedGridAngleGeoKey   3096
 
 #define ModelTypeProjected             1
@@ -55,6 +56,7 @@
 #define CT_LambertConfConic_1SP        9
 #define CT_LambertAzimEqualArea        10
 #define CT_AlbersEqualArea             11
+#define CT_PolarStereographic          15
 
 
 #define IS_SET(map, key) \
@@ -261,6 +263,7 @@ bool GeoTIFF::readKeys(TIFFFile &file, Ctx &ctx, QMap<quint16, Value> &kv) const
 			case ProjCenterNorthingGeoKey:
 			case ProjFalseOriginEastingGeoKey:
 			case ProjFalseOriginNorthingGeoKey:
+			case ProjStraightVertPoleLongGeoKey:
 				if (!readGeoValue(file, ctx.values, entry.ValueOffset,
 				  value.DOUBLE))
 					return false;
@@ -342,6 +345,8 @@ Projection::Method GeoTIFF::method(QMap<quint16, Value> &kv)
 			return Projection::Method(9820);
 		case CT_AlbersEqualArea:
 			return Projection::Method(9822);
+		case CT_PolarStereographic:
+			return Projection::Method(9829);
 		default:
 			_errorString = QString("%1: unknown coordinate transformation method")
 			  .arg(kv.value(ProjCoordTransGeoKey).SHORT);
@@ -407,6 +412,8 @@ bool GeoTIFF::projectedModel(QMap<quint16, Value> &kv)
 			lon0 = au.toDegrees(kv.value(ProjCenterLongGeoKey).DOUBLE);
 		else if (kv.contains(ProjFalseOriginLongGeoKey))
 			lon0 = au.toDegrees(kv.value(ProjFalseOriginLongGeoKey).DOUBLE);
+		else if (kv.contains(ProjStraightVertPoleLongGeoKey))
+			lon0 = au.toDegrees(kv.value(ProjStraightVertPoleLongGeoKey).DOUBLE);
 		else
 			lon0 = NAN;
 		if (kv.contains(ProjScaleAtNatOriginGeoKey))

src/map/pcs.cpp

@@ -32,6 +32,7 @@ static bool parameter(int key, double val, int units, Projection::Setup &setup)
 		case 8801:
 		case 8811:
 		case 8821:
+		case 8832:
 			{AngularUnits au(units);
 			if (au.isNull())
 				return false;

src/map/polarstereographic.cpp

@@ -0,0 +1,192 @@
+/*
+ * 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 "polarstereographic.h"
+
+
+#define POLAR_POW(EsSin) pow((1.0 - EsSin) / (1.0 + EsSin), _es_OVER_2)
+
+PolarStereographic::PolarStereographic(const Ellipsoid *ellipsoid,
+  double latitudeOrigin, double longitudeOrigin,
+  double falseEasting, double falseNorthing)
+{
+	_two_a = 2.0 * ellipsoid->radius();
+
+	if (longitudeOrigin > M_PI)
+		longitudeOrigin -= 2*M_PI;
+
+	if (latitudeOrigin < 0) {
+		_southernHemisphere = 1;
+		_originLatitude = -deg2rad(latitudeOrigin);
+		_originLongitude = -deg2rad(longitudeOrigin);
+	} else {
+		_southernHemisphere = 0;
+		_originLatitude = deg2rad(latitudeOrigin);
+		_originLongitude = deg2rad(longitudeOrigin);
+	}
+	_falseEasting = falseEasting;
+	_falseNorthing = falseNorthing;
+
+	double es2 = ellipsoid->es();
+	_es = sqrt(es2);
+	_es_OVER_2 = _es / 2.0;
+
+	if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10) {
+		double slat = sin(_originLatitude);
+		double essin = _es * slat;
+		double pow_es = POLAR_POW(essin);
+		double clat = cos(_originLatitude);
+		_mc = clat / sqrt(1.0 - essin * essin);
+		_a_mc = ellipsoid->radius() * _mc;
+		_tc = tan(M_PI_4 - _originLatitude / 2.0) / pow_es;
+	} else {
+		double one_PLUS_es = 1.0 + _es;
+		double one_MINUS_es = 1.0 - _es;
+		_e4 = sqrt(pow(one_PLUS_es, one_PLUS_es) * pow(one_MINUS_es,
+		  one_MINUS_es));
+	}
+}
+
+PointD PolarStereographic::ll2xy(const Coordinates &c) const
+{
+	double Easting, Northing;
+	double Longitude = deg2rad(c.lon());
+	double Latitude = deg2rad(c.lat());
+
+
+	if (fabs(fabs(Latitude) - M_PI_2) < 1.0e-10) {
+		Easting = 0.0;
+		Northing = 0.0;
+	} else {
+		if (_southernHemisphere != 0) {
+			Longitude *= -1.0;
+			Latitude *= -1.0;
+		}
+
+		double dlam = Longitude - _originLongitude;
+		if (dlam > M_PI)
+			dlam -= 2*M_PI;
+		if (dlam < -M_PI)
+			dlam += 2*M_PI;
+
+		double slat = sin(Latitude);
+		double essin = _es * slat;
+		double pow_es = POLAR_POW(essin);
+		double t = tan(M_PI_4 - Latitude / 2.0) / pow_es;
+
+		double rho;
+		if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10)
+			rho = _a_mc * t / _tc;
+		else
+			rho = _two_a * t / _e4;
+
+		Easting = rho * sin(dlam) + _falseEasting;
+
+		if (_southernHemisphere != 0) {
+			Easting *= -1.0;
+			Northing = rho * cos(dlam) + _falseNorthing;
+		} else
+			Northing = -rho * cos(dlam) + _falseNorthing;
+	}
+
+	return PointD(Easting, Northing);
+}
+
+Coordinates PolarStereographic::xy2ll(const PointD &p) const
+{
+	double Latitude, Longitude;
+	double dy = p.y() - _falseNorthing;
+	double dx = p.x() - _falseEasting;
+
+	if ((dy == 0.0) && (dx == 0.0)) {
+		Latitude = M_PI_2;
+		Longitude = _originLongitude;
+	} else {
+		if (_southernHemisphere != 0) {
+			dy *= -1.0;
+			dx *= -1.0;
+		}
+
+		double rho = sqrt(dx * dx + dy * dy);
+		double t;
+		if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10)
+			t = rho * _tc / (_a_mc);
+		else
+			t = rho * _e4 / (_two_a);
+		double PHI = M_PI_2 - 2.0 * atan(t);
+
+		double tempPHI = 0.0;
+		while (fabs(PHI - tempPHI) > 1.0e-10) {
+			tempPHI = PHI;
+			double sin_PHI = sin(PHI);
+			double essin =  _es * sin_PHI;
+			double pow_es = POLAR_POW(essin);
+			PHI = M_PI_2 - 2.0 * atan(t * pow_es);
+		}
+		Latitude = PHI;
+		Longitude = _originLongitude + atan2(dx, -dy);
+
+		if (Longitude > M_PI)
+			Longitude -= 2*M_PI;
+		else if (Longitude < -M_PI)
+			Longitude += 2*M_PI;
+
+		if (Latitude > M_PI_2)
+			Latitude = M_PI_2;
+		else if (Latitude < -M_PI_2)
+			Latitude = -M_PI_2;
+
+		if (Longitude > M_PI)
+			Longitude = M_PI;
+		else if (Longitude < -M_PI)
+			Longitude = -M_PI;
+	}
+
+	if (_southernHemisphere != 0) {
+		Latitude *= -1.0;
+		Longitude *= -1.0;
+	}
+
+	return Coordinates(rad2deg(Longitude), rad2deg(Latitude));
+}

src/map/polarstereographic.h

@@ -0,0 +1,34 @@
+#ifndef POLARSTEREOGRAPHIC_H
+#define POLARSTEREOGRAPHIC_H
+
+#include "ct.h"
+
+class Ellipsoid;
+
+class PolarStereographic  : public CT
+{
+public:
+	PolarStereographic(const Ellipsoid *ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double falseEasting, double falseNorthing);
+	virtual CT *clone() const {return new PolarStereographic(*this);}
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _originLatitude;
+	double _originLongitude;
+	double _falseEasting;
+	double _falseNorthing;
+
+	double _a_mc;
+	double _es;
+	double _es_OVER_2;
+	double _two_a;
+	double _mc;
+	double _tc;
+	double _e4;
+	int _southernHemisphere;
+};
+
+#endif // POLARSTEREOGRAPHIC_H

src/map/projection.cpp

@@ -6,6 +6,7 @@
 #include "albersequal.h"
 #include "lambertazimuthal.h"
 #include "krovak.h"
+#include "polarstereographic.h"
 #include "latlon.h"
 #include "gcs.h"
 #include "pcs.h"
@@ -25,6 +26,7 @@ Projection::Method::Method(int id)
 		case 9819:
 		case 9820:
 		case 9822:
+		case 9829:
 			_id = id;
 			break;
 		default:
@@ -87,6 +89,11 @@ Projection::Projection(const PCS *pcs) : _gcs(pcs->gcs()), _units(pcs->units()),
 			  setup.longitudeOrigin(), setup.falseEasting(),
 			  setup.falseNorthing());
 			break;
+		case 9829:
+			_ct = new PolarStereographic(ellipsoid, setup.latitudeOrigin(),
+			  setup.longitudeOrigin(), setup.falseEasting(),
+			  setup.falseNorthing());
+			break;
 		default:
 			_ct = 0;
 	}

src/map/rectd.cpp

@@ -0,0 +1,53 @@
+#include "common/rectc.h"
+#include "projection.h"
+#include "rectd.h"
+
+
+#define SAMPLE_POINTS 100
+
+static void growRect(const Projection &proj, const Coordinates &c, RectD &rect)
+{
+	if (c.isNull())
+		return;
+
+	PointD p(proj.ll2xy(c));
+
+	if (rect.isNull())
+		rect = RectD(p, p);
+	else {
+		if (p.x() < rect.left())
+			rect.setLeft(p.x());
+		if (p.x() > rect.right())
+			rect.setRight(p.x());
+		if (p.y() < rect.bottom())
+			rect.setBottom(p.y());
+		if (p.y() > rect.top())
+			rect.setTop(p.y());
+	}
+}
+
+RectD::RectD(const RectC &rect, const Projection &proj)
+{
+	RectD prect;
+	double dx = (rect.right() - rect.left()) / SAMPLE_POINTS;
+	double dy = (rect.top() - rect.bottom()) / SAMPLE_POINTS;
+
+	growRect(proj, rect.topLeft(), prect);
+
+	if (dx > 0) {
+		for (int i = 0; i <= SAMPLE_POINTS; i++) {
+			double x = rect.left() + i * dx;
+			growRect(proj, Coordinates(x, rect.bottom()), prect);
+			growRect(proj, Coordinates(x, rect.top()), prect);
+		}
+	}
+	if (dy > 0) {
+		for (int i = 0; i <= SAMPLE_POINTS; i++ ) {
+			double y = rect.bottom() + i * dy;
+			growRect(proj, Coordinates(rect.left(), y), prect);
+			growRect(proj, Coordinates(rect.right(), y), prect);
+		}
+	}
+
+	*this = prect;
+}

src/map/rectd.h

@@ -3,12 +3,16 @@
 
 #include "pointd.h"
 
+class RectC;
+class Projection;
+
 class RectD
 {
 public:
 	RectD() {}
 	RectD(const PointD &topLeft, const PointD &bottomRight)
 	  : _tl(topLeft), _br(bottomRight) {}
+	RectD(const RectC &rect, const Projection &proj);
 
 	PointD topLeft() const {return _tl;}
 	PointD bottomRight() const {return _br;}
@@ -18,6 +22,11 @@ public:
 	double top() const {return _tl.y();}
 	double bottom() const {return _br.y();}
 
+	void setLeft(double val) {_tl.rx() = val;}
+	void setRight(double val) {_br.rx() = val;}
+	void setTop(double val) {_tl.ry() = val;}
+	void setBottom(double val) {_br.ry() = val;}
+
 	double width() const {return (right() - left());}
 	double height() const {return (top() - bottom());}
 
@@ -26,6 +35,7 @@ public:
 	  && p.y() >= bottom());}
 
 	bool isNull() const {return _tl.isNull() && _br.isNull();}
+	bool isValid() const {return !(_tl.isNull() || _br.isNull());}
 
 private:
 	PointD _tl, _br;

src/map/wmsmap.cpp

@@ -65,12 +65,8 @@ void WMSMap::updateTransform()
 	double pixelSpan = sd2res(_zooms.at(_zoom));
 	if (_projection.isGeographic())
 		pixelSpan /= deg2rad(WGS84_RADIUS);
-	double sx = _bbox.width() / pixelSpan;
+	_transform = Transform(ReferencePoint(PointD(0, 0),
-	double sy = _bbox.height() / pixelSpan;
+	  _projection.ll2xy(_bbox.topLeft())), PointD(pixelSpan, pixelSpan));
-
-	ReferencePoint tl(PointD(0, 0), _bbox.topLeft());
-	ReferencePoint br(PointD(sx, sy), _bbox.bottomRight());
-	_transform = Transform(tl, br);
 }
 
 bool WMSMap::loadWMS()
@@ -84,8 +80,8 @@ bool WMSMap::loadWMS()
 	}
 
 	_projection = wms.projection();
-	_bbox = RectD(_projection.ll2xy(wms.boundingBox().topLeft()),
+	_bbox = wms.boundingBox();
-	  _projection.ll2xy(wms.boundingBox().bottomRight()));
+	_bounds = RectD(_bbox, _projection);
 	_tileLoader->setUrl(tileUrl(wms.version()));
 
 	if (wms.version() >= "1.3.0") {
@@ -124,8 +120,8 @@ void WMSMap::clearCache()
 
 QRectF WMSMap::bounds()
 {
-	return QRectF(_transform.proj2img(_bbox.topLeft()) / _mapRatio,
+	return QRectF(_transform.proj2img(_bounds.topLeft()) / _mapRatio,
-	  _transform.proj2img(_bbox.bottomRight()) / _mapRatio);
+	  _transform.proj2img(_bounds.bottomRight()) / _mapRatio);
 }
 
 int WMSMap::zoomFit(const QSize &size, const RectC &rect)

src/map/wmsmap.h

@@ -55,7 +55,8 @@ private:
 	Transform _transform;
 	CoordinateSystem _cs;
 	QVector<double> _zooms;
-	RectD _bbox;
+	RectC _bbox;
+	RectD _bounds;
 	int _zoom;
 	qreal _mapRatio;
 

src/map/wmtsmap.cpp

@@ -22,10 +22,10 @@ bool WMTSMap::loadWMTS()
 		return false;
 	}
 
-	_bounds = wmts.bounds();
 	_zooms = wmts.zooms();
 	_projection = wmts.projection();
 	_tileLoader->setUrl(wmts.tileUrl());
+	_bounds = RectD(wmts.bounds(), _projection);
 
 	if (_setup.coordinateSystem().axisOrder() == CoordinateSystem::Unknown)
 		_cs = _projection.coordinateSystem();
@@ -77,14 +77,8 @@ void WMTSMap::updateTransform()
 	double pixelSpan = sd2res(z.scaleDenominator());
 	if (_projection.isGeographic())
 		pixelSpan /= deg2rad(WGS84_RADIUS);
-	PointD tileSpan(z.tile().width() * pixelSpan, z.tile().height() * pixelSpan);
+	_transform = Transform(ReferencePoint(PointD(0, 0), topLeft),
-	PointD bottomRight(topLeft.x() + tileSpan.x() * z.matrix().width(),
+	  PointD(pixelSpan, pixelSpan));
-	  topLeft.y() - tileSpan.y() * z.matrix().height());
-
-	ReferencePoint tl(PointD(0, 0), topLeft);
-	ReferencePoint br(PointD(z.tile().width() * z.matrix().width(),
-	  z.tile().height() * z.matrix().height()), bottomRight);
-	_transform = Transform(tl, br);
 }
 
 QRectF WMTSMap::bounds()
@@ -99,9 +93,11 @@ QRectF WMTSMap::bounds()
 	  * z.tile().height()), QSize(z.tile().width() * z.limits().width(),
 	  z.tile().height() * z.limits().height()));
 
-	bounds = _bounds.isValid() ? QRectF(ll2xy(_bounds.topLeft()),
+	if (_bounds.isValid())
-	  ll2xy(_bounds.bottomRight())) : QRectF();
+		bounds = QRectF(_transform.proj2img(_bounds.topLeft())
-	return _bounds.isValid() ? tileBounds.intersected(bounds) : tileBounds;
+		  / coordinatesRatio(), _transform.proj2img(_bounds.bottomRight())
+		  / coordinatesRatio());
+	return bounds.isValid() ? tileBounds.intersected(bounds) : tileBounds;
 }
 
 int WMTSMap::zoomFit(const QSize &size, const RectC &rect)

src/map/wmtsmap.h

@@ -4,6 +4,7 @@
 #include "transform.h"
 #include "projection.h"
 #include "map.h"
+#include "rectd.h"
 #include "wmts.h"
 
 class TileLoader;
@@ -50,7 +51,7 @@ private:
 	QString _name;
 	WMTS::Setup _setup;
 	TileLoader *_tileLoader;
-	RectC _bounds;
+	RectD _bounds;
 	QList<WMTS::Zoom> _zooms;
 	Projection _projection;
 	Transform _transform;