Merge pull request #3 from Ritchie333/france

Merge support for French projections

CMakeLists.txt

@@ -114,6 +114,11 @@ add_executable(clearbounds
     src/gbos1936/Gbos1936.cpp
 )
 
+add_executable(testcas
+    src/TestCas.cpp
+    src/gbos1936/Gbos1936.cpp
+)
+
 target_include_directories(warp-sabre PUBLIC
     ${PROJECT_SOURCE_DIR}/libmorph
     ${PROJECT_SOURCE_DIR}/newmat-10

src/CopyPixels.cpp

@@ -2,6 +2,7 @@
 #include "ImgMagick.h"
 #include <vector>
 #include <iostream>
+#include <string>
 #include <stdio.h>
 #include "CopyPixels.h"
 #include "Tile.h"
@@ -11,36 +12,33 @@ using namespace std;
 
 class HelmertConverter gConverter;
 
+#define BEGIN_MASK_MAP(type)	string sType = type;
+#define	MASK_ENTRY(name,class)	if( sType == name ) { mask = new class; }
+#define END_MASK_MAP()
+
 CopyPixels *CopyPixels::Create(const char *type)
 {
 	CopyPixels *mask = NULL;
-	switch (type[0])
-	{
-	case 'G':
-		mask = new CopyPixelsWithOsMask();
-		break;
-	case 'O':
-		mask = new CopyPixelsWithRawMask();
-		break;
-	case 'C':
-		mask = new CopyPixelsWithCassini();
-		break;
-	case 'B':
-		mask = new CopyPixelsWithBonne();
-		break;
-	case 'I':
-		mask = new CopyPixelsWithOsI();
-		break;
-	case 'R':
-		mask = new CopyPixelsWithIrishBonne();
-		break;
-	case 'M':
-		mask = new CopyPixelsWithMercator();
-		break;
-	default:
+
+	BEGIN_MASK_MAP(type)
+	MASK_ENTRY("G", CopyPixelsWithOsMask )
+	MASK_ENTRY("O", CopyPixelsWithRawMask )
+	MASK_ENTRY("OS", CopyPixelsWithRawMask )
+	MASK_ENTRY("C", CopyPixelsWithCassini )
+	MASK_ENTRY("CAS", CopyPixelsWithCassini )
+	MASK_ENTRY("B", CopyPixelsWithBonne )
+	MASK_ENTRY("SB", CopyPixelsWithBonne )
+	MASK_ENTRY("I", CopyPixelsWithOsI )
+	MASK_ENTRY("R", CopyPixelsWithIrishBonne )
+	MASK_ENTRY("IB", CopyPixelsWithIrishBonne )
+	MASK_ENTRY("FB", CopyPixelsWithFrenchBonne )
+	MASK_ENTRY("PM", CopyPixelsWithParisMercator )
+	MASK_ENTRY("M", CopyPixelsWithMercator )
+	END_MASK_MAP()
+
+	if( !mask ) {
 		cout << "Invalid projection type " << type << endl;
 		throw 3001;
-		break;
 	}
 	return mask;
 }
@@ -180,6 +178,42 @@ void CopyPixelsWithMercator::UpdateBoundingBox(const char *mapref)
 	}
 }
 
+void CopyPixelsWithParisMercator::UpdateBoundingBox(const char *mapref)
+{
+	double dgLng = 0, dgLat = 0, dLng = 0, dLat = 0;
+	if (2 == sscanf(mapref, "%lf:%lf", &dgLat, &dgLng))
+	{
+		gConverter.ConvertParisToWgs84( dgLat * 0.9, dgLng * 0.9, dLat, dLng );
+
+		if (!boxset || gsouth > dLat)
+			gsouth = dLat;
+		if (!boxset || gnorth < dLat)
+			gnorth = dLat;
+		if (!boxset || geast < dLng)
+			geast = dLng;
+		if (!boxset || gwest > dLng)
+			gwest = dLng;
+		gVertx.clear();
+		gVerty.clear();
+
+		gVertx.push_back(gwest);
+		gVertx.push_back(geast);
+		gVertx.push_back(geast);
+		gVertx.push_back(gwest);
+
+		gVerty.push_back(gnorth);
+		gVerty.push_back(gnorth);
+		gVerty.push_back(gsouth);
+		gVerty.push_back(gsouth);
+
+		boxset = 1;
+	}
+	else
+	{
+		ThrowError("Invalid mercator reference", mapref);
+	}
+}
+
 int CopyPixelsWithOsMask::CheckIfInBox(double lat, double lon)
 {
 	double pnorth, peast, palt;
@@ -257,6 +291,21 @@ int CopyPixelsWithIrishBonne::CheckIfInBox(double lat, double lon)
 	return 1;
 }
 
