Merge pull request #4364 from rouault/etrf_to_etrf

Improve ETRFxxx to ETRFyyy, and WGS 84 (xxx) to WGS 84 (yyy)

Author: rouault

data/CMakeLists.txt

@@ -31,7 +31,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 "5e23b08c318cecbd9cace26233f8f5d4")
+set(PROJ_DB_SQL_EXPECTED_MD5 "8145b5734f4007f81b188ac3b690e8b3")
 
 add_custom_command(
   OUTPUT ${PROJ_DB}

data/sql/wgs84_realizations_concatenated_operations.sql

@@ -58,6 +58,7 @@ list(APPEND SQL_FILES
   "${SQL_DIR}/grid_alternatives.sql"
   "${SQL_DIR}/grid_alternatives_generated_noaa.sql"
   "${SQL_DIR}/nadcon5_concatenated_operations.sql"
+  "${SQL_DIR}/wgs84_realizations_concatenated_operations.sql"
   "${SQL_DIR}/customizations.sql"
   "${SQL_DIR}/nkg_post_customizations.sql"
   "${SQL_DIR}/commit.sql"

data/sql_filelist.cmake

@@ -0,0 +1,113 @@
+#!/usr/bin/env python
+###############################################################################
+# $Id$
+#
+#  Project:  PROJ
+#  Purpose:  Populate data/sql/wgs84_realizations_concatenated_operations.sql
+#  Author:   Even Rouault <even.rouault at spatialys.com>
+#
+###############################################################################
+#  Copyright (c) 2025, Even Rouault <even.rouault at spatialys.com>
+#
+# SPDX-License-Identifier: MIT
+###############################################################################
+
+import os
+
+script_dir_name = os.path.dirname(os.path.realpath(__file__))
+sql_dir_name = os.path.join(os.path.dirname(script_dir_name), 'data', 'sql')
+
+out_filename = os.path.join(sql_dir_name, 'wgs84_realizations_concatenated_operations') + '.sql'
+#print(out_filename)
+
+def sanitize_crs_name_for_code(name):
+    return name.replace('(', '_').replace(')', '').replace(' ','_').replace('__','_').upper()
+
+def gen_transformations(sql, transformations, crs_dict):
+    for i in range(len(transformations)):
+        if transformations[i][0] == "WGS 84 (G1150)" and transformations[i][1] == "WGS 84 (G1762)":
+            continue
+        for j in range(i+1,len(transformations)):
+            if transformations[j][0] == "WGS 84 (G1150)" and transformations[j][1] == "WGS 84 (G1762)":
+                continue
+            source_crs = transformations[i][0]
+            target_crs = transformations[j][1]
+
+            source_crs_week = 0 if source_crs == "WGS 84 (Transit)" else int(source_crs[len("WGS 84 (G"):-1])
+            target_crs_week = int(target_crs[len("WGS 84 (G"):-1])
+
+            transfm_code = sanitize_crs_name_for_code(source_crs) + "_TO_" + sanitize_crs_name_for_code(target_crs)
+            transfm_name = f"{source_crs} to {target_crs}"
+
+            source_crs_code = crs_dict[source_crs]
+            target_crs_code = crs_dict[target_crs]
+
+            step = 0
+            steps_sql = []
+            total_acc = 0
+            for k in range(i,j+1):
+                source_crs_name = transformations[k][0]
+                target_crs_name = transformations[k][1]
+                if source_crs_week <= 1150 and target_crs_week >= 1762 and \
+                   ((source_crs_name == "WGS 84 (G1150)" and target_crs_name == "WGS 84 (G1674)") or \
+                    (source_crs_name == "WGS 84 (G1674)" and target_crs_name == "WGS 84 (G1762)")):
+                    continue
+                step += 1
+                step_code = transformations[k][2]
+                acc = transformations[k][3]
+                total_acc += acc
+                steps_sql.append(f"INSERT INTO concatenated_operation_step VALUES('PROJ','{transfm_code}',{step},'EPSG','{step_code}','forward'); -- {source_crs_name} to {target_crs_name} (EPSG:{step_code}), {acc} m\n")
+
+            if len(steps_sql) <= 1:
+                continue
+
+            total_acc = round(total_acc * 100) / 100.0
+
+            sql += f"INSERT INTO concatenated_operation VALUES('PROJ','{transfm_code}','{transfm_name}','Transformation based on concatenation of transformations between WGS 84 realizations','EPSG','{source_crs_code}','EPSG','{target_crs_code}',{total_acc},NULL,0);\n"
+
+            for step_sql in steps_sql:
+                sql += step_sql
+
+            sql += f"INSERT INTO usage VALUES('PROJ','{transfm_code}_USAGE','concatenated_operation','PROJ','{transfm_code}',\n"
+            sql += "    'EPSG','1262', -- extent: World\n"
+            sql += "    'EPSG','1027' -- scope: Geodesy\n"
+            sql += ");\n"
+
+            sql += "\n"
+
+    return sql
+
+crs_dict = {
+    "WGS 84 (Transit)": 7815,
+    "WGS 84 (G730)": 7656,
+    "WGS 84 (G873)": 7658,
+    "WGS 84 (G1150)": 7660,
+    "WGS 84 (G1674)": 7662,
+    "WGS 84 (G1762)": 7664,
+    "WGS 84 (G2139)": 9753,
+    "WGS 84 (G2296)": 10604,
+}
+
+f = open(out_filename, 'wb')
+f.write("--- This file has been generated by scripts/build_wgs84_realizations_concatenated_operations.py. DO NOT EDIT !\n\n".encode('UTF-8'))
+
+f.write("-- Concatenated accuracy is sum of accuracies\n\n".encode('UTF-8'))
+
+sql = ""
+
+transformations = [
+    ("WGS 84 (Transit)", "WGS 84 (G730)", 9960, 0.7),  # regular Helmert
+    ("WGS 84 (G730)", "WGS 84 (G873)", 9961, 0.04),    # regular Helmert
+    ("WGS 84 (G873)", "WGS 84 (G1150)", 9962, 0.03),   # time-dependent Helmert
+    ("WGS 84 (G1150)", "WGS 84 (G1674)", 9963, 0.02),  # time-dependent Helmert
+    ("WGS 84 (G1150)", "WGS 84 (G1762)", 7668, 0.02),  # regular Helmert. Note that it doesn't go through G1674
+    ("WGS 84 (G1674)", "WGS 84 (G1762)", 7667, 0.01),  # regular Helmert
+    ("WGS 84 (G1762)", "WGS 84 (G2139)", 9756, 0.01),  # regular Helmert
+    ("WGS 84 (G2139)", "WGS 84 (G2296)", 10607, 0.01), # regular Helmert
+]
+
+sql = gen_transformations(sql, transformations, crs_dict)
+
+f.write(sql.encode('UTF-8'))
+
+f.close()

scripts/build_wgs84_realizations_concatenated_operations.py

@@ -7308,7 +7308,8 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
                                           const std::string &code) {
         return !(d->run("SELECT 1 FROM coordinate_operation_view WHERE "
                         "(source_crs_auth_name = ? AND source_crs_code = ?) OR "
-                        "(target_crs_auth_name = ? AND target_crs_code = ?)",
+                        "(target_crs_auth_name = ? AND target_crs_code = ?) "
+                        "LIMIT 1",
                         {auth_name, code, auth_name, code})
                      .empty());
     };
@@ -7407,36 +7408,50 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
         "AND v2.target_crs_auth_name = ? AND v2.target_crs_code = ? ");
     std::string minDate;
     std::string criterionOnIntermediateCRS;
