Merge pull request #1641 from rouault/fix_proj_create_crs_to_crs_if_only_ballpark_transf

proj_create_crs_to_crs(): fix when there are only transformations with ballpark steps

2 files changed, 44 insertions, 3 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index db1d9516..4b47c4e9 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1137,6 +1137,8 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
 
     try
     {
+        bool skipDefaultTransforms = true;
+
         // Iterate over source->target candidate transformations and reproject
         // their long-lat bounding box into the source CRS.
         for( int i = 0; i < op_count; i++ )
@@ -1149,12 +1151,25 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
             double north_lat = 0.0;
 
             const char* name = proj_get_name(op);
-            if( name && (strstr(name, "Ballpark geographic offset") ||
+            bool canUseOp = true;
+            if( skipDefaultTransforms &&
+                name && (strstr(name, "Ballpark geographic offset") ||
                          strstr(name, "Ballpark geocentric translation")) )
             {
-                // Skip default transformations
+                // Skip default transformations unless there is already one at
+                // the beginning (in which case all of them will have one)
+                if( i == 0 )
+                {
+                    skipDefaultTransforms = false;
+                }
+                else
+                {
+                    canUseOp = false;
+                }
             }
-            else if( proj_get_area_of_use(ctx, op,
+
+            if( canUseOp &&
+                proj_get_area_of_use(ctx, op,
                         &west_lon, &south_lat, &east_lon, &north_lat, nullptr) )
             {
                 if( west_lon <= east_lon )
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index bd59358e..483ecdbd 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -3997,4 +3997,30 @@ TEST_F(CApi, proj_crs_create_bound_vertical_crs_to_WGS84) {
               "+vunits=m +no_defs +type=crs");
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_create_crs_to_crs_with_only_ballpark_transformations) {
+    // ETRS89 / UTM zone 31N + EGM96 height to WGS 84 (G1762)
+    auto P =
+        proj_create_crs_to_crs(m_ctxt, "EPSG:25831+5773", "EPSG:7665", nullptr);
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto Pnormalized = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper keeper_Pnormalized(Pnormalized);
+    ASSERT_NE(Pnormalized, nullptr);
+
+    PJ_COORD coord;
+    coord.xyzt.x = 500000;
+    coord.xyzt.y = 4500000;
+    coord.xyzt.z = 0;
+    coord.xyzt.t = 0;
+    coord = proj_trans(Pnormalized, PJ_FWD, coord);
+    EXPECT_NEAR(coord.xyzt.x, 3.0, 1e-9);
+    EXPECT_NEAR(coord.xyzt.y, 40.65085651660555, 1e-9);
+    if (coord.xyzt.z != 0) {
+        // z will depend if the egm96_15.gtx grid is there or not
+        EXPECT_NEAR(coord.xyzt.z, 47.04784081844435, 1e-3);
+    }
+}
+
 } // namespace