Added support for oblique stereographic projections

src/map/geotiff.cpp

@@ -57,6 +57,7 @@
 #define CT_LambertAzimEqualArea        10
 #define CT_AlbersEqualArea             11
 #define CT_PolarStereographic          15
+#define CT_ObliqueStereographic        16
 
 
 #define IS_SET(map, key) \
@@ -347,6 +348,8 @@ Projection::Method GeoTIFF::method(QMap<quint16, Value> &kv)
 			return Projection::Method(9822);
 		case CT_PolarStereographic:
 			return Projection::Method(9829);
+		case CT_ObliqueStereographic:
+			return Projection::Method(9809);
 		default:
 			_errorString = QString("%1: unknown coordinate transformation method")
 			  .arg(kv.value(ProjCoordTransGeoKey).SHORT);

src/map/obliquestereographic.cpp

@@ -0,0 +1,82 @@
+#include "ellipsoid.h"
+#include "obliquestereographic.h"
+
+
+#define S1(x) ((1.0 + x) / (1.0 - x))
+#define S2(x) ((1.0 - _e * x) / (1.0 + _e * x))
+
+ObliqueStereographic::ObliqueStereographic(const Ellipsoid *ellipsoid,
+  double latitudeOrigin, double longitudeOrigin, double scale,
+  double falseEasting, double falseNorthing)
+  : _FE(falseEasting), _FN(falseNorthing)
+{
+	double lat0 = deg2rad(latitudeOrigin);
+	double sinPhi0 = sin(lat0);
+	double cosPhi0 = cos(lat0);
+	double sinLat0s = sinPhi0 * sinPhi0;
+
+	double a = ellipsoid->radius();
+	_es = ellipsoid->es();
+	_e = sqrt(_es);
+
+	double rho0 = a * (1.0 - _es) / pow(1.0 - _es * sinLat0s, 1.5);
+	double nu0 = a / sqrt(1.0 - _es * sinLat0s);
+
+	_n = sqrt(1.0 + ((_es * pow(cosPhi0, 4)) / (1.0 - _es)));
+	double w1 = pow(S1(sinPhi0) * pow(S2(sinPhi0), _e), _n);
+	double sinChi0 = (w1 - 1) / (w1 + 1);
+	_c = (_n + sinPhi0) * (1.0 - sinChi0) / ((_n - sinPhi0) * (1.0 + sinChi0));
+	double w2 = _c * w1;
+
+	_chi0 = asin((w2 - 1.0)/(w2 + 1.0));
+	_sinChi0 = sin(_chi0); _cosChi0 = cos(_chi0);
+	_lambda0 = deg2rad(longitudeOrigin);
+	_twoRk0 = 2.0 * sqrt(rho0 * nu0) * scale;
+	_g = _twoRk0 * tan(M_PI_4 - _chi0 / 2.0);
+	_h = 2.0 * _twoRk0 * tan(_chi0) + _g;
+}
+
+PointD ObliqueStereographic::ll2xy(const Coordinates &c) const
+{
+	double lon = deg2rad(c.lon());
+	double lat = deg2rad(c.lat());
+
+	double sinPhi = sin(lat);
+	double w = _c * pow(S1(sinPhi) * pow(S2(sinPhi), _e), _n);
+	double chi = asin((w - 1.0) / (w + 1.0));
+	double lambda = _n * (lon - _lambda0) + _lambda0;
+	double B = (1.0 + sin(chi) * _sinChi0 + cos(chi) * _cosChi0
+	  * cos(lambda - _lambda0));
+
+	return PointD(_FE + _twoRk0 * cos(chi) * sin(lambda - _lambda0) / B,
+	  _FN + _twoRk0 * (sin(chi) * _cosChi0 - cos(chi) * _sinChi0
+	  * cos(lambda - _lambda0)) / B);
+}
+
+Coordinates ObliqueStereographic::xy2ll(const PointD &p) const
+{
+	double i = atan((p.x() - _FE) / (_h + (p.y() - _FN)));
+	double j = atan((p.x() - _FE) / (_g - (p.y() - _FN))) - i;
+
+	double chi = _chi0 + 2.0 * atan(((p.y() - _FN) - (p.x() - _FE) * tan(j/2.0))
+	  / _twoRk0);
+	double lambda = j + 2.0 * i + _lambda0;
+
+	double psi0 = 0.5 * log((1.0 + sin(chi)) / (_c * (1.0 - sin(chi)))) / _n;
+	double phi1 = 2.0 * atan(pow(M_E, psi0)) - M_PI_2;
+
+	double psi, phi = phi1, prev = phi;
+	for (int i = 0; i < 8; i++) {
+		double sinPhi = sin(phi);
+		psi = log((tan(phi/2.0 + M_PI_4)) * pow((1.0 - _e * sinPhi)
+		  / (1.0 + _e * sinPhi), _e/2.0));
+		phi = phi - (psi - psi0) * cos(phi) * (1.0 - _es * sinPhi * sinPhi)
+		  / (1.0 - _es);
+		if (fabs(phi - prev) < 1e-10)
+			break;
+		prev = phi;
+	}
+
+	return Coordinates(rad2deg((lambda - _lambda0) / _n + _lambda0),
+	  rad2deg(phi));
+}

src/map/obliquestereographic.h

@@ -0,0 +1,29 @@
+#ifndef OBLIQUESTEREOGRAPHIC_H
+#define OBLIQUESTEREOGRAPHIC_H
+
+#include "ct.h"
+
+class Ellipsoid;
+
+class ObliqueStereographic : public CT
+{
+public:
+	ObliqueStereographic(const Ellipsoid *ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double scale, double falseEasting,
+	  double falseNorthing);
+	virtual CT *clone() const {return new ObliqueStereographic(*this);}
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _e, _es;
+	double _chi0, _sinChi0, _cosChi0;
+	double _lambda0;
+	double _n;
+	double _c;
+	double _FE, _FN;
+	double _twoRk0,	_g, _h;
+};
+
+#endif // OBLIQUESTEREOGRAPHIC_H

src/map/projection.cpp

@@ -7,6 +7,7 @@
 #include "lambertazimuthal.h"
 #include "krovak.h"
 #include "polarstereographic.h"
+#include "obliquestereographic.h"
 #include "latlon.h"
 #include "gcs.h"
 #include "pcs.h"
@@ -22,6 +23,7 @@ Projection::Method::Method(int id)
 		case 9802:
 		case 9804:
 		case 9807:
+		case 9809:
 		case 9815:
 		case 9819:
 		case 9820:
@@ -72,6 +74,11 @@ Projection::Projection(const PCS *pcs) : _gcs(pcs->gcs()), _units(pcs->units()),
 			  setup.longitudeOrigin(), setup.scale(), setup.falseEasting(),
 			  setup.falseNorthing());
 			break;
+		case 9809:
+			_ct = new ObliqueStereographic(ellipsoid, setup.latitudeOrigin(),
+			  setup.longitudeOrigin(), setup.scale(), setup.falseEasting(),
+			  setup.falseNorthing());
+			break;
 		case 9819:
 			_ct = new Krovak(ellipsoid, setup.standardParallel1(),
 			  setup.standardParallel2(), setup.scale(), setup.latitudeOrigin(),