+
+    const auto sourceCRS = d->createFactory(sourceCRSAuthName)
+                               ->createCoordinateReferenceSystem(sourceCRSCode);
+    const auto targetCRS = d->createFactory(targetCRSAuthName)
+                               ->createCoordinateReferenceSystem(targetCRSCode);
+
+    const bool ETRFtoETRF = starts_with(sourceCRS->nameStr(), "ETRF") &&
+                            starts_with(targetCRS->nameStr(), "ETRF");
+
     if (allowedIntermediateObjectType == ObjectType::GEOGRAPHIC_CRS) {
-        auto sourceCRS = d->createFactory(sourceCRSAuthName)
-                             ->createGeodeticCRS(sourceCRSCode);
-        auto targetCRS = d->createFactory(targetCRSAuthName)
-                             ->createGeodeticCRS(targetCRSCode);
-        const auto &sourceDatum = sourceCRS->datum();
-        const auto &targetDatum = targetCRS->datum();
-        if (sourceDatum && sourceDatum->publicationDate().has_value() &&
-            targetDatum && targetDatum->publicationDate().has_value()) {
-            const auto sourceDate(sourceDatum->publicationDate()->toString());
-            const auto targetDate(targetDatum->publicationDate()->toString());
-            minDate = std::min(sourceDate, targetDate);
-            // Check that the datum of the intermediateCRS has a publication
-            // date most recent that the one of the source and the target CRS
-            // Except when using the usual WGS84 pivot which happens to have a
-            // NULL publication date.
-            criterionOnIntermediateCRS =
-                "AND EXISTS(SELECT 1 FROM geodetic_crs x "
-                "JOIN geodetic_datum y "
-                "ON "
-                "y.auth_name = x.datum_auth_name AND "
-                "y.code = x.datum_code "
-                "WHERE "
-                "x.auth_name = v1.target_crs_auth_name AND "
-                "x.code = v1.target_crs_code AND "
-                "x.type IN ('geographic 2D', 'geographic 3D') AND "
-                "(y.publication_date IS NULL OR "
-                "(y.publication_date >= '" +
-                minDate + "'))) ";
-        } else {
+        const auto &sourceGeogCRS =
+            dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
+        const auto &targetGeogCRS =
+            dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
+        if (sourceGeogCRS && targetGeogCRS) {
+            const auto &sourceDatum = sourceGeogCRS->datum();
+            const auto &targetDatum = targetGeogCRS->datum();
+            if (sourceDatum && sourceDatum->publicationDate().has_value() &&
+                targetDatum && targetDatum->publicationDate().has_value()) {
+                const auto sourceDate(
+                    sourceDatum->publicationDate()->toString());
+                const auto targetDate(
+                    targetDatum->publicationDate()->toString());
+                minDate = std::min(sourceDate, targetDate);
+                // Check that the datum of the intermediateCRS has a publication
+                // date most recent that the one of the source and the target
+                // CRS Except when using the usual WGS84 pivot which happens to
+                // have a NULL publication date.
+                criterionOnIntermediateCRS =
+                    "AND EXISTS(SELECT 1 FROM geodetic_crs x "
+                    "JOIN geodetic_datum y "
+                    "ON "
+                    "y.auth_name = x.datum_auth_name AND "
+                    "y.code = x.datum_code "
+                    "WHERE "
+                    "x.auth_name = v1.target_crs_auth_name AND "
+                    "x.code = v1.target_crs_code AND "
+                    "x.type IN ('geographic 2D', 'geographic 3D') AND "
+                    "(y.publication_date IS NULL OR "
+                    "(y.publication_date >= '" +
+                    minDate + "'))) ";
+            }
+        }
+        if (criterionOnIntermediateCRS.empty()) {
             criterionOnIntermediateCRS =
                 "AND EXISTS(SELECT 1 FROM geodetic_crs x WHERE "
                 "x.auth_name = v1.target_crs_auth_name AND "
@@ -7586,6 +7601,33 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
     if (discardSuperseded) {
         res = filterOutSuperseded(std::move(res));
     }
+
+    const auto checkPivotForETRFToETRF =
+        [ETRFtoETRF, &sourceCRS,
+         &targetCRS](const crs::CRSPtr &intermediateCRS) {
+            // Make sure that ETRF2000 to ETRF2014 doesn't go throught ITRF9x or
+            // ITRF>2014
+            if (ETRFtoETRF && intermediateCRS &&
+                starts_with(intermediateCRS->nameStr(), "ITRF")) {
+                const auto normalizeDate = [](int v) {
+                    return (v >= 80 && v <= 99) ? v + 1900 : v;
+                };
+                const int srcDate = normalizeDate(
+                    atoi(sourceCRS->nameStr().c_str() + strlen("ETRF")));
+                const int tgtDate = normalizeDate(
+                    atoi(targetCRS->nameStr().c_str() + strlen("ETRF")));
+                const int intermDate = normalizeDate(
+                    atoi(intermediateCRS->nameStr().c_str() + strlen("ITRF")));
+                if (srcDate > 0 && tgtDate > 0 && intermDate > 0) {
+                    if (intermDate < std::min(srcDate, tgtDate) ||
+                        intermDate > std::max(srcDate, tgtDate)) {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        };
+
     for (const auto &row : res) {
         const auto &table1 = row[0];
         const auto &auth_name1 = row[1];
@@ -7604,6 +7646,9 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
                                    targetCRSAuthName, targetCRSCode)) {
                 continue;
             }
+            if (!checkPivotForETRFToETRF(op1->targetCRS())) {
+                continue;
+            }
             auto op2 =
                 d->createFactory(auth_name2)
                     ->createCoordinateOperation(
@@ -7642,6 +7687,7 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
     if (discardSuperseded) {
         res = filterOutSuperseded(std::move(res));
     }
+
     for (const auto &row : res) {
         const auto &table1 = row[0];
         const auto &auth_name1 = row[1];
@@ -7660,6 +7706,9 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
                                    targetCRSAuthName, targetCRSCode)) {
                 continue;
             }
+            if (!checkPivotForETRFToETRF(op1->targetCRS())) {
+                continue;
+            }
             auto op2 =
                 d->createFactory(auth_name2)
                     ->createCoordinateOperation(
@@ -7737,6 +7786,9 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
                                    targetCRSAuthName, targetCRSCode)) {
                 continue;
             }
+            if (!checkPivotForETRFToETRF(op1->sourceCRS())) {
+                continue;
+            }
             auto op2 =
                 d->createFactory(auth_name2)
                     ->createCoordinateOperation(
@@ -7793,6 +7845,9 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
                                    targetCRSAuthName, targetCRSCode)) {
                 continue;
             }
+            if (!checkPivotForETRFToETRF(op1->sourceCRS())) {
+                continue;
+            }
             auto op2 =
                 d->createFactory(auth_name2)
                     ->createCoordinateOperation(

src/iso19111/factory.cpp

@@ -1867,3 +1867,26 @@ tests:
       +step +inv +proj=cart +ellps=WGS84
       +step +proj=unitconvert +xy_in=rad +xy_out=deg
       +step +proj=axisswap +order=2,1
+- comment: >
+   Test that ETRF2000 to ETRF2014 only goes through ITRF2000, ITRF2008 or ITRF2014, and not ITRF9x or ITRF>2020
+  args: -s ETRF2000 -t ETRF2014 --3d --summary --authority EPSG
+  out: |
+    Candidate operations found: 6
+    unknown id, Conversion from ETRF2000 (geog3D) to ETRF2000 (geocentric) + Inverse of ITRF2014 to ETRF2000 (1) + ITRF2014 to ETRF2014 (2) + Conversion from ETRF2014 (geocentric) to ETRF2014 (geog3D), 0 m, Europe - onshore and offshore - ETRF extent - approximately 16°W to 33°E and 33°N to 84°N., time-dependent operation
+    unknown id, Conversion from ETRF2000 (geog3D) to ETRF2000 (geocentric) + Inverse of ITRF2014 to ETRF2000 (1) + ITRF2014 to ETRF2014 (1) + Conversion from ETRF2014 (geocentric) to ETRF2014 (geog3D), 0 m, Europe - onshore and offshore - ETRF extent - approximately 16°W to 33°E and 33°N to 84°N., time-dependent operation
+    unknown id, Conversion from ETRF2000 (geog3D) to ETRF2000 (geocentric) + Inverse of ITRF2008 to ETRF2000 (1) + ITRF2008 to ETRF2014 (1) + Conversion from ETRF2014 (geocentric) to ETRF2014 (geog3D), 0 m, Europe - onshore and offshore - ETRF extent - approximately 16°W to 33°E and 33°N to 84°N., time-dependent operation
+    unknown id, Conversion from ETRF2000 (geog3D) to ETRF2000 (geocentric) + Inverse of ITRF2005 to ETRF2000 (1) + ITRF2005 to ETRF2014 (1) + Conversion from ETRF2014 (geocentric) to ETRF2014 (geog3D), 0 m, Europe - onshore and offshore - ETRF extent - approximately 16°W to 33°E and 33°N to 84°N., time-dependent operation
+    unknown id, Conversion from ETRF2000 (geog3D) to ETRF2000 (geocentric) + Inverse of ITRF2000 to ETRF2000 (2) + ITRF2000 to ETRF2014 (1) + Conversion from ETRF2014 (geocentric) to ETRF2014 (geog3D), 0 m, Europe - onshore and offshore - ETRF extent - approximately 16°W to 33°E and 33°N to 84°N., time-dependent operation
+    unknown id, Conversion from ETRF2000 (geog3D) to ETRF2000 (geocentric) + Inverse of ITRF2000 to ETRF2000 (1) + ITRF2000 to ETRF2014 (1) + Conversion from ETRF2014 (geocentric) to ETRF2014 (geog3D), 0 m, Europe - onshore and offshore - ETRF extent - approximately 16°W to 33°E and 33°N to 84°N., time-dependent operation
+- comment: >
+   Test that ETRF2000-PL to ETRF2014 is not negatively affected by https://github.com/OSGeo/PROJ/pull/4364 changes (the current actual result is not golden in any way, using ETRF2000 to ETRF2014 could probably make sense)
+  args: -s ETRF2000-PL -t ETRF2014 --summary
+  out: |
+    Candidate operations found: 1
+    unknown id, ETRF2000-PL to ETRS89 (1) + ETRS89 to ETRF2014, 0.1 m, Poland - onshore and offshore.
+- comment: >
+   Test WGS 84 (G1150) to WGS 84 (G2296)
+  args: -s "WGS 84 (G1150)" -t "WGS 84 (G2296)" --summary
+  out: |
+    Candidate operations found: 1
+    unknown id, Conversion from WGS 84 (G1150) (geog2D) to WGS 84 (G1150) (geocentric) + WGS 84 (G1150) to WGS 84 (G1762) (1) + WGS 84 (G1762) to WGS 84 (G2139) (1) + WGS 84 (G2139) to WGS 84 (G2296) (1) + Conversion from WGS 84 (G2296) (geocentric) to WGS 84 (G2296) (geog2D), 0.04 m, World