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

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 "648183209a519b679cca32c4cec988e2")
+set(PROJ_DB_SQL_EXPECTED_MD5 "8145b5734f4007f81b188ac3b690e8b3")
 
 add_custom_command(
   OUTPUT ${PROJ_DB}

data/sql/customizations.sql

@@ -581,3 +581,51 @@ FROM grid_transformation gt
 JOIN usage u ON u.object_auth_name = gt.auth_name AND u.object_code = gt.code AND u.object_table_name = 'grid_transformation'
 WHERE method_auth_name = 'EPSG' AND method_name LIKE 'Geog3D to Geog2D+%'
 AND EXISTS (SELECT 1 FROM grid_transformation gt2 WHERE gt2.auth_name = 'PROJ' AND gt2.code = gt.auth_name || '_' || gt.code || '_RESTRICTED_TO_VERTCRS');
+
+-- Add records corresponding to EGM2008 grid for WGS 84 realizations
+
+INSERT INTO "grid_transformation"
+SELECT
+    'PROJ' AS auth_name,
+    replace(replace(replace(gcrs.name, ' ', '_'), '(', ''), ')', '') || '_TO_EGM2008',
+    gcrs.name || ' to EGM2008 height (from ' || gt.name || ')' AS name,
+    gt.description,
+    gt.method_auth_name,
+    gt.method_code,
+    gt.method_name,
+    gcrs.auth_name,
+    gcrs.code,
+    gt.target_crs_auth_name,
+    gt.target_crs_code,
+    gt.accuracy,
+    gt.grid_param_auth_name,
+    gt.grid_param_code,
+    gt.grid_param_name,
+    gt.grid_name,
+    gt.grid2_param_auth_name,
+    gt.grid2_param_code,
+    gt.grid2_param_name,
+    gt.grid2_name,
+    gt.interpolation_crs_auth_name,
+    gt.interpolation_crs_code,
+    gt.operation_version,
+    gt.deprecated
+FROM grid_transformation gt, geodetic_crs gcrs
+WHERE gt.name = 'WGS 84 to EGM2008 height (1)'
+AND gcrs.auth_name = 'EPSG' AND gcrs.name LIKE 'WGS 84 (G%' AND gcrs.type='geographic 3D' and gcrs.deprecated=0;
+
+INSERT INTO "usage"
+SELECT
+    'PROJ' AS auth_name,
+    'USAGE_' || replace(replace(replace(gcrs.name, ' ', '_'), '(', ''), ')', '') || '_TO_EGM2008' AS code,
+    'grid_transformation' AS object_table_name,
+    'PROJ' AS object_auth_name,
+    replace(replace(replace(gcrs.name, ' ', '_'), '(', ''), ')', '') || '_TO_EGM2008' AS object_code,
+    u.extent_auth_name,
+    u.extent_code,
+    u.scope_auth_name,
+    u.scope_code
+FROM grid_transformation gt, geodetic_crs gcrs
+JOIN usage u ON u.object_auth_name = gt.auth_name AND u.object_code = gt.code AND u.object_table_name = 'grid_transformation'
+WHERE gt.name = 'WGS 84 to EGM2008 height (1)'
+AND gcrs.auth_name = 'EPSG' AND gcrs.name LIKE 'WGS 84 (G%' AND gcrs.type='geographic 3D' and gcrs.deprecated=0;

data/sql_filelist.cmake

@@ -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"

scripts/build_wgs84_realizations_concatenated_operations.py

@@ -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()

src/iso19111/factory.cpp

@@ -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/operation/singleoperation.cpp

@@ -3482,8 +3482,28 @@ bool SingleOperation::exportToPROJStringGeneric(
 
         auto sourceCRSGeod =
             dynamic_cast<const crs::GeodeticCRS *>(sourceCRS().get());
+        if (!sourceCRSGeod) {
+            auto sourceCRSCompound =
+                dynamic_cast<const crs::CompoundCRS *>(sourceCRS().get());
+            if (sourceCRSCompound) {
+                sourceCRSGeod = dynamic_cast<const crs::GeodeticCRS *>(
+                    sourceCRSCompound->componentReferenceSystems()
+                        .front()
+                        .get());
+            }
+        }
         auto targetCRSGeod =
             dynamic_cast<const crs::GeodeticCRS *>(targetCRS().get());
+        if (!targetCRSGeod) {
+            auto targetCRSCompound =
+                dynamic_cast<const crs::CompoundCRS *>(targetCRS().get());
+            if (targetCRSCompound) {
+                targetCRSGeod = dynamic_cast<const crs::GeodeticCRS *>(
+                    targetCRSCompound->componentReferenceSystems()
+                        .front()
+                        .get());
+            }
+        }
         if (sourceCRSGeod && targetCRSGeod) {
             auto sourceCRSGeog =
                 dynamic_cast<const crs::GeographicCRS *>(sourceCRSGeod);