From 926bbd2064d17bfbc46854dde8750e18bbf6f766 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 5 Jun 2024 19:21:15 +0200 Subject: [PATCH] BoundCRS::exportToPROJ(): handle case of NADCON conus grid Fixes #4154 --- include/proj/coordinateoperation.hpp | 2 +- include/proj/crs.hpp | 3 +- src/iso19111/crs.cpp | 17 ++++++-- src/iso19111/operation/singleoperation.cpp | 35 +++++++++++++++- test/unit/test_crs.cpp | 48 ++++++++++++++++++++++ 5 files changed, 97 insertions(+), 8 deletions(-) diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp index 3d288ddd82..d4b666e207 100644 --- a/include/proj/coordinateoperation.hpp +++ b/include/proj/coordinateoperation.hpp @@ -1702,7 +1702,7 @@ class PROJ_GCC_DLL Transformation : public SingleOperation { PROJ_PRIVATE : //! @cond Doxygen_Suppress PROJ_INTERNAL const std::string & - getNTv2Filename() const; + getPROJ4NadgridsCompatibleFilename() const; PROJ_FOR_TEST std::vector getTOWGS84Parameters() const; // throw(io::FormattingException) diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp index d9c3846573..a6efcd4fe1 100644 --- a/include/proj/crs.hpp +++ b/include/proj/crs.hpp @@ -1064,7 +1064,8 @@ class PROJ_GCC_DLL BoundCRS final : public CRS, PROJ_INTERNAL BoundCRSNNPtr shallowCloneAsBoundCRS() const; PROJ_INTERNAL bool isTOWGS84Compatible() const; - PROJ_INTERNAL std::string getHDatumPROJ4GRIDS() const; + PROJ_INTERNAL std::string + getHDatumPROJ4GRIDS(const io::DatabaseContextPtr &databaseContext) const; PROJ_INTERNAL std::string getVDatumPROJ4GRIDS(const crs::GeographicCRS *geogCRSOfCompoundCRS, const char **outGeoidCRSValue) const; diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 6f3c636b55..5fd9ba1b2b 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -5886,9 +5886,16 @@ bool BoundCRS::isTOWGS84Compatible() const { // --------------------------------------------------------------------------- -std::string BoundCRS::getHDatumPROJ4GRIDS() const { +std::string BoundCRS::getHDatumPROJ4GRIDS( + const io::DatabaseContextPtr &databaseContext) const { if (ci_equal(d->hubCRS()->nameStr(), "WGS 84")) { - return d->transformation()->getNTv2Filename(); + if (databaseContext) { + return d->transformation() + ->substitutePROJAlternativeGridNames( + NN_NO_CHECK(databaseContext)) + ->getPROJ4NadgridsCompatibleFilename(); + } + return d->transformation()->getPROJ4NadgridsCompatibleFilename(); } return std::string(); } @@ -5945,7 +5952,8 @@ void BoundCRS::_exportToWKT(io::WKTFormatter *formatter) const { return; } - auto hdatumProj4GridName = getHDatumPROJ4GRIDS(); + auto hdatumProj4GridName = + getHDatumPROJ4GRIDS(formatter->databaseContext()); if (!hdatumProj4GridName.empty()) { formatter->setHDatumExtension(hdatumProj4GridName); d->baseCRS()->_exportToWKT(formatter); @@ -6039,7 +6047,8 @@ void BoundCRS::_exportToPROJString( crs_exportable->_exportToPROJString(formatter); formatter->setVDatumExtension(std::string(), std::string()); } else { - auto hdatumProj4GridName = getHDatumPROJ4GRIDS(); + auto hdatumProj4GridName = + getHDatumPROJ4GRIDS(formatter->databaseContext()); if (!hdatumProj4GridName.empty()) { formatter->setHDatumExtension(hdatumProj4GridName); crs_exportable->_exportToPROJString(formatter); diff --git a/src/iso19111/operation/singleoperation.cpp b/src/iso19111/operation/singleoperation.cpp index 19dd67da15..79b5815a9d 100644 --- a/src/iso19111/operation/singleoperation.cpp +++ b/src/iso19111/operation/singleoperation.cpp @@ -1867,9 +1867,40 @@ static const std::string &_getNTv2Filename(const SingleOperation *op, // --------------------------------------------------------------------------- //! @cond Doxygen_Suppress -const std::string &Transformation::getNTv2Filename() const { +const std::string &Transformation::getPROJ4NadgridsCompatibleFilename() const { - return _getNTv2Filename(this, false); + const std::string &filename = _getNTv2Filename(this, false); + if (!filename.empty()) { + return filename; + } + + if (method()->getEPSGCode() == EPSG_CODE_METHOD_NADCON) { + const auto &latitudeFileParameter = + parameterValue(EPSG_NAME_PARAMETER_LATITUDE_DIFFERENCE_FILE, + EPSG_CODE_PARAMETER_LATITUDE_DIFFERENCE_FILE); + const auto &longitudeFileParameter = + parameterValue(EPSG_NAME_PARAMETER_LONGITUDE_DIFFERENCE_FILE, + EPSG_CODE_PARAMETER_LONGITUDE_DIFFERENCE_FILE); + if (latitudeFileParameter && + latitudeFileParameter->type() == ParameterValue::Type::FILENAME && + longitudeFileParameter && + longitudeFileParameter->type() == ParameterValue::Type::FILENAME) { + return latitudeFileParameter->valueFile(); + } + } + + if (ci_equal(method()->nameStr(), + PROJ_WKT2_NAME_METHOD_HORIZONTAL_SHIFT_GTIFF)) { + const auto &fileParameter = parameterValue( + EPSG_NAME_PARAMETER_LATITUDE_LONGITUDE_DIFFERENCE_FILE, + EPSG_CODE_PARAMETER_LATITUDE_LONGITUDE_DIFFERENCE_FILE); + if (fileParameter && + fileParameter->type() == ParameterValue::Type::FILENAME) { + return fileParameter->valueFile(); + } + } + + return nullString; } //! @endcond diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp index 8ccafb199e..31a3551079 100644 --- a/test/unit/test_crs.cpp +++ b/test/unit/test_crs.cpp @@ -5209,6 +5209,54 @@ TEST(crs, boundCRS_projectedCRS_to_PROJ_string) { // --------------------------------------------------------------------------- +TEST(crs, boundCRS_nadcon_to_PROJ_string_usePROJAlternativeGridNames_false) { + + auto dbContext = DatabaseContext::create(); + auto factoryEPSG = AuthorityFactory::create(dbContext, "EPSG"); + + // NAD27 to WGS 84 (79) + auto op = factoryEPSG->createCoordinateOperation( + "15851", /* usePROJAlternativeGridNames = */ false); + auto transf = nn_dynamic_pointer_cast(op); + + auto crs = BoundCRS::create(GeographicCRS::EPSG_4267, // NAD27 + GeographicCRS::EPSG_4326, // WGS84 + NN_CHECK_THROW(transf)); + + EXPECT_EQ(crs->exportToPROJString( + PROJStringFormatter::create( + PROJStringFormatter::Convention::PROJ_4, dbContext) + .get()), + "+proj=longlat +ellps=clrk66 +nadgrids=us_noaa_conus.tif " + "+no_defs +type=crs"); +} + +// --------------------------------------------------------------------------- + +TEST(crs, boundCRS_nadcon_to_PROJ_string_usePROJAlternativeGridNames_true) { + + auto dbContext = DatabaseContext::create(); + auto factoryEPSG = AuthorityFactory::create(dbContext, "EPSG"); + + // NAD27 to WGS 84 (79) + auto op = factoryEPSG->createCoordinateOperation( + "15851", /* usePROJAlternativeGridNames = */ true); + auto transf = nn_dynamic_pointer_cast(op); + + auto crs = BoundCRS::create(GeographicCRS::EPSG_4267, // NAD27 + GeographicCRS::EPSG_4326, // WGS84 + NN_CHECK_THROW(transf)); + + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + "+proj=longlat +ellps=clrk66 +nadgrids=us_noaa_conus.tif +no_defs " + "+type=crs"); +} + +// --------------------------------------------------------------------------- + TEST(crs, boundCRS_identify_db) { auto dbContext = DatabaseContext::create(); auto factoryEPSG = AuthorityFactory::create(dbContext, "EPSG");