Added support for EPSG 21781 PCS (Swiss grid)

Fixed swapped LCC1 and LCC2 GeoTIFF coordinate transformation codes

pkg/csv/gcs.csv

@@ -21,7 +21,7 @@ Canton Astro 1966,4716,6716,9122,7022,8901,9603,298,-304,-375
 Cape,4222,6222,9122,7012,8901,9603,-136,-108,-292
 Cape Canaveral,4717,6717,9122,7008,8901,9603,-2,150,181
 Carthage,4223,6223,9122,7012,8901,9603,-263,6,431
-CH-1903,4149,6149,9122,7004,8901,9603,674,15,405
+CH-1903,4149,6149,9122,7004,8901,9603,674.374,15.056,405.343
 Chatham 1971,4672,6672,9122,7022,8901,9603,175,-38,113
 Chua Astro,4224,6224,9122,7022,8901,9603,-134,229,-29
 Corrego Alegre,4225,6225,9122,7022,8901,9603,-206,172,-6

src/map/geotiff.cpp

@@ -33,14 +33,32 @@
 #define ProjNatOriginLatGeoKey       3081
 #define ProjFalseEastingGeoKey       3082
 #define ProjFalseNorthingGeoKey      3083
+#define ProjCenterLongGeoKey         3088
+#define ProjCenterLatGeoKey          3089
 #define ProjScaleAtNatOriginGeoKey   3092
+#define ProjScaleAtCenterGeoKey      3093
+#define ProjAzimuthAngleGeoKey       3094
+#define ProjRectifiedGridAngleGeoKey 3096
 
 #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 IS_SET(map, key) \
 	((map).contains(key) && (map).value(key).SHORT != 32767)
+#define TO_M(lu, val) \
+	(std::isnan(val) ? val : (lu).toMeters(val))
+#define TO_DEG(au, val) \
+	(std::isnan(val) ? val : (au).toDegrees(val))
 
 #define ARRAY_SIZE(a) \
 	(sizeof(a) / sizeof(*a))
@@ -235,6 +253,11 @@ bool GeoTIFF::readKeys(TIFFFile &file, Ctx &ctx, QMap<quint16, Value> &kv) const
 			case ProjFalseNorthingGeoKey:
 			case ProjStdParallel1GeoKey:
 			case ProjStdParallel2GeoKey:
+			case ProjCenterLongGeoKey:
+			case ProjCenterLatGeoKey:
+			case ProjScaleAtCenterGeoKey:
+			case ProjAzimuthAngleGeoKey:
+			case ProjRectifiedGridAngleGeoKey:
 				if (!readGeoValue(file, ctx.values, entry.ValueOffset,
 				  value.DOUBLE))
 					return false;
@@ -297,17 +320,19 @@ Projection::Method GeoTIFF::method(QMap<quint16, Value> &kv)
 	}
 
 	switch (kv.value(ProjCoordTransGeoKey).SHORT) {
-		case 1:
+		case CT_TransverseMercator:
 			return Projection::Method(9807);
-		case 7:
+		case CT_ObliqueMercator:
+			return Projection::Method(9815);
+		case CT_Mercator:
 			return Projection::Method(1024);
-		case 8:
+		case CT_LambertConfConic_2SP:
-			return Projection::Method(9801);
-		case 9:
 			return Projection::Method(9802);
-		case 10:
+		case CT_LambertConfConic_1SP:
+			return Projection::Method(9801);
+		case CT_LambertAzimEqualArea:
 			return Projection::Method(9820);
-		case 11:
+		case CT_AlbersEqualArea:
 			return Projection::Method(9822);
 		default:
 			_errorString = QString("%1: unknown coordinate transformation method")
@@ -341,6 +366,8 @@ bool GeoTIFF::projectedModel(QMap<quint16, Value> &kv)
 		_projection = Projection(pcs->gcs(), pcs->method(), pcs->setup(),
 		  pcs->units());
 	} else {
+		double lat0, lon0, scale, fe, fn, sp1, sp2;
+
 		const GCS *g = gcs(kv);
 		if (!g)
 			return false;
@@ -358,15 +385,47 @@ bool GeoTIFF::projectedModel(QMap<quint16, Value> &kv)
 			return false;
 		}
 
