Merge pull request #4168 from rouault/fix_4154

BoundCRS::exportToPROJ(): handle case of NADCON conus grid
OSGeo · Jun 10, 2024 · 2ebf62c · 2ebf62c
2 parents d467ecb + 926bbd2
commit 2ebf62c
Show file tree

Hide file tree

Showing 5 changed files with 97 additions and 8 deletions.
diff --git 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<double>
     getTOWGS84Parameters() const; // throw(io::FormattingException)

diff --git 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
@@ -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
@@ -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
@@ -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<Transformation>(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<Transformation>(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");