From 3cf09379f47d20246f51a854c5b6f0828b694316 Mon Sep 17 00:00:00 2001 From: Mikko Partio Date: Wed, 20 Nov 2024 09:26:55 +0200 Subject: [PATCH] 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 | 4 ++ himan-lib/source/grid.cpp | 45 +++++++++++++++++++ himan-lib/source/lambert_conformal_grid.cpp | 23 ++-------- himan-lib/source/lambert_equal_area_grid.cpp | 20 ++------- himan-lib/source/transverse_mercator_grid.cpp | 20 ++------- 5 files changed, 58 insertions(+), 54 deletions(-) diff --git a/himan-lib/include/grid.h b/himan-lib/include/grid.h index f3f9b3364..609c5b914 100644 --- a/himan-lib/include/grid.h +++ b/himan-lib/include/grid.h @@ -276,6 +276,10 @@ class regular_grid : public grid virtual std::unique_ptr SpatialReference() const; virtual std::vector GridPointsInProjectionSpace() const override; + static std::pair 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 = ""); diff --git a/himan-lib/source/grid.cpp b/himan-lib/source/grid.cpp index 84719179f..16db85a5f 100644 --- a/himan-lib/source/grid.cpp +++ b/himan-lib/source/grid.cpp @@ -4,6 +4,7 @@ */ #include "grid.h" +#include "util.h" #include #include #include @@ -593,6 +594,50 @@ std::vector regular_grid::XY(const regular_grid& target) const return sourceXY; } +std::pair 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) diff --git a/himan-lib/source/lambert_conformal_grid.cpp b/himan-lib/source/lambert_conformal_grid.cpp index 4ad8a0f04..fe0b5a425 100644 --- a/himan-lib/source/lambert_conformal_grid.cpp +++ b/himan-lib/source/lambert_conformal_grid.cpp @@ -1,6 +1,5 @@ #include "lambert_conformal_grid.h" #include "info.h" -#include "util.h" #include #include @@ -74,23 +73,10 @@ void lambert_conformal_grid::CreateCoordinateTransformations(const point& firstP itsXYToLatLonTransformer = std::unique_ptr( 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(itsSpatialReference->CloneGeogCS()); if (itsSpatialReference->SetLCC(parallel1, parallel2, parallel1, orientation, fe, fn) != OGRERR_NONE) diff --git a/himan-lib/source/lambert_equal_area_grid.cpp b/himan-lib/source/lambert_equal_area_grid.cpp index e355e9781..53b00552c 100644 --- a/himan-lib/source/lambert_equal_area_grid.cpp +++ b/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 #include @@ -66,29 +65,16 @@ void lambert_equal_area_grid::CreateCoordinateTransformations(const point& first itsXYToLatLonTransformer = std::unique_ptr( 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(itsSpatialReference->CloneGeogCS()); if (itsSpatialReference->SetLAEA(parallel, orientation, fe, fn) != OGRERR_NONE) diff --git a/himan-lib/source/transverse_mercator_grid.cpp b/himan-lib/source/transverse_mercator_grid.cpp index 36931294d..0844e2607 100644 --- a/himan-lib/source/transverse_mercator_grid.cpp +++ b/himan-lib/source/transverse_mercator_grid.cpp @@ -1,6 +1,5 @@ #include "transverse_mercator_grid.h" #include "info.h" -#include "util.h" #include #include @@ -69,21 +68,10 @@ void transverse_mercator_grid::CreateCoordinateTransformations(const point& firs itsXYToLatLonTransformer = std::unique_ptr( 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; } @@ -91,8 +79,6 @@ void transverse_mercator_grid::CreateCoordinateTransformations(const point& firs 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(itsSpatialReference->CloneGeogCS()); if (itsSpatialReference->SetTM(parallel, orientation, scale, fe, fn) != OGRERR_NONE)