diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql
index 858a4b774c..d47212bc44 100644
--- a/data/sql/customizations.sql
+++ b/data/sql/customizations.sql
@@ -563,7 +563,7 @@ JOIN compound_crs c ON gt.target_crs_code = c.code AND gt.target_crs_auth_name =
 JOIN geodetic_crs gcrs ON gt.source_crs_auth_name = gcrs.auth_name AND gt.source_crs_code = gcrs.code
 JOIN vertical_crs vcrs on vcrs.auth_name = c.vertical_crs_auth_name AND vcrs.code = c.vertical_crs_code
 WHERE method_auth_name = 'EPSG' AND method_name LIKE 'Geog3D to Geog2D+%'
-AND NOT EXISTS (SELECT 1 FROM grid_transformation gt2 WHERE gt2.method_name LIKE 'Geographic3D to%' AND gt2.source_crs_auth_name = gt.source_crs_auth_name AND gt2.source_crs_code = gt.source_crs_code AND gt2.target_crs_auth_name = vcrs.auth_name AND gt2.target_crs_code = vcrs.code)
+AND NOT EXISTS (SELECT 1 FROM grid_transformation gt2 WHERE gt2.method_name LIKE 'Geographic3D to%' AND gt2.source_crs_auth_name = gt.source_crs_auth_name AND gt2.source_crs_code = gt.source_crs_code AND gt2.target_crs_auth_name = vcrs.auth_name AND gt2.target_crs_code = vcrs.code AND gt2.grid_name = gt.grid_name)
 AND gt.deprecated = 0;
 
 INSERT INTO "usage"
diff --git a/data/sql/grid_alternatives.sql b/data/sql/grid_alternatives.sql
index 7ad3af99ab..8decfefeb9 100644
--- a/data/sql/grid_alternatives.sql
+++ b/data/sql/grid_alternatives.sql
@@ -102,6 +102,9 @@ VALUES
 ('chgeo2004_ETRS.agr','ch_swisstopo_chgeo2004_ETRS89_LHN95.tif',NULL,'GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ch_swisstopo_chgeo2004_ETRS89_LHN95.tif',1,1,NULL),
 ('chgeo2004_htrans_ETRS.agr','ch_swisstopo_chgeo2004_ETRS89_LN02.tif',NULL,'GTiff','geoid_like',0,NULL,'https://cdn.proj.org/ch_swisstopo_chgeo2004_ETRS89_LN02.tif',1,1,NULL),
 
+-- cz_cuzk - ČÚZK
+('cz_cuzk_CR-2005.tif','cz_cuzk_CR-2005.tif',NULL,'GTiff','geoid_like',0,NULL,'https://cdn.proj.org/cz_cuzk_CR-2005.tif',1,1,NULL),
+
 -- de_adv - Arbeitsgemeinschaft der Vermessungsverwaltungender der Länder der Bundesrepublik Deutschland (AdV)
 ('BETA2007.gsb','de_adv_BETA2007.tif','BETA2007.gsb','GTiff','hgridshift',0,NULL,'https://cdn.proj.org/de_adv_BETA2007.tif',1,1,NULL),
 