+int CopyPixelsWithFrenchBonne::CheckIfInBox(double lat, double lon)
+{
+	double pnorth, peast;
+	gConverter.ConvertWgs84ToBnF(lat, lon, 0.0, peast, pnorth);
+	if (pnorth < gsouth)
+		return 0;
+	if (pnorth > gnorth)
+		return 0;
+	if (peast < gwest)
+		return 0;
+	if (peast > geast)
+		return 0;
+	return 1;
+}
+
 int CopyPixelsWithMercator::CheckIfInBox(double lat, double lon)
 {
 	int c = 0, i = 0, j = 0;

src/CopyPixels.h

@@ -68,14 +68,28 @@ class CopyPixelsWithIrishBonne : public CopyPixelsWithBonne
 	virtual int CheckIfInBox(double lat, double lon);
 };
 
+class CopyPixelsWithFrenchBonne : public CopyPixelsWithBonne
+{
+public:
+	CopyPixelsWithFrenchBonne() {}
+	virtual int CheckIfInBox(double lat, double lon);
+};
+
 class CopyPixelsWithMercator : public CopyPixels
 {
 public:
 	CopyPixelsWithMercator() {}
 	virtual void UpdateBoundingBox(const char *mapref);
 	virtual int CheckIfInBox(double lat, double lon);
 
-private:
+protected:
 	std::vector<double> gVertx;
 	std::vector<double> gVerty;
 };
+
+class CopyPixelsWithParisMercator : public CopyPixelsWithMercator
+{
+public:
+	CopyPixelsWithParisMercator() {}
+	virtual void UpdateBoundingBox(const char *mapref);
+};

src/TestCas.cpp

@@ -36,5 +36,22 @@ int main()
   gConverter.ConvertWgs84ToOsi(53.5,-8,0.0, e, n, h);
   std::cout << e << "," << n << std::endl;
 
+  std::cout << "French Bonne" << std::endl;
+  //gConverter.ConvertWgs84ToBnF(48.86, 2.33722222, 0.0, e, n );
+  gConverter.ConvertWgs84ToBnF(48.836439, 2.336506, 0.0, e, n );
+  std::cout << "Point zero " << e << "," << n << std::endl;
+  gConverter.ConvertWgs84ToBnF(47.957873, -0.543229, 0.0, e, n );
+  std::cout << "47.957873, -0.543229 " << e << "," << n << std::endl;
+  gConverter.ConvertBnFToWgs84( -215, -93.5, 0.0, lat, lon, alt );
+  std::cout << "-215, -93.75 " << lat << "," << lon << std::endl;
+
+  gConverter.ConvertWgs84ToBnF(47.955808, 5.214928, 0.0, e, n );
+  std::cout << "47.955808, 5.214928 " << e << "," << n << std::endl;
+
+  gConverter.ConvertBnFToWgs84( 215, -93.5, 0.0, lat, lon, alt );
+  std::cout << "215, -93.75 " << lat << "," << lon << std::endl;
+
+
+
   return 0;
 }

src/WarpOs.cpp

@@ -45,6 +45,7 @@ class PolyProjectArgs
 		Cassini,
 		BonneS,
 		BonneI,
+		BonneF,
 		OSI
 	} ProjType;
 
