Fix handling of false easting and northing from proj4 string

Previous code did not respect the given falsings. Also includes
minor refactoring where the falsing-code was moved to parent class.

himan-lib/include/grid.h

@@ -276,6 +276,10 @@ class regular_grid : public grid
 	virtual std::unique_ptr<OGRSpatialReference> SpatialReference() const;
 	virtual std::vector<point> GridPointsInProjectionSpace() const override;
 
+	static std::pair<double, double> FalseEastingAndNorthing(const OGRSpatialReference* spRef,
+	                                                         OGRCoordinateTransformation* llToProjXForm,
+	                                                         const point& firstPoint, bool firstPointIsProjected);
+
    protected:
 	regular_grid(HPGridType gridType, HPScanningMode scMode, double di, double dj, size_t ni, size_t nj,
 	             bool uvRelativeToGrid, const std::string& theName = "");

himan-lib/source/grid.cpp

@@ -4,6 +4,7 @@
  */
 
 #include "grid.h"
+#include "util.h"
 #include <cpl_conv.h>
 #include <ogr_geometry.h>
 #include <ogr_spatialref.h>
@@ -593,6 +594,50 @@ std::vector<point> regular_grid::XY(const regular_grid& target) const
 	return sourceXY;
 }
 
+std::pair<double, double> regular_grid::FalseEastingAndNorthing(const OGRSpatialReference* spRef,
+                                                                OGRCoordinateTransformation* llToProjXForm,
+                                                                const point& firstPoint, bool firstPointIsProjected)
+{
+	// First check if spatial reference already has false easting and northing,
+	// if not, calculate them from the first point of the grid
+	double fe = spRef->GetProjParm(SRS_PP_FALSE_EASTING, 0.0);
+	double fn = spRef->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0);
+
+	if (fe != 0.0 || fn != 0.0)
+	{
+		return std::make_pair(fe, fn);
+	}
+
+	// First get first point coordinates in projected space. Transform from
+	// latlon if necessary.
+	double lat = firstPoint.Y(), lon = firstPoint.X();
+
+	if (firstPointIsProjected == false)
+	{
+		if (!llToProjXForm->Transform(1, &lon, &lat))
+		{
+			logger logr("regular_grid");
+			logr.Error("Failed to get false easting and northing");
+			return std::make_pair(MissingDouble(), MissingDouble());
+		}
+	}
+
+	// HIMAN-336: limit false easting/northing accuracy to four decimal places (millimeters)
+
+	lon = util::round(lon, 4);
+	lat = util::round(lat, 4);
+
+	if (fabs(lon) < 1e-4 and fabs(lat) < 1e-4)
+	{
+		return std::make_pair(0.0, 0.0);
+	}
+
+	fe = spRef->GetProjParm(SRS_PP_FALSE_EASTING, 0.0) - lon;
+	fn = spRef->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0) - lat;
+
+	return std::make_pair(fe, fn);
+}
+
 //--------------- irregular grid
 
 irregular_grid::irregular_grid(HPGridType type) : grid(kIrregularGrid, type, false)

himan-lib/source/lambert_conformal_grid.cpp

@@ -1,6 +1,5 @@
 #include "lambert_conformal_grid.h"
 #include "info.h"
-#include "util.h"
 #include <functional>
 #include <ogr_spatialref.h>
 
