Merge pull request #3333 from rouault/fix_3328

Fix issue related to transforming to CGVD28 height (EPSG:5713)

data/sql/commit.sql

@@ -73,14 +73,14 @@ FOR EACH ROW BEGIN
                       AND g1.source_crs_code = g2.source_crs_code
                       AND g1.target_crs_auth_name = g2.target_crs_auth_name
                       AND g1.target_crs_code = g2.target_crs_code
-                      WHERE g1.auth_name = 'PROJ' AND g1.code NOT LIKE '%_RESTRICTED_TO_VERTCRS%' AND g2.auth_name = 'EPSG')
+                      WHERE g1.auth_name = 'PROJ' AND g1.code NOT LIKE '%_RESTRICTED_TO_VERTCRS%' AND g2.auth_name = 'EPSG' AND g2.deprecated = 0)
         OR EXISTS (SELECT 1 FROM grid_transformation g1
                       JOIN grid_transformation g2
                       ON g1.source_crs_auth_name = g2.target_crs_auth_name
                       AND g1.source_crs_code = g2.target_crs_code
                       AND g1.target_crs_auth_name = g1.source_crs_auth_name
                       AND g1.target_crs_code = g1.source_crs_code
-                      WHERE g1.auth_name = 'PROJ' AND g1.code NOT LIKE '%_RESTRICTED_TO_VERTCRS%' AND g2.auth_name = 'EPSG');
+                      WHERE g1.auth_name = 'PROJ' AND g1.code NOT LIKE '%_RESTRICTED_TO_VERTCRS%' AND g2.auth_name = 'EPSG' AND g2.deprecated = 0);
 
     SELECT RAISE(ABORT, 'Arg! there is now a EPSG:102100 object. Hack in createFromUserInput() will no longer work')
         WHERE EXISTS(SELECT 1 FROM crs_view WHERE auth_name = 'EPSG' AND code = '102100');

data/sql/customizations.sql

@@ -144,6 +144,50 @@ UPDATE grid_transformation SET accuracy = 2.0 WHERE auth_name = 'EPSG' AND code
 
 UPDATE grid_transformation SET accuracy = 2.0 WHERE auth_name = 'EPSG' AND code = '1462';
 
+
+-- Create a PROJ copy of EPSG:9123 "NAD83(CSRS) to CGVD28 height (1)" by
+-- removing the deprecation flag.
+-- Having a transformation from/to CGVD28 using the "generic" NAD83(CSRS) is much
+-- more convenient for low accuracy cases like /~https://github.com/OSGeo/PROJ/issues/3328
+INSERT INTO grid_transformation SELECT
+    'PROJ' AS auth_name,
+    'EPSG_9123' AS code,
+    name,
+    description || ' Imported from EPSG:9123 with deprecation flag removed by PROJ' AS description,
+    method_auth_name,
+    method_code,
+    method_name,
+    source_crs_auth_name,
+    source_crs_code,
+    target_crs_auth_name,
+    target_crs_code,
+    accuracy,
+    grid_param_auth_name,
+    grid_param_code,
+    grid_param_name,
+    grid_name,
+    grid2_param_auth_name,
+    grid2_param_code,
+    grid2_param_name,
+    grid2_name,
+    interpolation_crs_auth_name,
+    interpolation_crs_code,
+    operation_version,
+    0 AS deprecated
+    FROM grid_transformation WHERE auth_name = 'EPSG' AND code = '9123';
+
+INSERT INTO usage SELECT
+    'PROJ' AS auth_name,
+    'USAGE_PROJ_EPSG_9123' AS code,
+    object_table_name,
+    'PROJ' AS object_auth_name,
+    'EPSG_9123' AS object_code,
+    extent_auth_name,
+    extent_code,
+    scope_auth_name,
+    scope_code
+    FROM usage WHERE object_table_name = 'grid_transformation' AND object_auth_name = 'EPSG' AND object_code = '9123';
+
 -- Define the allowed authorities, and their precedence, when researching a
 -- coordinate operation
 

data/sql/grid_alternatives.sql

@@ -47,7 +47,8 @@ VALUES
 ('CGG2013an83.byn','ca_nrc_CGG2013an83.tif','CGG2013an83.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_CGG2013an83.tif',1,1,NULL),
 ('CGG2013i83.byn','ca_nrc_CGG2013i08.tif','CGG2013i08.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_CGG2013i08.tif',1,1,NULL),
 ('CGG2013n83.byn','ca_nrc_CGG2013n83.tif','CGG2013n83.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_CGG2013n83.tif',1,1,NULL),
