Fix import of WKT of concatenated operation with inverse conversion of a compound CRS of a projected CRS (fixes #3076)

2 files changed, 213 insertions, 6 deletions

diff --git a/src/iso19111/operation/concatenatedoperation.cpp b/src/iso19111/operation/concatenatedoperation.cpp
index 7da561b4..e5728c4c 100644
--- a/src/iso19111/operation/concatenatedoperation.cpp
+++ b/src/iso19111/operation/concatenatedoperation.cpp
@@ -276,14 +276,28 @@ void ConcatenatedOperation::fixStepsDirection(
         }
     }
 
+    const auto extractDerivedCRS =
+        [](const crs::CRS *crs) -> const crs::DerivedCRS * {
+        auto derivedCRS = dynamic_cast<const crs::DerivedCRS *>(crs);
+        if (derivedCRS)
+            return derivedCRS;
+        auto compoundCRS = dynamic_cast<const crs::CompoundCRS *>(crs);
+        if (compoundCRS) {
+            derivedCRS = dynamic_cast<const crs::DerivedCRS *>(
+                compoundCRS->componentReferenceSystems().front().get());
+            if (derivedCRS)
+                return derivedCRS;
+        }
+        return nullptr;
+    };
+
     for (size_t i = 0; i < operationsInOut.size(); ++i) {
         auto &op = operationsInOut[i];
         auto l_sourceCRS = op->sourceCRS();
         auto l_targetCRS = op->targetCRS();
         auto conv = dynamic_cast<const Conversion *>(op.get());
         if (conv && i == 0 && !l_sourceCRS && !l_targetCRS) {
-            if (auto derivedCRS = dynamic_cast<const crs::DerivedCRS *>(
-                    concatOpSourceCRS.get())) {
+            if (auto derivedCRS = extractDerivedCRS(concatOpSourceCRS.get())) {
                 if (i + 1 < operationsInOut.size()) {
                     // use the sourceCRS of the next operation as our target CRS
                     l_targetCRS = operationsInOut[i + 1]->sourceCRS();
@@ -323,8 +337,7 @@ void ConcatenatedOperation::fixStepsDirection(
             }
         } else if (conv && i + 1 == operationsInOut.size() && !l_sourceCRS &&
                    !l_targetCRS) {
-            auto derivedCRS =
-                dynamic_cast<const crs::DerivedCRS *>(concatOpTargetCRS.get());
+            auto derivedCRS = extractDerivedCRS(concatOpTargetCRS.get());
             if (derivedCRS) {
                 if (i >= 1) {
                     // use the sourceCRS of the previous operation as our source
@@ -350,8 +363,7 @@ void ConcatenatedOperation::fixStepsDirection(
             } else if (i >= 1) {
                 l_sourceCRS = operationsInOut[i - 1]->targetCRS();
                 if (l_sourceCRS) {
-                    derivedCRS = dynamic_cast<const crs::DerivedCRS *>(
-                        l_sourceCRS.get());
+                    derivedCRS = extractDerivedCRS(l_sourceCRS.get());
                     if (derivedCRS &&
                         conv->isEquivalentTo(
                             derivedCRS->derivingConversion().get(),
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 69cab65c..f883ddb7 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -3827,6 +3827,201 @@ TEST(
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, CONCATENATEDOPERATION_with_inverse_conversion_of_compound) {
+
+    auto wkt =
+        "CONCATENATEDOPERATION[\"Inverse of RD New + Amersfoort to ETRS89 (9) "
+        "+ Inverse of ETRS89 to NAP height (2) + ETRS89 to WGS 84 (1)\",\n"
+        "    SOURCECRS[\n"
+        "        COMPOUNDCRS[\"Amersfoort / RD New + NAP height\",\n"
+        "            PROJCRS[\"Amersfoort / RD New\",\n"
+        "                BASEGEOGCRS[\"Amersfoort\",\n"
+        "                    DATUM[\"Amersfoort\",\n"
+        "                        ELLIPSOID[\"Bessel "
+        "1841\",6377397.155,299.1528128,\n"
+        "                            LENGTHUNIT[\"metre\",1]]],\n"
+        "                    PRIMEM[\"Greenwich\",0,\n"
+        "                        ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "                    ID[\"EPSG\",4289]],\n"
+        "                CONVERSION[\"RD New\",\n"
+        "                    METHOD[\"Oblique Stereographic\",\n"
+        "                        ID[\"EPSG\",9809]],\n"
+        "                    PARAMETER[\"Latitude of natural "
+        "origin\",52.1561605555556,\n"
+        "                        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                        ID[\"EPSG\",8801]],\n"
+        "                    PARAMETER[\"Longitude of natural "
+        "origin\",5.38763888888889,\n"
+        "                        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                        ID[\"EPSG\",8802]],\n"
+        "                    PARAMETER[\"Scale factor at natural "
+        "origin\",0.9999079,\n"
+        "                        SCALEUNIT[\"unity\",1],\n"
+        "                        ID[\"EPSG\",8805]],\n"
+        "                    PARAMETER[\"False easting\",155000,\n"
+        "                        LENGTHUNIT[\"metre\",1],\n"
+        "                        ID[\"EPSG\",8806]],\n"
+        "                    PARAMETER[\"False northing\",463000,\n"
+        "                        LENGTHUNIT[\"metre\",1],\n"
+        "                        ID[\"EPSG\",8807]]],\n"
+        "                CS[Cartesian,2],\n"
+        "                    AXIS[\"easting (X)\",east,\n"
+        "                        ORDER[1],\n"
+        "                        LENGTHUNIT[\"metre\",1]],\n"
+        "                    AXIS[\"northing (Y)\",north,\n"
+        "                        ORDER[2],\n"
+        "                        LENGTHUNIT[\"metre\",1]]],\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\",7415]]],\n"
+        "    TARGETCRS[\n"
+        "        GEOGCRS[\"WGS 84 (3D)\",\n"
+        "            ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n"
+        "                MEMBER[\"World Geodetic System 1984 (Transit)\"],\n"
+        "                MEMBER[\"World Geodetic System 1984 (G730)\"],\n"
+        "                MEMBER[\"World Geodetic System 1984 (G873)\"],\n"
+        "                MEMBER[\"World Geodetic System 1984 (G1150)\"],\n"
+        "                MEMBER[\"World Geodetic System 1984 (G1674)\"],\n"
+        "                MEMBER[\"World Geodetic System 1984 (G1762)\"],\n"
+        "                MEMBER[\"World Geodetic System 1984 (G2139)\"],\n"
+        "                ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                    LENGTHUNIT[\"metre\",1]],\n"
+        "                ENSEMBLEACCURACY[2.0]],\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 minute second "
+        "hemisphere\",0.0174532925199433]],\n"
+        "                AXIS[\"geodetic longitude (Long)\",east,\n"
+        "                    ORDER[2],\n"
+        "                    ANGLEUNIT[\"degree minute second "
+        "hemisphere\",0.0174532925199433]],\n"
+        "                AXIS[\"ellipsoidal height (h)\",up,\n"
+        "                    ORDER[3],\n"
+        "                    LENGTHUNIT[\"metre\",1]],\n"
+        "            ID[\"EPSG\",4329]]],\n"
+        "    STEP[\n"
+        "        CONVERSION[\"Inverse of RD New\",\n"
+        "            METHOD[\"Inverse of Oblique Stereographic\",\n"
+        "                ID[\"INVERSE(EPSG)\",9809]],\n"
+        "            PARAMETER[\"Latitude of natural "
+        "origin\",52.1561605555556,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8801]],\n"
+        "            PARAMETER[\"Longitude of natural "
+        "origin\",5.38763888888889,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8802]],\n"
+        "            PARAMETER[\"Scale factor at natural origin\",0.9999079,\n"
+        "                SCALEUNIT[\"unity\",1],\n"
+        "                ID[\"EPSG\",8805]],\n"
+        "            PARAMETER[\"False easting\",155000,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8806]],\n"
+        "            PARAMETER[\"False northing\",463000,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8807]],\n"
+        "            ID[\"INVERSE(EPSG)\",19914]]],\n"
+        "    STEP[\n"
+        "        COORDINATEOPERATION[\"PROJ-based coordinate operation\",\n"
+        "            SOURCECRS[\n"
+        "                COMPOUNDCRS[\"Amersfoort + NAP height\",\n"
+        "                    GEOGCRS[\"Amersfoort\",\n"
+        "                        DATUM[\"Amersfoort\",\n"
+        "                            ELLIPSOID[\"Bessel "
+        "1841\",6377397.155,299.1528128,\n"
+        "                                LENGTHUNIT[\"metre\",1]]],\n"
+        "                        PRIMEM[\"Greenwich\",0,\n"
+        "                            "
+        "ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "                        CS[ellipsoidal,2],\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"
+        "                        ID[\"EPSG\",4289]],\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"
+        "            TARGETCRS[\n"
+        "                GEOGCRS[\"WGS 84 (3D)\",\n"
+        "                    ENSEMBLE[\"World Geodetic System 1984 "
+        "ensemble\",\n"
+        "                        MEMBER[\"World Geodetic System 1984 "
+        "(Transit)\"],\n"
+        "                        MEMBER[\"World Geodetic System 1984 "
+        "(G730)\"],\n"
+        "                        MEMBER[\"World Geodetic System 1984 "
+        "(G873)\"],\n"
+        "                        MEMBER[\"World Geodetic System 1984 "
+        "(G1150)\"],\n"
+        "                        MEMBER[\"World Geodetic System 1984 "
+        "(G1674)\"],\n"
+        "                        MEMBER[\"World Geodetic System 1984 "
+        "(G1762)\"],\n"
+        "                        MEMBER[\"World Geodetic System 1984 "
+        "(G2139)\"],\n"
+        "                        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                            LENGTHUNIT[\"metre\",1]],\n"
+        "                        ENSEMBLEACCURACY[2.0]],\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 minute second "
+        "hemisphere\",0.0174532925199433]],\n"
+        "                        AXIS[\"geodetic longitude (Long)\",east,\n"
+        "                            ORDER[2],\n"
+        "                            ANGLEUNIT[\"degree minute second "
+        "hemisphere\",0.0174532925199433]],\n"
+        "                        AXIS[\"ellipsoidal height (h)\",up,\n"
+        "                            ORDER[3],\n"
+        "                            LENGTHUNIT[\"metre\",1]],\n"
+        "                    ID[\"EPSG\",4329]]],\n"
+        "            METHOD[\"PROJ-based operation method: +proj=pipeline "
+        "+step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg "
+        "+xy_out=rad +step +proj=hgridshift +grids=nl_nsgi_rdtrans2018.tif "
+        "+step +proj=vgridshift +grids=nl_nsgi_nlgeo2018.tif +multiplier=1 "
+        "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+        "+order=2,1\"],\n"
+        "            OPERATIONACCURACY[1.002]]],\n"
+        "    USAGE[\n"
+        "        SCOPE[\"unknown\"],\n"
+        "        AREA[\"Netherlands - onshore, including Waddenzee, Dutch "
+        "Wadden Islands and 12-mile offshore coastal zone.\"],\n"
+        "        BBOX[50.75,3.2,53.7,7.22]]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto concat = nn_dynamic_pointer_cast<ConcatenatedOperation>(obj);
+    ASSERT_TRUE(concat != nullptr);
+
+    EXPECT_EQ(concat->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +inv +proj=sterea +lat_0=52.1561605555556 "
+              "+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 "
+              "+ellps=bessel "
+              "+step +proj=hgridshift +grids=nl_nsgi_rdtrans2018.tif "
+              "+step +proj=vgridshift +grids=nl_nsgi_nlgeo2018.tif "
+              "+multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, BOUNDCRS_transformation_from_names) {
 
     auto projcrs = ProjectedCRS::create(