@@ -74,23 +73,10 @@ void lambert_conformal_grid::CreateCoordinateTransformations(const point& firstP
 	itsXYToLatLonTransformer = std::unique_ptr<OGRCoordinateTransformation>(
 	    OGRCreateCoordinateTransformation(itsSpatialReference.get(), geogCS.get()));
 
-	double lat = firstPoint.Y(), lon = firstPoint.X();
+	const auto [fe, fn] = regular_grid::FalseEastingAndNorthing(
+	    itsSpatialReference.get(), itsLatLonToXYTransformer.get(), firstPoint, firstPointIsProjected);
 
-	if (firstPointIsProjected == false)
-	{
-		if (!itsLatLonToXYTransformer->Transform(1, &lon, &lat))
-		{
-			itsLogger.Error("Failed to get false easting and northing");
-			return;
-		}
-	}
-
-	// HIMAN-336: limit false easting/northing accuracy to four decimal places (millimeters)
-
-	lon = util::round(lon, 4);
-	lat = util::round(lat, 4);
-
-	if (fabs(lon) < 1e-4 and fabs(lat) < 1e-4)
+	if (fabs(fe) < 1e-4 && fabs(fn) < 1e-4)
 	{
 		return;
 	}
@@ -99,9 +85,6 @@ void lambert_conformal_grid::CreateCoordinateTransformations(const point& firstP
 	const double parallel1 = StandardParallel1();
 	const double parallel2 = StandardParallel2();
 
-	const double fe = itsSpatialReference->GetProjParm(SRS_PP_FALSE_EASTING, 0.0) - lon;
-	const double fn = itsSpatialReference->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0) - lat;
-
 	// need to recreate spatial reference to get SetLCC to accept new falsings
 	itsSpatialReference = std::unique_ptr<OGRSpatialReference>(itsSpatialReference->CloneGeogCS());
 	if (itsSpatialReference->SetLCC(parallel1, parallel2, parallel1, orientation, fe, fn) != OGRERR_NONE)

himan-lib/source/lambert_equal_area_grid.cpp

@@ -1,6 +1,5 @@
 #include "lambert_equal_area_grid.h"
 #include "info.h"
-#include "util.h"
 #include <functional>
 #include <ogr_spatialref.h>
 
@@ -66,29 +65,16 @@ void lambert_equal_area_grid::CreateCoordinateTransformations(const point& first
 	itsXYToLatLonTransformer = std::unique_ptr<OGRCoordinateTransformation>(
 	    OGRCreateCoordinateTransformation(itsSpatialReference.get(), geogCS.get()));
 
-	double lat = firstPoint.Y(), lon = firstPoint.X();
+	const auto [fe, fn] = regular_grid::FalseEastingAndNorthing(
+	    itsSpatialReference.get(), itsLatLonToXYTransformer.get(), firstPoint, firstPointIsProjected);
 
-	if (firstPointIsProjected == false)
-	{
-		if (!itsLatLonToXYTransformer->Transform(1, &lon, &lat))
-		{
-			itsLogger.Fatal("Failed to get false easting and northing");
-			himan::Abort();
-		}
-	}
-
-	lon = util::round(lon, 4);
-	lat = util::round(lat, 4);
-
-	if (fabs(lon) < 1e-4 and fabs(lat) < 1e-4)
+	if (fabs(fe) < 1e-4 && fabs(fn) < 1e-4)
 	{
 		return;
 	}
 
 	const double orientation = Orientation();
 	const double parallel = StandardParallel();
-	const double fe = itsSpatialReference->GetProjParm(SRS_PP_FALSE_EASTING, 0.0) - lon;
-	const double fn = itsSpatialReference->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0) - lat;
 
 	itsSpatialReference = std::unique_ptr<OGRSpatialReference>(itsSpatialReference->CloneGeogCS());
 	if (itsSpatialReference->SetLAEA(parallel, orientation, fe, fn) != OGRERR_NONE)

himan-lib/source/transverse_mercator_grid.cpp

@@ -1,6 +1,5 @@
 #include "transverse_mercator_grid.h"
 #include "info.h"
-#include "util.h"
 #include <functional>
 #include <ogr_spatialref.h>
 
@@ -69,30 +68,17 @@ void transverse_mercator_grid::CreateCoordinateTransformations(const point& firs
 	itsXYToLatLonTransformer = std::unique_ptr<OGRCoordinateTransformation>(
 	    OGRCreateCoordinateTransformation(itsSpatialReference.get(), geogCS.get()));
 
-	double lat = firstPoint.Y(), lon = firstPoint.X();
+	const auto [fe, fn] = regular_grid::FalseEastingAndNorthing(
+	    itsSpatialReference.get(), itsLatLonToXYTransformer.get(), firstPoint, firstPointIsProjected);
 
-	if (firstPointIsProjected == false)
-	{
-		if (!itsLatLonToXYTransformer->Transform(1, &lon, &lat))
-		{
-			itsLogger.Fatal("Failed to get false easting and northing");
-			himan::Abort();
-		}
-	}
-
-	lon = util::round(lon, 4);
-	lat = util::round(lat, 4);
-
-	if (fabs(lon) < 1e-4 and fabs(lat) < 1e-4)
+	if (fabs(fe) < 1e-4 && fabs(fn) < 1e-4)
 	{
 		return;
 	}
 
 	const double orientation = Orientation();
 	const double parallel = StandardParallel();
 	const double scale = Scale();
-	const double fe = itsSpatialReference->GetProjParm(SRS_PP_FALSE_EASTING, 0.0) - lon;
-	const double fn = itsSpatialReference->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0) - lat;
 
 	itsSpatialReference = std::unique_ptr<OGRSpatialReference>(itsSpatialReference->CloneGeogCS());
 	if (itsSpatialReference->SetTM(parallel, orientation, scale, fe, fn) != OGRERR_NONE)