@@ -108,6 +109,8 @@ const Point ProjRefToOutImg(const Point& ref, PolyProjectArgs::ProjType projType
 	case PolyProjectArgs::BonneI:
 		gFallbackConverter.ConvertBnIToWgs84(ref.x, ref.y, 0.0, lat, lon, alt);
 		break;
+	case PolyProjectArgs::BonneF:
+		gFallbackConverter.ConvertBnFToWgs84(ref.x, ref.y, 0.0, lat, lon, alt );
 	}
 	if (projType != PolyProjectArgs::OSGB && !args->mercatorOut)
 	{
@@ -295,6 +298,10 @@ int main(int argc, char *argv[])
 		{
 			projType = PolyProjectArgs::BonneI;
 		}
+		if( inproj == "bonnef" )
+		{
+			projType = PolyProjectArgs::BonneF;
+		}
 		if (inproj == "osi")
 		{
 			projType = PolyProjectArgs::OSI;
@@ -363,6 +370,35 @@ int main(int argc, char *argv[])
 				if (iy > east || !setBox)
 					east = iy;
 
+				if (projType != PolyProjectArgs::Mercator)
+				{
+					cout << "Cannot take lat lon input when using OS as input projection" << endl;
+					exit(0);
+				}
+				else
+					srcImgToRef.AddPoint(imgX, imgY, ix, iy);
+
+				setBox = 1;
+			}
+			if (strcmp(line[0].GetVals(), "g") == 0)	// French maps are in grads
+			{
+				double imgX = line[3].GetVald();
+				double imgY = line[4].GetVald();
+				double ix = line[1].GetVald() / 0.9;
+				double iy = line[2].GetVald() / 0.9;
+				double lat = 0.0, lon = 0.0;
+				gConverter.ConvertParisToWgs84( ix, iy, lat, lon );
+
+				if (lat < south || !setBox)
+					south = lat;
+				if (lat > north || !setBox)
+					north = lat;
+				if (lon < west || !setBox)
+					west = lon;
+				if (lon > east || !setBox)
+					east = lon;
+				setBox = 1;
+
 				if (projType != PolyProjectArgs::Mercator)
 				{
 					cout << "Cannot take lat lon input when using OS as input projection" << endl;
@@ -465,6 +501,12 @@ int main(int argc, char *argv[])
 				dNorthing = line[2].GetVald();
 				gFallbackConverter.ConvertBnIToWgs84(dEasting, dNorthing, 0.0, lat, lon, alt);
 			}
+			if (strcmp(line[0].GetVals(), "bnf") == 0)
+			{
+				dEasting = line[1].GetVald();
+				dNorthing = line[2].GetVald();
+				gFallbackConverter.ConvertBnFToWgs84(dEasting, dNorthing, 0.0, lat, lon, alt);
+			}
 			if (lat != -1.0 && lon != -1.0)
 			{
 
@@ -503,6 +545,61 @@ int main(int argc, char *argv[])
 				// cout << "gbos " << dEasting << "," << dNorthing << "\t" << lat << "," << lon << endl;
 			}
 		}
+
+		if (line.NumVals() == 6) // Read in UTM data
+		{
+			int zone = 0;
+			double dEasting = 0, dNorthing = 0;
+			double imgX = line[4].GetVald();
+			double imgY = line[5].GetVald();
+			double lat = -1.0, lon = -1.0;
+			if (strcmp(line[0].GetVals(), "utm") == 0)
+			{
+				zone = line[2].GetVali();
+				dEasting = line[3].GetVald();
+				dNorthing = line[4].GetVald();
+
+				gConverter.ConvertUTM50ToWgs84( zone, dEasting, dNorthing, lat, lon);
+			}
+			if (lat != -1.0 && lon != -1.0)
+			{
+
+				cout << "conv " << zone << ":" << dEasting << "," << dNorthing << "->" << lat << "," << lon << endl;
+
+				// Add point to transform constraints
+				if (projType != PolyProjectArgs::Mercator)
+					srcImgToRef.AddPoint(imgX, imgY, dEasting, dNorthing);
+				else
+					srcImgToRef.AddPoint(imgX, imgY, lat, lon);
+
+				if (mercatorOut)
+				{
+					if (lat < south || !setBox)
+						south = lat;
+					if (lat > north || !setBox)
+						north = lat;
+					if (lon < west || !setBox)
+						west = lon;
+					if (lon > east || !setBox)
+						east = lon;
+					setBox = 1;
+				}
+				if (gbosOut)
+				{
+					if (dNorthing < south || !setBox)
+						south = dNorthing;
+					if (dNorthing > north || !setBox)
+						north = dNorthing;
+					if (dEasting < west || !setBox)
+						west = dEasting;
+					if (dEasting > east || !setBox)
+						east = dEasting;
+					setBox = 1;
+				}
+				// cout << "gbos " << dEasting << "," << dNorthing << "\t" << lat << "," << lon << endl;
+			}
+		}
+
 	}
 
 	if (corners.size() == 0)