diff --git a/data/sql/transformations_czechia.sql b/data/sql/transformations_czechia.sql
new file mode 100644
index 0000000000..94c9f1ee39
--- /dev/null
+++ b/data/sql/transformations_czechia.sql
@@ -0,0 +1,24 @@
+-- This file is hand generated.
+
+-- Czechia transformations
+
+-- Geoid transformation
+
+INSERT INTO "grid_transformation" VALUES(
+    'PROJ','ETRS89_TO_BALTIC_HEIGHT_CZECHIA','ETRS89 to Baltic 1957 height (Czechia)',
+    NULL,
+    'EPSG','9665','Geographic3D to GravityRelatedHeight (gtx)',
+    'EPSG','4937', -- source CRS (ETRS 89)
+    'EPSG','8357', -- target CRS (Baltic 1957 height)
+    0.05, -- guessed...
+    'EPSG','8666','Geoid (height correction) model file','cz_cuzk_CR-2005.tif',
+    NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
+
+INSERT INTO "usage" VALUES(
+    'PROJ',
+    'ETRS89_TO_BALTIC_HEIGHT_CZECHIA_USAGE',
+    'grid_transformation',
+    'PROJ',
+    'ETRS89_TO_BALTIC_HEIGHT_CZECHIA',
+    'EPSG','1079','EPSG','1189'
+);
diff --git a/data/sql/transformations_czechia_extra.sql b/data/sql/transformations_czechia_extra.sql
new file mode 100644
index 0000000000..a14b8de452
--- /dev/null
+++ b/data/sql/transformations_czechia_extra.sql
@@ -0,0 +1,170 @@
+-- This file is hand generated.
+
+-- WARNING: for now this file is not integrated in the default build, due to
+-- licensing of grids cz_cuzk_table_yx_3_v1710_east_north.tif and cz_cuzk_CR-2005.tif
+-- being not clarified.
+
+-- Czechia transformations
+
+-- 2d transformations between S-JTSK / Krovak and S-JTSK/05 / Modified Krovak using a grid with easting,northing offsets
+
+INSERT INTO other_transformation VALUES(
+    'PROJ','S_JTSK_E_N_TO_S_JTSK05_E_N','S-JTSK / Krovak East North (EPSG:5514) to S-JTSK/05 / Modified Krovak East North (EPSG:5516)',
+    'Transformation based on grid table_yx_3_v1710.dat',
+    'PROJ','PROJString',
+    '+proj=gridshift +grids=cz_cuzk_table_yx_3_v1710_east_north.tif',
+    'EPSG','5514','EPSG','5516',0.035,
+    NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
+INSERT INTO "usage" VALUES('PROJ','S_JTSK_E_N_TO_S_JTSK05_E_N_USAGE','other_transformation','PROJ','S_JTSK_E_N_TO_S_JTSK05_E_N','EPSG','1079','EPSG','1189');
+
+
+INSERT INTO other_transformation VALUES(
+    'PROJ','S_JTSK_E_N_TO_S_JTSK05','S-JTSK / Krovak East North (EPSG:5514) to S-JTSK/05 / Modified Krovak (EPSG:5515)',
+    'Transformation based on grid table_yx_3_v1710.dat',
+    'PROJ','PROJString',
+    '+proj=pipeline ' ||
+        '+step +proj=gridshift +grids=cz_cuzk_table_yx_3_v1710_east_north.tif ' ||
+        '+step +proj=axisswap +order=-2,-1',
+    'EPSG','5514','EPSG','5515',0.035,
+    NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
+INSERT INTO "usage" VALUES('PROJ','S_JTSK_E_N_TO_S_JTSK05_USAGE','other_transformation','PROJ','S_JTSK_E_N_TO_S_JTSK05','EPSG','1079','EPSG','1189');
+
+
+INSERT INTO other_transformation VALUES(
+    'PROJ','S_JTSK_TO_S_JTSK05_E_N','S-JTSK / Krovak (EPSG:5513) to S-JTSK/05 / Modified Krovak East North (EPSG:5516)',
+    'Transformation based on grid table_yx_3_v1710.dat',
+    'PROJ','PROJString',
+    '+proj=pipeline ' ||
+        '+step +proj=axisswap +order=-2,-1 ' ||
+        '+step +proj=gridshift +grids=cz_cuzk_table_yx_3_v1710_east_north.tif',
+    'EPSG','5513','EPSG','5516',0.035,
+    NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
+INSERT INTO "usage" VALUES('PROJ','S_JTSK_TO_S_JTSK05_E_N_USAGE','other_transformation','PROJ','S_JTSK_TO_S_JTSK05_E_N','EPSG','1079','EPSG','1189');
+
+
+INSERT INTO other_transformation VALUES(
+    'PROJ','S_JTSK_TO_S_JTSK05','S-JTSK / Krovak (EPSG:5513) to S-JTSK/05 / Modified Krovak (EPSG:5515)',
+    'Transformation based on grid table_yx_3_v1710.dat',
+    'PROJ','PROJString',
+    '+proj=pipeline ' ||
+        '+step +proj=axisswap +order=-2,-1 ' ||
+        '+step +proj=gridshift +grids=cz_cuzk_table_yx_3_v1710_east_north.tif ' ||
+        '+step +proj=axisswap +order=-2,-1',
+    'EPSG','5513','EPSG','5515',0.035,
+    NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
+INSERT INTO "usage" VALUES('PROJ','S_JTSK_TO_S_JTSK05_USAGE','other_transformation','PROJ','S_JTSK_TO_S_JTSK05','EPSG','1079','EPSG','1189');
+
+--- Geographic transformation: S/JTSK <--> S/JTSK/05
+
+UPDATE other_transformation SET accuracy = 0.1 WHERE name = 'S-JTSK to S-JTSK/05 (1)';
+
+INSERT INTO "concatenated_operation" VALUES(
+    'PROJ','S_JTSK_GEOG_TO_S_JTSK05_GEOG','S-JTSK (EPSG:4156) to S-JTSK/05 (EPSG:5228)',
+    'Transformation based on grid table_yx_3_v1710.dat','EPSG','4156','EPSG','5228',NULL,NULL,0);
+INSERT INTO "concatenated_operation_step" VALUES('PROJ','S_JTSK_GEOG_TO_S_JTSK05_GEOG',1,'EPSG','5510');
+INSERT INTO "concatenated_operation_step" VALUES('PROJ','S_JTSK_GEOG_TO_S_JTSK05_GEOG',2,'PROJ','S_JTSK_E_N_TO_S_JTSK05_E_N');
+INSERT INTO "concatenated_operation_step" VALUES('PROJ','S_JTSK_GEOG_TO_S_JTSK05_GEOG',3,'EPSG','5512');
+INSERT INTO "usage" VALUES('PROJ','S_JTSK_GEOG_TO_S_JTSK05_GEOG_USAGE','concatenated_operation','PROJ','S_JTSK_GEOG_TO_S_JTSK05_GEOG','EPSG','1079','EPSG','1189');
+
+
+-- CUZK recommands to do ETRS89 to S-JTSK / Krovak by doing:
+-- 1) ETRS89 to S-JTSK/05 using a Helmert transformation (EPSG:5226)
+-- 2) S-JTSK/05 to S-JTSK/05 / Modified Krovak, using modified Krovak projection (EPSG:5512)
+-- 3) S-JTSK/05 / Modified Krovak to S-JTSK / Krovak using (reverse) grid table_yx_3_v1710.dat
+-- Cf https://www.cuzk.cz/Zememerictvi/Geodeticke-zaklady-na-uzemi-CR/GNSS/Nova-realizace-systemu-ETRS89-v-CR/Metodika-prevodu-ETRF2000-vs-S-JTSK-var2(101208).aspx (Metodika převodu mezi ETRF2000 a S-JTSK varianta 2)
+INSERT INTO "concatenated_operation" VALUES(
+    'PROJ','ETRS89_TO_S_JTSK_E_N','ETRS89 to S-JTSK / Krovak East North (EPSG:5514)',
+    'Transformation based on grid table_yx_3_v1710.dat','EPSG','4258','EPSG','5514',NULL,NULL,0);
+INSERT INTO "concatenated_operation_step" VALUES('PROJ','ETRS89_TO_S_JTSK_E_N',1,'EPSG','5226');
+INSERT INTO "concatenated_operation_step" VALUES('PROJ','ETRS89_TO_S_JTSK_E_N',2,'EPSG','5512');
+INSERT INTO "concatenated_operation_step" VALUES('PROJ','ETRS89_TO_S_JTSK_E_N',3,'PROJ','S_JTSK_E_N_TO_S_JTSK05_E_N'); -- in reverse direction
+INSERT INTO "usage" VALUES('PROJ','ETRS89_TO_S_JTSK_E_N_USAGE','concatenated_operation','PROJ','ETRS89_TO_S_JTSK_E_N','EPSG','1079','EPSG','1189');
+
+
+-- 3d transformations between ETRS 89 and compound CRS S-JTSK / Krovak (or S-JTSK/05 / Modified Krovak) + Baltic 1957 height
+
+INSERT INTO "compound_crs" VALUES('PROJ','S_JTSK05_E_N_BALTIC_HEIGHT','S-JTSK/05 / Modified Krovak East North + Baltic 1957 height',NULL,'EPSG','5516','EPSG','8357',0);
+INSERT INTO "usage" VALUES('PROJ','S_JTSK05_E_N_BALTIC_HEIGHT_USAGE','compound_crs','PROJ','S_JTSK05_E_N_BALTIC_HEIGHT','EPSG','1079','EPSG','1189');
+
+-- Simulate a variable
+CREATE TEMPORARY TABLE temp_var(k,v);
+
+INSERT INTO temp_var VALUES(
+    'ETRS89_3D_TO_S_JTSK05_E_N_BALTIC_HEIGHT',
+    '+proj=pipeline ' ||
+        '+step +proj=axisswap +order=2,1 ' ||
+        '+step +proj=unitconvert +xy_in=deg +xy_out=rad ' ||
+        '+step +proj=push +v_3 +omit_inv ' ||
+        '+step +inv +proj=vgridshift +grids=cz_cuzk_CR-2005.tif +multiplier=1 ' ||
+        '+step +proj=push +v_3 +omit_fwd ' ||  -- on reverse path, restore initial Baltic height
+        '+step +proj=push +v_4 ' ||
+        '+step +proj=set +v_4=0 +omit_inv ' ||
+        '+step +proj=axisswap +order=1,2,4,3 +omit_inv ' ||  -- on forward path, save Baltic height in v_4 component...
+        '+step +proj=pop +v_3 +omit_inv ' || -- on forward parth, restore initial ellipsoidal height
+        '+step +proj=cart +ellps=GRS80 ' ||
+        '+step +inv +proj=helmert +x=572.213 +y=85.334 +z=461.94 +rx=-4.9732 +ry=-1.529 +rz=-5.2484 +s=3.5378 +convention=coordinate_frame ' ||
+        '+step +inv +proj=cart +ellps=bessel ' ||
+        '+step +proj=mod_krovak +lat_0=49.5 +lon_0=24.8333333333333 +alpha=30.2881397222222 +k=0.9999 +x_0=5000000 +y_0=5000000 +ellps=bessel ' ||
+        '+step +proj=axisswap +order=1,2,4,3 +omit_inv ' || -- on forward path, restore Baltic height from v_4 component...
+        '+step +proj=set +v_4=0 +omit_inv ' ||
+        '+step +proj=pop +v_4 ' ||
+        '+step +proj=pop +v_3 +omit_fwd '  -- on reverse path, save initial Baltic height
+);
+
+INSERT INTO other_transformation VALUES(
+    'PROJ','ETRS89_3D_TO_S_JTSK05_E_N_BALTIC_HEIGHT','ETRS89 to S-JTSK/05 / Modified Krovak East North + Baltic 1957 height',
+    'Transformation based on grid CR-2005.gtx',
+    'PROJ','PROJString',
+    (SELECT v FROM temp_var WHERE k = 'ETRS89_3D_TO_S_JTSK05_E_N_BALTIC_HEIGHT'),
+    'EPSG','4937','PROJ','S_JTSK05_E_N_BALTIC_HEIGHT',
+    0.05, -- guessed...
+    NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
+INSERT INTO "usage" VALUES('PROJ','ETRS89_3D_TO_S_JTSK05_E_N_BALTIC_HEIGHT_USAGE','other_transformation','PROJ','ETRS89_3D_TO_S_JTSK05_E_N_BALTIC_HEIGHT','EPSG','1079','EPSG','1189');
+
+
+INSERT INTO "compound_crs" VALUES('PROJ','S_JTSK05_BALTIC_HEIGHT','S-JTSK/05 / Modified Krovak + Baltic 1957 height',NULL,'EPSG','5515','EPSG','8357',0);
+INSERT INTO "usage" VALUES('PROJ','S_JTSK05_BALTIC_HEIGHT_USAGE','compound_crs','PROJ','S_JTSK05_BALTIC_HEIGHT','EPSG','1079','EPSG','1095');
+
+INSERT INTO other_transformation VALUES(
+    'PROJ','ETRS89_3D_TO_S_JTSK05_BALTIC_HEIGHT','ETRS89 to S-JTSK/05 / Modified Krovak + Baltic 1957 height',
+    'Transformation based on grid CR-2005.gtx',
+    'PROJ','PROJString',
+    (SELECT v FROM temp_var WHERE k = 'ETRS89_3D_TO_S_JTSK05_E_N_BALTIC_HEIGHT') ||
+        '+step +proj=axisswap +order=-2,-1', -- East North --> Southing Westing
+    'EPSG','4937','PROJ','S_JTSK05_BALTIC_HEIGHT',
+    0.05, -- guessed...
+    NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
+INSERT INTO "usage" VALUES('PROJ','ETRS89_3D_TO_S_JTSK05_BALTIC_HEIGHT_USAGE','other_transformation','PROJ','ETRS89_3D_TO_S_JTSK05_BALTIC_HEIGHT','EPSG','1079','EPSG','1095');
+
+
+INSERT INTO "compound_crs" VALUES('PROJ','S_JTSK_E_N_BALTIC_HEIGHT','S-JTSK / Krovak East North + Baltic 1957 height',NULL,'EPSG','5514','EPSG','8357',0);
+INSERT INTO "usage" VALUES('PROJ','S_JTSK_E_N_BALTIC_HEIGHT_USAGE','compound_crs','PROJ','S_JTSK_E_N_BALTIC_HEIGHT','EPSG','1079','EPSG','1189');
+
+INSERT INTO other_transformation VALUES(
+    'PROJ','ETRS89_3D_TO_S_JTSK_E_N_BALTIC_HEIGHT','ETRS89 to S-JTSK / Krovak East North + Baltic 1957 height',
+    'Transformation based on grids CR-2005.gtx and table_yx_3_v1710.dat ',
+    'PROJ','PROJString',
+    (SELECT v FROM temp_var WHERE k = 'ETRS89_3D_TO_S_JTSK05_E_N_BALTIC_HEIGHT') ||
+        '+step +inv +proj=gridshift +grids=cz_cuzk_table_yx_3_v1710_east_north.tif ',
+    'EPSG','4937','PROJ','S_JTSK_E_N_BALTIC_HEIGHT',
+    0.05, -- guessed...
+    NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
+INSERT INTO "usage" VALUES('PROJ','ETRS89_3D_TO_S_JTSK_E_N_BALTIC_HEIGHT_USAGE','other_transformation','PROJ','ETRS89_3D_TO_S_JTSK_E_N_BALTIC_HEIGHT','EPSG','1079','EPSG','1189');
+
+
+INSERT INTO "compound_crs" VALUES('PROJ','S_JTSK_BALTIC_HEIGHT','S-JTSK / Krovak + Baltic 1957 height',NULL,'EPSG','5513','EPSG','8357',0);
+INSERT INTO "usage" VALUES('PROJ','S_JTSK_BALTIC_HEIGHT_USAGE','compound_crs','PROJ','S_JTSK_BALTIC_HEIGHT','EPSG','1079','EPSG','1095');
+
+INSERT INTO other_transformation VALUES(
+    'PROJ','ETRS89_3D_TO_S_JTSK_BALTIC_HEIGHT','ETRS89 to S-JTSK / Krovak + Baltic 1957 height',
+    'Transformation based on grids CR-2005.gtx and table_yx_3_v1710.dat ',
+    'PROJ','PROJString',
+    (SELECT v FROM temp_var WHERE k = 'ETRS89_3D_TO_S_JTSK05_E_N_BALTIC_HEIGHT') ||
+        '+step +inv +proj=gridshift +grids=cz_cuzk_table_yx_3_v1710_east_north.tif ' ||
+        '+step +proj=axisswap +order=-2,-1', -- East North --> Southing Westing
+    'EPSG','4937','PROJ','S_JTSK_BALTIC_HEIGHT',
+    0.05, -- guessed...
+    NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
+INSERT INTO "usage" VALUES('PROJ','ETRS89_3D_TO_S_JTSK_BALTIC_HEIGHT_USAGE','other_transformation','PROJ','ETRS89_3D_TO_S_JTSK_BALTIC_HEIGHT','EPSG','1079','EPSG','1095');
+
+DROP TABLE temp_var;
diff --git a/data/sql_filelist.cmake b/data/sql_filelist.cmake
index b88451ecdb..099b0771db 100644
--- a/data/sql_filelist.cmake
+++ b/data/sql_filelist.cmake
@@ -36,6 +36,9 @@ set(SQL_FILES
   "${SQL_DIR}/nkg.sql"
   "${SQL_DIR}/iau.sql"
   "${SQL_DIR}/nrcan.sql"
+  "${SQL_DIR}/transformations_czechia.sql"
+  # Below file not yet integrated to unclear licensing of referenced grid file
+  # "${SQL_DIR}/transformations_czechia_extra.sql"
   "${SQL_DIR}/grid_alternatives.sql"
   "${SQL_DIR}/grid_alternatives_generated_noaa.sql"
   "${SQL_DIR}/nadcon5_concatenated_operations.sql"
diff --git a/src/iso19111/operation/concatenatedoperation.cpp b/src/iso19111/operation/concatenatedoperation.cpp
index 6fe3462c0d..6345cb9ce9 100644
--- a/src/iso19111/operation/concatenatedoperation.cpp
+++ b/src/iso19111/operation/concatenatedoperation.cpp
@@ -309,6 +309,36 @@ void ConcatenatedOperation::fixStepsDirection(
         }
     }
 
+    // If the first operation is a transformation whose target CRS matches the
+    // source CRS of the concatenated operation, then reverse it.
+    if (operationsInOut.size() >= 2) {
+        auto &op = operationsInOut.front();
+        auto l_sourceCRS = op->sourceCRS();
+        auto l_targetCRS = op->targetCRS();
+        if (l_sourceCRS && l_targetCRS &&
+            !areCRSMoreOrLessEquivalent(l_sourceCRS.get(),
+                                        concatOpSourceCRS.get()) &&
+            areCRSMoreOrLessEquivalent(l_targetCRS.get(),
+                                       concatOpSourceCRS.get())) {
+            op = op->inverse();
+        }
+    }
+
+    // If the last operation is a transformation whose source CRS matches the
+    // target CRS of the concatenated operation, then reverse it.
+    if (operationsInOut.size() >= 2) {
+        auto &op = operationsInOut.back();
+        auto l_sourceCRS = op->sourceCRS();
+        auto l_targetCRS = op->targetCRS();
+        if (l_sourceCRS && l_targetCRS &&
+            !areCRSMoreOrLessEquivalent(l_targetCRS.get(),
+                                        concatOpTargetCRS.get()) &&
+            areCRSMoreOrLessEquivalent(l_sourceCRS.get(),
+                                       concatOpTargetCRS.get())) {
+            op = op->inverse();
+        }
+    }
+
     const auto extractDerivedCRS =
         [](const crs::CRS *crs) -> const crs::DerivedCRS * {
         auto derivedCRS = dynamic_cast<const crs::DerivedCRS *>(crs);
diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp
index 115d7e0ef8..8f8854257c 100644
--- a/test/unit/test_factory.cpp
+++ b/test/unit/test_factory.cpp
@@ -2362,6 +2362,69 @@ TEST(factory, AuthorityFactory_getAvailableGeoidmodels) {
 
 // ---------------------------------------------------------------------------
 
+TEST_F(FactoryWithTmpDatabase,
+       AuthorityFactory_test_inversion_first_and_last_steps_of_concat_op) {
+    createStructure();
+    populateWithFakeEPSG();
+
+    // Completely dummy, to test proper inversion of first and last
+    // steps in ConcatenatedOperation, when it is needed
+    ASSERT_TRUE(execute("INSERT INTO geodetic_datum "
+                        "VALUES('EPSG','OTHER_DATUM','Other datum','',"
+                        "'EPSG','7030','EPSG','8901',NULL,NULL,NULL,"
+                        "'my anchor',NULL,0);"))
+        << last_error();
+    ASSERT_TRUE(
+        execute("INSERT INTO geodetic_crs VALUES('EPSG','OTHER_GEOG_CRS',"
+                "'OTHER_GEOG_CRS',NULL,'geographic 2D','EPSG','6422',"
+                "'EPSG','OTHER_DATUM',NULL,0);"))
+        << last_error();
+    ASSERT_TRUE(
+        execute("INSERT INTO other_transformation "
+                "VALUES('EPSG','4326_TO_OTHER_GEOG_CRS','name',NULL,"
+                "'EPSG','9601','Longitude rotation',"
+                "'EPSG','4326','EPSG','OTHER_GEOG_CRS',0.0,'EPSG'"
+                ",'8602','Longitude "
+                "offset',-17.4,'EPSG','9110',NULL,NULL,NULL,NULL,NULL,NULL,"
+                "NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,"
+                "NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,"
+                "NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);"))
+        << last_error();
+    ASSERT_TRUE(
+        execute("INSERT INTO other_transformation "
+                "VALUES('EPSG','OTHER_GEOG_CRS_TO_4326','name',NULL,"
+                "'EPSG','9601','Longitude rotation',"
+                "'EPSG','OTHER_GEOG_CRS','EPSG','4326',0.0,'EPSG'"
+                ",'8602','Longitude "
+                "offset',17.4,'EPSG','9110',NULL,NULL,NULL,NULL,NULL,NULL,"
+                "NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,"
+                "NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,"
+                "NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);"))
+        << last_error();
+    ASSERT_TRUE(execute("INSERT INTO concatenated_operation "
+                        "VALUES('EPSG','DUMMY_CONCATENATED_2','name',NULL,"
+                        "'EPSG','4326','EPSG'"
+                        ",'4326',NULL,NULL,0);"))
+        << last_error();
+    ASSERT_TRUE(execute("INSERT INTO concatenated_operation_step "
+                        "VALUES('EPSG','DUMMY_CONCATENATED_2',1,"
+                        "'EPSG','OTHER_GEOG_CRS_TO_4326');"))
+        << last_error();
+
+    ASSERT_TRUE(execute("INSERT INTO concatenated_operation_step "
+                        "VALUES('EPSG','DUMMY_CONCATENATED_2',2,"
+                        "'EPSG','4326_TO_OTHER_GEOG_CRS');"))
+        << last_error();
+
+    auto factoryEPSG = AuthorityFactory::create(DatabaseContext::create(m_ctxt),
+                                                std::string("EPSG"));
+    EXPECT_TRUE(nn_dynamic_pointer_cast<ConcatenatedOperation>(
+                    factoryEPSG->createObject("DUMMY_CONCATENATED_2")) !=
+                nullptr);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST_F(FactoryWithTmpDatabase,
        AuthorityFactory_test_with_fake_EPSG_and_OTHER_database) {
     createStructure();
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index 0b52a87b34..f6b5f469eb 100644
--- a/test/unit/test_operationfactory.cpp
+++ b/test/unit/test_operationfactory.cpp
@@ -813,13 +813,23 @@ TEST(operation, vertCRS_to_geogCRS_context) {
             authFactory->createCoordinateReferenceSystem("8357"),
             // ETRS89
             authFactory->createCoordinateReferenceSystem("4937"), ctxt);
-        ASSERT_EQ(list.size(), 2U);
+        ASSERT_EQ(list.size(), 3U);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
             "+proj=pipeline "
             "+step +proj=axisswap +order=2,1 "
             "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
             "+step +proj=vgridshift "
+            "+grids=cz_cuzk_CR-2005.tif "
+            "+multiplier=1 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+            "+step +proj=axisswap +order=2,1");
+        EXPECT_EQ(
+            list[1]->exportToPROJString(PROJStringFormatter::create().get()),
+            "+proj=pipeline "
+            "+step +proj=axisswap +order=2,1 "
+            "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+            "+step +proj=vgridshift "
             "+grids=sk_gku_Slovakia_ETRS89h_to_Baltic1957.tif "
             "+multiplier=1 "
             "+step +proj=unitconvert +xy_in=rad +xy_out=deg "