ConcatenatedOperation::fixStepsDirection(): add even more hacks to de…

…al with EPSG:10675 'BES2020 to Saba height'
rouault · Dec 17, 2024 · 4064893 · 4064893
1 parent 40170a8
commit 4064893
Show file tree

Hide file tree

Showing 7 changed files with 128 additions and 42 deletions.
diff --git a/data/CMakeLists.txt b/data/CMakeLists.txt
@@ -45,7 +45,7 @@ set(ALL_SQL_IN "${CMAKE_CURRENT_BINARY_DIR}/all.sql.in")
 set(PROJ_DB "${CMAKE_CURRENT_BINARY_DIR}/proj.db")
 include(sql_filelist.cmake)
 
-set(PROJ_DB_SQL_EXPECTED_MD5 "c2184d9fca631e88f21e4c0d8899d391")
+set(PROJ_DB_SQL_EXPECTED_MD5 "99e510ac698e4961c38cc2a09d07ccb0")
 
 add_custom_command(
   OUTPUT ${PROJ_DB}

diff --git a/data/sql/concatenated_operation.sql b/data/sql/concatenated_operation.sql
@@ -495,3 +495,5 @@ INSERT INTO "concatenated_operation" VALUES('EPSG','10496','ETRS89 + DVR90(2013)
 INSERT INTO "usage" VALUES('EPSG','20516','concatenated_operation','EPSG','10496','EPSG','1080','EPSG','1273');
 INSERT INTO "concatenated_operation" VALUES('EPSG','10616','SRGI2013 + INAGeoid2020 v1 height to SRGI2013 + INAGeoid v2 height (1)','In central Java INAGeoid2020 v2 height minus INAGeoid2020 v1 height is approximately +0.2m (v1 surface is above the v2 surface). This difference varies significantly across Indonesia.','EPSG','9529','EPSG','20043',0.2,'BIG-Idn 2022',0);
 INSERT INTO "usage" VALUES('EPSG','21321','concatenated_operation','EPSG','10616','EPSG','1122','EPSG','1178');
+INSERT INTO "concatenated_operation" VALUES('EPSG','10675','BES2020 to Saba height (1)','Although documented with source and target CRSs in the geog2D domain, step 1 is a 3D transformation in which implicit concatenated operations through geocentric Cartesian coordinates should be used. Refer to EPSG Guidance Notes 7-1 and 7-2.','EPSG','10638','EPSG','10642',0.1,'NSDI-Bes 2023',0);
+INSERT INTO "usage" VALUES('EPSG','21877','concatenated_operation','EPSG','10675','EPSG','4757','EPSG','1133');
diff --git a/data/sql/concatenated_operation_step.sql b/data/sql/concatenated_operation_step.sql
@@ -504,3 +504,5 @@ INSERT INTO "concatenated_operation_step" VALUES('EPSG','10496',1,'EPSG','10492'
 INSERT INTO "concatenated_operation_step" VALUES('EPSG','10496',2,'EPSG','10494');
 INSERT INTO "concatenated_operation_step" VALUES('EPSG','10616',1,'EPSG','9629');
 INSERT INTO "concatenated_operation_step" VALUES('EPSG','10616',2,'EPSG','10145');
+INSERT INTO "concatenated_operation_step" VALUES('EPSG','10675',1,'EPSG','10646');
+INSERT INTO "concatenated_operation_step" VALUES('EPSG','10675',2,'EPSG','10657');
diff --git a/data/sql/consistency_checks_triggers.sql b/data/sql/consistency_checks_triggers.sql
@@ -611,6 +611,25 @@ FOR EACH ROW BEGIN
             AND concat_op.source_crs_auth_name = step_op.target_crs_auth_name
             AND concat_op.source_crs_code = step_op.target_crs_code)
 
+        -- same as above, but check by CRS name, and only for geodetic CRS.
+        -- For example the concatenated operation EPSG:10675 ("BES2020 to Saba height (1)")
+        -- has EPSG:10638 "BES2020" (geographic 3D) as source CRS
+        -- but its first step is EPSG:10646 ("Saba to BES2020 Saba (1)") which has
+        -- EPSG:10639 "BES2020" (geographic 2D) as target CRS !
+        AND NOT EXISTS (
+            SELECT 1 FROM coordinate_operation_view step_op
+            LEFT JOIN concatenated_operation concat_op ON
+            concat_op.auth_name = NEW.operation_auth_name AND concat_op.code = NEW.operation_code
+            LEFT JOIN geodetic_crs concat_op_source_crs ON
+            concat_op_source_crs.auth_name = concat_op.source_crs_auth_name
+            AND concat_op_source_crs.code = concat_op.source_crs_code
+            LEFT JOIN geodetic_crs step_op_target_crs ON
+            step_op_target_crs.auth_name = step_op.target_crs_auth_name
+            AND step_op_target_crs.code = step_op.target_crs_code
+            WHERE concat_op.deprecated = 0
+            AND step_op.auth_name = NEW.step_auth_name AND step_op.code = NEW.step_code
+            AND concat_op_source_crs.name = step_op_target_crs.name)
+
         -- or if source_crs of step 1 is a conversion
         AND NOT EXISTS (
             SELECT 1 FROM conversion_table step_op

diff --git a/scripts/build_db.py b/scripts/build_db.py
@@ -930,13 +930,6 @@ def fill_concatenated_operation(proj_db_cursor):
         expected_order = 1
         steps_code = []
 
-        # FIXME: https://epsg.org/concatenated-operation_10675/BES2020-to-Saba-height-1.html ill defined in EPSG 11.023
-        # due to first step referencing BES2020 Saba geographic 2D (EPSG:10639), but source CRS of concatenated
-        # operation referencing BES2020 Saba geographic 3D (EPSG:10638)
-        if code == 10675:
-            print("FIXME! Skipping EPSG:10675 'BES2020 to Saba height (1)' for now")
-            continue
-
         iterator = proj_db_cursor.execute("SELECT op_path_step, single_operation_code FROM epsg_coordoperationpath WHERE concat_operation_code = ? ORDER BY op_path_step", (code,))
         for (order, single_operation_code) in iterator:
             assert order == expected_order

diff --git a/src/iso19111/operation/concatenatedoperation.cpp b/src/iso19111/operation/concatenatedoperation.cpp
@@ -569,36 +569,73 @@ void ConcatenatedOperation::fixStepsDirection(
                 auto newOp(Conversion::createGeographicGeocentric(
                     NN_NO_CHECK(prevOpTarget), NN_NO_CHECK(l_targetCRS)));
                 operationsInOut.insert(operationsInOut.begin() + i, newOp);
-                // Particular case for https://github.com/OSGeo/PROJ/issues/3819
-                // where the antepenultimate transformation goes to A
-                // (geographic 3D) // and the last transformation is a NADCON 3D
-                // but between A (geographic 2D) to B (geographic 2D), and the
-                // concatenated transformation target CRS is B (geographic 3D)
-                // This is due to an oddity of the EPSG database that registers
-                // the NADCON 3D transformation between the 2D geographic CRS
-                // and not the 3D ones.
-            } else if (i + 1 == operationsInOut.size() &&
-                       l_sourceCRS->nameStr() == prevOpTarget->nameStr() &&
-                       l_targetCRS->nameStr() == concatOpTargetCRS->nameStr() &&
-                       isGeographic(l_targetCRS.get()) &&
-                       isGeographic(concatOpTargetCRS.get()) &&
-                       isGeographic(l_sourceCRS.get()) &&
-                       isGeographic(prevOpTarget.get()) &&
-                       dynamic_cast<const crs::GeographicCRS *>(
-                           prevOpTarget.get())
-                               ->coordinateSystem()
-                               ->axisList()
-                               .size() == 3 &&
-                       dynamic_cast<const crs::GeographicCRS *>(
-                           l_sourceCRS.get())
-                               ->coordinateSystem()
-                               ->axisList()
-                               .size() == 2 &&
-                       dynamic_cast<const crs::GeographicCRS *>(
-                           l_targetCRS.get())
-                               ->coordinateSystem()
-                               ->axisList()
-                               .size() == 2) {
+            } else if (i == 0 &&
+                       concatOpSourceCRS->nameStr() == l_targetCRS->nameStr() &&
+                       isGeographic(concatOpSourceCRS.get()) &&
+                       isGeographic(l_sourceCRS.get())) {
+                // Particular case for EPSG:10675 "BES2020 to Saba height"
+                op = op->inverse();
+
+                // This logic to convert between Geog2D <--> Geog3D could
+                // potentiall by generalized...
+                operationsInOut.insert(
+                    operationsInOut.begin(),
+                    NN_CHECK_ASSERT(
+                        CoordinateOperationFactory::create()->createOperation(
+                            concatOpSourceCRS, NN_NO_CHECK(l_targetCRS))));
+
+                // This logic to convert between Geog2D <--> Geog3D could
+                // potentiall by generalized...
+                if (operationsInOut.size() >= 3 &&
+                    operationsInOut[1]->targetCRS() &&
+                    operationsInOut[2]->sourceCRS() &&
+                    !areCRSMoreOrLessEquivalent(
+                        operationsInOut[1]->targetCRS().get(),
+                        operationsInOut[2]->sourceCRS().get()) &&
+                    operationsInOut[1]->targetCRS()->nameStr() ==
+                        operationsInOut[2]->sourceCRS()->nameStr() &&
+                    isGeographic(operationsInOut[1]->targetCRS().get()) &&
+                    isGeographic(operationsInOut[2]->sourceCRS().get())) {
+                    operationsInOut.insert(
+                        operationsInOut.begin() + 2,
+                        NN_CHECK_ASSERT(
+                            CoordinateOperationFactory::create()
+                                ->createOperation(
+                                    NN_NO_CHECK(
+                                        operationsInOut[1]->targetCRS()),
+                                    NN_NO_CHECK(
+                                        operationsInOut[2]->sourceCRS()))));
+                }
+            }
+            // Particular case for https://github.com/OSGeo/PROJ/issues/3819
+            // where the antepenultimate transformation goes to A
+            // (geographic 3D)
+            // and the last transformation is a NADCON 3D
+            // but between A (geographic 2D) to B (geographic 2D), and the
+            // concatenated transformation target CRS is B (geographic 3D)
+            // This is due to an oddity of the EPSG database that registers
+            // the NADCON 3D transformation between the 2D geographic CRS
+            // and not the 3D ones.
+            else if (i + 1 == operationsInOut.size() &&
+                     l_sourceCRS->nameStr() == prevOpTarget->nameStr() &&
+                     l_targetCRS->nameStr() == concatOpTargetCRS->nameStr() &&
+                     isGeographic(l_targetCRS.get()) &&
+                     isGeographic(concatOpTargetCRS.get()) &&
+                     isGeographic(l_sourceCRS.get()) &&
+                     isGeographic(prevOpTarget.get()) &&
+                     dynamic_cast<const crs::GeographicCRS *>(
+                         prevOpTarget.get())
+                             ->coordinateSystem()
+                             ->axisList()
+                             .size() == 3 &&
+                     dynamic_cast<const crs::GeographicCRS *>(l_sourceCRS.get())
+                             ->coordinateSystem()
+                             ->axisList()
+                             .size() == 2 &&
+                     dynamic_cast<const crs::GeographicCRS *>(l_targetCRS.get())
+                             ->coordinateSystem()
+                             ->axisList()
+                             .size() == 2) {
                 const auto transf =
                     dynamic_cast<const Transformation *>(op.get());
                 if (transf &&
@@ -618,10 +655,20 @@ void ConcatenatedOperation::fixStepsDirection(
         auto l_sourceCRS = operationsInOut.front()->sourceCRS();
         if (l_sourceCRS && !areCRSMoreOrLessEquivalent(
                                l_sourceCRS.get(), concatOpSourceCRS.get())) {
-            throw InvalidOperation("The source CRS of the first step of "
-                                   "concatenated operation is not the same "
-                                   "as the source CRS of the concatenated "
-                                   "operation itself");
+            if (l_sourceCRS->nameStr() == concatOpSourceCRS->nameStr() &&
+                ((isGeographic(l_sourceCRS.get()) &&
+                  isGeocentric(concatOpSourceCRS.get())) ||
+                 (isGeocentric(l_sourceCRS.get()) &&
+                  isGeographic(concatOpSourceCRS.get())))) {
+                auto newOp(Conversion::createGeographicGeocentric(
+                    NN_NO_CHECK(l_sourceCRS), concatOpSourceCRS));
+                operationsInOut.push_back(newOp);
+            } else {
+                throw InvalidOperation("The source CRS of the first step of "
+                                       "concatenated operation is not the same "
+                                       "as the source CRS of the concatenated "
+                                       "operation itself");
+            }
         }
 
         auto l_targetCRS = operationsInOut.back()->targetCRS();

diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
@@ -6164,3 +6164,26 @@ TEST(operation, CoordinateFrameRotationFullMatrixGeog2D) {
               "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
               "+step +proj=axisswap +order=2,1");
 }
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, ConcatenatedOperationEPSG_10675) {
+    auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto op = factory->createCoordinateOperation("10675", false);
+    auto concatOp = nn_dynamic_pointer_cast<ConcatenatedOperation>(op);
+    ASSERT_TRUE(concatOp != nullptr);
+    EXPECT_EQ(concatOp->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=push +v_3 "
+              "+step +proj=cart +ellps=GRS80 "
+              "+step +inv +proj=helmert +exact "
+              "+x=1138.7432 +y=-2064.4761 +z=110.7016 "
+              "+rx=-214.615206 +ry=479.360036 +rz=-164.703951 +s=-402.32073 "
+              "+convention=coordinate_frame "
+              "+step +inv +proj=cart +ellps=intl "
+              "+step +proj=pop +v_3 "
+              "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+              "+step +proj=axisswap +order=2,1");
+}