-('HT2_0.byn','ca_nrc_HT2_2010v70.tif','HT2_2010v70.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70.tif',1,1,NULL), -- deprecated name
+('HT2_0.byn','ca_nrc_HT2_2010v70.tif','HT2_2010v70.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70.tif',1,1,NULL), -- deprecated name. dubious mapping !!!
+('HT2_1997.byn','ca_nrc_HT2_2010v70.tif',NULL,'GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70.tif',1,1,NULL), -- dubious mapping !!!
 ('HT2_2010v70.byn','ca_nrc_HT2_2010v70.tif','HT2_2010v70.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70.tif',1,1,NULL),
 ('HT2_2010v70_CGG2013a.byn','ca_nrc_HT2_2010v70_CGG2013a.tif',NULL,'GTiff','vgridshift',0,NULL,'https://cdn.proj.org/ca_nrc_HT2_2010v70_CGG2013a.tif',1,1,NULL),
 -- the PROJ grid is the reverse way of the EPSG one

src/iso19111/io.cpp

@@ -8618,7 +8618,8 @@ const std::string &PROJStringFormatter::toString() const {
                         iterNext = std::next(iterNext);
                     }
                 }
-                if (ok && iterNext != steps.end()) {
+                ok &= iterNext != steps.end();
+                if (ok) {
                     ok = false;
                     auto &nextStep = *iterNext;
                     if (nextStep.name == "axisswap" &&
@@ -8628,7 +8629,8 @@ const std::string &PROJStringFormatter::toString() const {
                         iterNext = std::next(iterNext);
                     }
                 }
-                if (ok && iterNext != steps.end()) {
+                ok &= iterNext != steps.end();
+                if (ok) {
                     ok = false;
                     auto &nextStep = *iterNext;
                     if (nextStep.name == "unitconvert" &&
@@ -8640,15 +8642,17 @@ const std::string &PROJStringFormatter::toString() const {
                     }
                 }
                 auto iterVgridshift = iterNext;
-                if (ok && iterNext != steps.end()) {
+                ok &= iterNext != steps.end();
+                if (ok) {
                     ok = false;
                     auto &nextStep = *iterNext;
                     if (nextStep.name == "vgridshift") {
                         ok = true;
                         iterNext = std::next(iterNext);
                     }
                 }
-                if (ok && iterNext != steps.end()) {
+                ok &= iterNext != steps.end();
+                if (ok) {
                     ok = false;
                     auto &nextStep = *iterNext;
                     if (nextStep.name == "unitconvert" &&
@@ -8659,7 +8663,8 @@ const std::string &PROJStringFormatter::toString() const {
                         iterNext = std::next(iterNext);
                     }
                 }
-                if (ok && iterNext != steps.end()) {
+                ok &= iterNext != steps.end();
+                if (ok) {
                     ok = false;
                     auto &nextStep = *iterNext;
                     if (nextStep.name == "axisswap" &&
@@ -8669,14 +8674,16 @@ const std::string &PROJStringFormatter::toString() const {
                         iterNext = std::next(iterNext);
                     }
                 }
-                if (ok && iterNext != steps.end()) {
+                ok &= iterNext != steps.end();
+                if (ok) {
                     ok = false;
                     auto &nextStep = *iterNext;
                     if (nextStep.name == "pop" &&
                         nextStep.paramValues.size() == 2 &&
                         nextStep.paramValues[0].keyEquals("v_1") &&
                         nextStep.paramValues[1].keyEquals("v_2")) {
                         ok = true;
+                        // iterNext = std::next(iterNext);
                     }
                 }
                 if (ok) {
@@ -8688,6 +8695,107 @@ const std::string &PROJStringFormatter::toString() const {
                 }
             }
 
+            // +step +proj=axisswap +order=2,1
+            // +step +proj=unitconvert +xy_in=deg +xy_out=rad
+            // +step +proj=vgridshift ...
+            // +step +proj=unitconvert +xy_in=rad +xy_out=deg
+            // +step +proj=axisswap +order=2,1
+            // +step +proj=push +v_1 +v_2
+            // +step +proj=axisswap +order=2,1
+            // +step +proj=unitconvert +xy_in=deg +xy_out=rad
+            // ==>
+            // +step +proj=push +v_1 +v_2
+            // +step +proj=axisswap +order=2,1
+            // +step +proj=unitconvert +xy_in=deg +xy_out=rad
+            // +step +proj=vgridshift ...
+
+            if (prevStep.name == "axisswap" && prevStepParamCount == 1 &&
+                prevStep.paramValues[0].equals("order", "2,1") &&
+                curStep.name == "unitconvert" && curStepParamCount == 2 &&
+                !curStep.inverted &&
+                curStep.paramValues[0].equals("xy_in", "deg") &&
+                curStep.paramValues[1].equals("xy_out", "rad")) {
+                auto iterNext = std::next(iterCur);
+                bool ok = false;
+                auto iterVgridshift = iterNext;
+                if (iterNext != steps.end()) {
+                    auto &nextStep = *iterNext;
+                    if (nextStep.name == "vgridshift") {
+                        ok = true;
+                        iterNext = std::next(iterNext);
+                    }
+                }
+                ok &= iterNext != steps.end();
+                if (ok) {
+                    ok = false;
+                    auto &nextStep = *iterNext;
+                    if (nextStep.name == "unitconvert" && !nextStep.inverted &&
+                        nextStep.paramValues.size() == 2 &&
+                        nextStep.paramValues[0].equals("xy_in", "rad") &&
+                        nextStep.paramValues[1].equals("xy_out", "deg")) {
+                        ok = true;
+                        iterNext = std::next(iterNext);
+                    }
+                }
+                ok &= iterNext != steps.end();
+                if (ok) {
+                    ok = false;
+                    auto &nextStep = *iterNext;
+                    if (nextStep.name == "axisswap" &&
+                        nextStep.paramValues.size() == 1 &&
+                        nextStep.paramValues[0].equals("order", "2,1")) {
+                        ok = true;
+                        iterNext = std::next(iterNext);
+                    }
+                }
+                auto iterPush = iterNext;
+                ok &= iterNext != steps.end();
+                if (ok) {
+                    ok = false;
+                    auto &nextStep = *iterNext;
+                    if (nextStep.name == "push" &&
+                        nextStep.paramValues.size() == 2 &&
+                        nextStep.paramValues[0].keyEquals("v_1") &&
+                        nextStep.paramValues[1].keyEquals("v_2")) {
+                        ok = true;
+                        iterNext = std::next(iterNext);
+                    }
+                }
+                ok &= iterNext != steps.end();
+                if (ok) {
+                    ok = false;
+                    auto &nextStep = *iterNext;
+                    if (nextStep.name == "axisswap" &&
+                        nextStep.paramValues.size() == 1 &&
+                        nextStep.paramValues[0].equals("order", "2,1")) {
+                        ok = true;
+                        iterNext = std::next(iterNext);
+                    }
+                }
+                ok &= iterNext != steps.end();
+                if (ok) {
+                    ok = false;
+                    auto &nextStep = *iterNext;
+                    if (nextStep.name == "unitconvert" &&
+                        nextStep.paramValues.size() == 2 &&
+                        !nextStep.inverted &&
+                        nextStep.paramValues[0].equals("xy_in", "deg") &&
+                        nextStep.paramValues[1].equals("xy_out", "rad")) {
+                        ok = true;
+                        // iterNext = std::next(iterNext);
+                    }
+                }
+
+                if (ok) {
+                    auto stepVgridshift = *iterVgridshift;
+                    steps.erase(iterPrev, iterPush);
+                    steps.insert(std::next(iterNext), stepVgridshift);
+                    iterPrev = iterPush;
+                    iterCur = std::next(iterPush);
+                    continue;
+                }
+            }
+
             ++iterCur;
         }
     }

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -2099,7 +2099,17 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
                               ->demoteTo2D(std::string(), nullptr)
                               .get(),
                           util::IComparable::Criterion::EQUIVALENT)) {
-            const int methodEPSGCode = transf->method()->getEPSGCode();
+            int methodEPSGCode = transf->method()->getEPSGCode();
+            if (methodEPSGCode == 0) {
+                // If the transformation is actually an inverse transformation,
+                // we will not get the EPSG code. So get the forward
+                // transformation.
+                const auto invTrans = transf->inverse();
+                const auto invTransAsTrans =
+                    dynamic_cast<Transformation *>(invTrans.get());
+                if (invTransAsTrans)
+                    methodEPSGCode = invTransAsTrans->method()->getEPSGCode();
+            }
 
             const bool bGeocentricTranslation =
                 methodEPSGCode ==

test/unit/test_io.cpp

@@ -8942,6 +8942,49 @@ TEST(io, projstringformatter_optim_hgridshift_vgridshift_hgridshift_inv) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projstringformatter_optim_as_uc_vgridshift_uc_as_push_as_uc) {
+    // Nominal case
+    {
+        auto fmt = PROJStringFormatter::create();
+        fmt->addStep("axisswap");
+        fmt->addParam("order", "2,1");
+
+        fmt->addStep("unitconvert");
+        fmt->addParam("xy_in", "deg");
+        fmt->addParam("xy_out", "rad");
+
+        fmt->addStep("vgridshift");
+        fmt->addParam("grids", "foo");
+
+        fmt->addStep("unitconvert");
+        fmt->addParam("xy_in", "rad");
+        fmt->addParam("xy_out", "deg");
+
+        fmt->addStep("axisswap");
+        fmt->addParam("order", "2,1");
+
+        fmt->addStep("push");
+        fmt->addParam("v_1");
+        fmt->addParam("v_2");
+
+        fmt->addStep("axisswap");
+        fmt->addParam("order", "2,1");
+
+        fmt->addStep("unitconvert");
+        fmt->addParam("xy_in", "deg");
+        fmt->addParam("xy_out", "rad");
+
+        EXPECT_EQ(fmt->toString(),
+                  "+proj=pipeline "
+                  "+step +proj=push +v_1 +v_2 "
+                  "+step +proj=axisswap +order=2,1 "
+                  "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+                  "+step +proj=vgridshift +grids=foo");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_longlat) {
 
     auto expected = "GEODCRS[\"unknown\",\n"

test/unit/test_operationfactory.cpp

@@ -4991,6 +4991,63 @@ TEST(operation, compoundCRS_to_compoundCRS_issue_2720) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_to_compoundCRS_issue_3328) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), std::string());
+    // "WGS 84 + EGM96 height"
+    auto srcObj = createFromUserInput("EPSG:4326+5773",
+                                      authFactory->databaseContext(), false);
+    auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+    ASSERT_TRUE(src != nullptr);
+
+    // "WGS 84 + CGVD28 height"
+    auto dstObj = createFromUserInput("EPSG:4326+5713",
+                                      authFactory->databaseContext(), false);
+    auto dst = nn_dynamic_pointer_cast<CRS>(dstObj);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt);
+    ASSERT_GE(list.size(), 1U);
+    EXPECT_EQ(list[0]->nameStr(), "Inverse of WGS 84 to EGM96 height (1) + "
+                                  "Inverse of NAD83(CSRS) to WGS 84 (2) + "
+                                  "NAD83(CSRS) to CGVD28 height (1) + "
+                                  "NAD83(CSRS) to WGS 84 (2)");
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=push +v_1 +v_2 "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1 "
+              "+step +proj=cart +ellps=WGS84 "
+              "+step +inv +proj=helmert +x=-0.991 +y=1.9072 +z=0.5129 "
+              "+rx=-0.0257899075194932 "
+              "+ry=-0.0096500989602704 +rz=-0.0116599432323421 +s=0 "
+              "+convention=coordinate_frame "
+              "+step +inv +proj=cart +ellps=GRS80 "
+              "+step +inv +proj=vgridshift +grids=ca_nrc_HT2_2010v70.tif "
+              "+multiplier=1 "
+              "+step +proj=push +v_3 "
+              "+step +proj=cart +ellps=GRS80 "
+              "+step +proj=helmert +x=-0.991 +y=1.9072 +z=0.5129 "
+              "+rx=-0.0257899075194932 "
+              "+ry=-0.0096500989602704 +rz=-0.0116599432323421 +s=0 "
+              "+convention=coordinate_frame "
+              "+step +inv +proj=cart +ellps=WGS84 "
+              "+step +proj=pop +v_3 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=pop +v_1 +v_2");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(
     operation,
     compoundCRS_to_compoundCRS_concatenated_operation_with_two_vert_transformation_and_ballpark_geog) {