Merge pull request #4363 from rouault/fix_4362

Database: replicate EGM2008 grid transformation record to WGS 84 realizations

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 "648183209a519b679cca32c4cec988e2")
+set(PROJ_DB_SQL_EXPECTED_MD5 "5e23b08c318cecbd9cace26233f8f5d4")
 
 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;

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);

test/cli/test_projinfo.yaml

@@ -1838,3 +1838,32 @@ tests:
     Test that "CGVD28 height" to "CGVD28(HTv2.0) height" is a no-op
   args: -s "CGVD28 height" -t "CGVD28(HTv2.0) height" -o PROJ -q
   out: +proj=noop
+- comment: >
+   Test custom records for WGS 84 realizations to EGM2008 (https://github.com/OSGeo/PROJ/issues/4362)
+  args: -s "WGS 84 (G1150) + EGM2008 height" -t "WGS 84 (G1674)" --3d -o PROJ -q
+  out: |
+    +proj=pipeline
+      +step +proj=axisswap +order=2,1
+      +step +proj=unitconvert +xy_in=deg +xy_out=rad
+      +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
+      +step +proj=cart +ellps=WGS84
+      +step +proj=helmert +x=-0.0024 +y=0.0016 +z=0.0232 +rx=-0.00027 +ry=0.00027
+            +rz=-0.00038 +s=0.00208 +dx=-0.0001 +dy=-0.0001 +dz=0.0018 +drx=0 +dry=0
+            +drz=0 +ds=-8e-05 +t_epoch=2005 +convention=coordinate_frame
+      +step +inv +proj=cart +ellps=WGS84
+      +step +proj=unitconvert +xy_in=rad +xy_out=deg
+      +step +proj=axisswap +order=2,1
+- comment: >
+   Test 'WGS 84 (Gxxx) + EGM2008 height' to 'WGS 84 (Gyyy) + EGM2008 height'
+  args: -s "WGS 84 (G1150) + EGM2008 height" -t "WGS 84 (G1674) + EGM2008 height" -o PROJ -q
+  out: |
+    +proj=pipeline
+      +step +proj=axisswap +order=2,1
+      +step +proj=unitconvert +xy_in=deg +xy_out=rad
+      +step +proj=cart +ellps=WGS84
+      +step +proj=helmert +x=-0.0024 +y=0.0016 +z=0.0232 +rx=-0.00027 +ry=0.00027
+            +rz=-0.00038 +s=0.00208 +dx=-0.0001 +dy=-0.0001 +dz=0.0018 +drx=0 +dry=0
+            +drz=0 +ds=-8e-05 +t_epoch=2005 +convention=coordinate_frame
+      +step +inv +proj=cart +ellps=WGS84
+      +step +proj=unitconvert +xy_in=rad +xy_out=deg
+      +step +proj=axisswap +order=2,1