Merge pull request #1369 from rouault/grid_fixes

Grid related fixes

3 files changed, 68 insertions, 3 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index d7f138a4..50a6a003 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -8778,13 +8778,25 @@ void Transformation::_exportToPROJString(
         if (fileParameter &&
             fileParameter->type() == ParameterValue::Type::FILENAME) {
             auto filename = fileParameter->valueFile();
-            if (isMethodInverseOf) {
+            bool doInversion = isMethodInverseOf;
+            if (!identifiers().empty() &&
+                *identifiers().front()->codeSpace() ==
+                    metadata::Identifier::EPSG &&
+                method()->nameStr() ==
+                    "Geographic3D to GravityRelatedHeight (US .gtx)" &&
+                ends_with(filename, ".gtx")) {
+                // gtx files, from straight EPSG definition, must be applied in
+                // reverse order for "Geographic3D to GravityRelatedHeight"
+                // method
+                doInversion = !doInversion;
+            }
+            if (doInversion) {
                 formatter->startInversion();
             }
             formatter->addStep("vgridshift");
             formatter->addParam("grids", filename);
             formatter->addParam("multiplier", 1.0);
-            if (isMethodInverseOf) {
+            if (doInversion) {
                 formatter->stopInversion();
             }
             return;
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index b893faff..14421022 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -1024,8 +1024,8 @@ bool DatabaseContext::lookForGridInfo(const std::string &projFilename,
         info.url = url;
         info.directDownload = directDownload;
         info.openLicense = openLicense;
-        info.gridAvailable = gridAvailable;
     }
+    info.gridAvailable = gridAvailable;
     info.found = ret;
     d->cache(projFilename, info);
     return ret;
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 848fb23b..0f719a6d 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -5758,6 +5758,59 @@ TEST(operation, transformation_height_to_PROJ_string) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, transformation_Geographic3D_to_GravityRelatedHeight_gtx) {
+    auto wkt =
+        "COORDINATEOPERATION[\"ETRS89 to NAP height (1)\",\n"
+        "    VERSION[\"RDNAP-Nld 2008\"],\n"
+        "    SOURCECRS[\n"
+        "        GEOGCRS[\"ETRS89\",\n"
+        "            DATUM[\"European Terrestrial Reference System 1989\",\n"
+        "                ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
+        "                    LENGTHUNIT[\"metre\",1]]],\n"
+        "            PRIMEM[\"Greenwich\",0,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "            CS[ellipsoidal,3],\n"
+        "                AXIS[\"geodetic latitude (Lat)\",north,\n"
+        "                    ORDER[1],\n"
+        "                    ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "                AXIS[\"geodetic longitude (Lon)\",east,\n"
+        "                    ORDER[2],\n"
+        "                    ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "                AXIS[\"ellipsoidal height (h)\",up,\n"
+        "                    ORDER[3],\n"
+        "                    LENGTHUNIT[\"metre\",1]],\n"
+        "            ID[\"EPSG\",4937]]],\n"
+        "    TARGETCRS[\n"
+        "        VERTCRS[\"NAP height\",\n"
+        "            VDATUM[\"Normaal Amsterdams Peil\"],\n"
+        "            CS[vertical,1],\n"
+        "                AXIS[\"gravity-related height (H)\",up,\n"
+        "                    LENGTHUNIT[\"metre\",1]],\n"
+        "            ID[\"EPSG\",5709]]],\n"
+        "    METHOD[\"Geographic3D to GravityRelatedHeight (US .gtx)\",\n"
+        "        ID[\"EPSG\",9665]],\n"
+        "    PARAMETERFILE[\"Geoid (height correction) model "
+        "file\",\"naptrans2008.gtx\"],\n"
+        "    OPERATIONACCURACY[0.01],\n"
+        "    USAGE[\n"
+        "        SCOPE[\"unknown\"],\n"
+        "        AREA[\"Netherlands - onshore\"],\n"
+        "        BBOX[50.75,3.2,53.7,7.22]],\n"
+        "    ID[\"EPSG\",7001]]";
+    ;
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
+    ASSERT_TRUE(transf != nullptr);
+
+    // Check that we correctly inverse files in the case of
+    // "Geographic3D to GravityRelatedHeight (US .gtx)"
+    EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +inv +proj=vgridshift "
+              "+grids=naptrans2008.gtx +multiplier=1");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, transformation_ntv2_to_PROJ_string) {
     auto transformation = Transformation::createNTv2(
         PropertyMap(), GeographicCRS::EPSG_4807, GeographicCRS::EPSG_4326,