-		Projection::Setup setup(
+		if (kv.contains(ProjNatOriginLatGeoKey))
-		  au.toDegrees(kv.value(ProjNatOriginLatGeoKey).DOUBLE),
+			lat0 = kv.value(ProjNatOriginLatGeoKey).DOUBLE;
-		  au.toDegrees(kv.value(ProjNatOriginLongGeoKey).DOUBLE),
+		else if (kv.contains(ProjCenterLatGeoKey))
-		  kv.value(ProjScaleAtNatOriginGeoKey).DOUBLE,
+			lat0 = kv.value(ProjCenterLatGeoKey).DOUBLE;
-		  lu.toMeters(kv.value(ProjFalseEastingGeoKey).DOUBLE),
+		else
-		  lu.toMeters(kv.value(ProjFalseNorthingGeoKey).DOUBLE),
+			lat0 = NAN;
-		  au.toDegrees(kv.value(ProjStdParallel1GeoKey).DOUBLE),
+		if (kv.contains(ProjNatOriginLongGeoKey))
-		  au.toDegrees(kv.value(ProjStdParallel2GeoKey).DOUBLE)
+			lon0 = kv.value(ProjNatOriginLongGeoKey).DOUBLE;
-		);
+		else if (kv.contains(ProjCenterLongGeoKey))
+			lon0 = kv.value(ProjCenterLongGeoKey).DOUBLE;
+		else
+			lon0 = NAN;
+		if (kv.contains(ProjScaleAtNatOriginGeoKey))
+			scale = kv.value(ProjScaleAtNatOriginGeoKey).DOUBLE;
+		else if (kv.contains(ProjScaleAtCenterGeoKey))
+			scale = kv.value(ProjScaleAtCenterGeoKey).DOUBLE;
+		else
+			scale = NAN;
+		if (kv.contains(ProjStdParallel1GeoKey))
+			sp1 = kv.value(ProjStdParallel1GeoKey).DOUBLE;
+		else if (kv.contains(ProjAzimuthAngleGeoKey))
+			sp1 = kv.value(ProjAzimuthAngleGeoKey).DOUBLE;
+		else
+			sp1 = NAN;
+		if (kv.contains(ProjStdParallel2GeoKey))
+			sp2 = kv.value(ProjStdParallel2GeoKey).DOUBLE;
+		else if (kv.contains(ProjRectifiedGridAngleGeoKey))
+			sp2 = kv.value(ProjRectifiedGridAngleGeoKey).DOUBLE;
+		else
+			sp2 = NAN;
+		if (kv.contains(ProjFalseNorthingGeoKey))
+			fn = kv.value(ProjFalseNorthingGeoKey).DOUBLE;
+		else
+			fn = NAN;
+		if (kv.contains(ProjFalseEastingGeoKey))
+			fe = kv.value(ProjFalseEastingGeoKey).DOUBLE;
+		else
+			fe = NAN;
+
+		Projection::Setup setup(TO_DEG(au, lat0), TO_DEG(au, lon0), scale,
+		  TO_M(lu, fe), TO_M(lu, fn), TO_DEG(au, sp1), TO_DEG(au, sp2));
 
 		_projection = Projection(g, m, setup, lu);
 	}

src/map/pcs.cpp

@@ -22,6 +22,7 @@ static bool parameter(int key, double val, int units, Projection::Setup &setup)
 {
 	switch (key) {
 		case 8801:
+		case 8811:
 		case 8821:
 			{AngularUnits au(units);
 			if (au.isNull())
@@ -29,6 +30,7 @@ static bool parameter(int key, double val, int units, Projection::Setup &setup)
 			setup.setLatitudeOrigin(au.toDegrees(val));}
 			return true;
 		case 8802:
+		case 8812:
 		case 8822:
 			{AngularUnits au(units);
 			if (au.isNull())
@@ -36,9 +38,11 @@ static bool parameter(int key, double val, int units, Projection::Setup &setup)
 			setup.setLongitudeOrigin(au.toDegrees(val));}
 			return true;
 		case 8805:
+		case 8815:
 			setup.setScale(val);
 			return true;
 		case 8806:
+		case 8816:
 		case 8826:
 			{LinearUnits lu(units);
 			if (lu.isNull())
@@ -46,18 +50,21 @@ static bool parameter(int key, double val, int units, Projection::Setup &setup)
 			setup.setFalseEasting(lu.toMeters(val));}
 			return true;
 		case 8807:
+		case 8817:
 		case 8827:
 			{LinearUnits lu(units);
 			if (lu.isNull())
 				return false;
 			setup.setFalseNorthing(lu.toMeters(val));}
 			return true;
+		case 8813:
 		case 8823:
 			{AngularUnits au(units);
 			if (au.isNull())
 				return false;
 			setup.setStandardParallel1(au.toDegrees(val));}
 			return true;
+		case 8814:
 		case 8824:
 			{AngularUnits au(units);
 			if (au.isNull())
@@ -74,7 +81,7 @@ static int projectionSetup(const QList<QByteArray> &list,
 {
 	bool r1, r2, r3;
 
-	for (int i = 6; i < 24; i += 3) {
+	for (int i = 6; i < 27; i += 3) {
 		QString ks = list[i].trimmed();
 		if (ks.isEmpty())
 			break;
@@ -129,7 +136,7 @@ void PCS::loadList(const QString &path)
 
 		QByteArray line = file.readLine();
 		QList<QByteArray> list = line.split(',');
-		if (list.size() != 24) {
+		if (list.size() != 27) {
 			qWarning("%s: %d: Format error", qPrintable(path), ln);
 			continue;
 		}

src/map/projection.cpp

@@ -16,6 +16,7 @@ Projection::Method::Method(int id)
 		case 9801:
 		case 9802:
 		case 9807:
+		case 9815:
 		case 9820:
 		case 9822:
 		case 9841:
@@ -32,23 +33,24 @@ Projection::Projection(const GCS *gcs, const Method &method, const Setup &setup,
 	const Ellipsoid *ellipsoid = _gcs->datum().ellipsoid();
 
 	switch (method.id()) {
-		case 9807:
-			_ct = new TransverseMercator(ellipsoid, setup.latitudeOrigin(),
-			  setup.longitudeOrigin(), setup.scale(), setup.falseEasting(),
-			  setup.falseNorthing());
-			break;
 		case 1024:
 		case 9841:
 			_ct = new Mercator();
 			break;
+		case 9801:
+		case 9815: // Oblique mercator aproximation using LCC1
+			_ct = new LambertConic1(ellipsoid, setup.latitudeOrigin(),
+			  setup.longitudeOrigin(), setup.scale(), setup.falseEasting(),
+			  setup.falseNorthing());
+			break;
 		case 9802:
 			_ct = new LambertConic2(ellipsoid, setup.standardParallel1(),
 			  setup.standardParallel2(), setup.latitudeOrigin(),
 			  setup.longitudeOrigin(), setup.falseEasting(),
 			  setup.falseNorthing());
 			break;
-		case 9801:
+		case 9807:
-			_ct = new LambertConic1(ellipsoid, setup.latitudeOrigin(),
+			_ct = new TransverseMercator(ellipsoid, setup.latitudeOrigin(),
 			  setup.longitudeOrigin(), setup.scale(), setup.falseEasting(),
 			  setup.falseNorthing());
 			break;
@@ -66,6 +68,8 @@ Projection::Projection(const GCS *gcs, const Method &method, const Setup &setup,
 		default:
 			_ct = 0;
 	}
+
+	Q_ASSERT(_ct != 0);
 }
 
 Projection::Projection(const GCS *gcs) : _gcs(gcs), _geographic(true)