proj_create_crs_to_crs(): avoid potential reprojection failures when reprojecting area of use to source and target CRS

Was found with https://github.com/OSGeo/PROJ/pull/1989
when using cs2cs EPSG:4937 EPSG:31258+5778
- We do not need to do vertical transformation in that context. This failed here
  because the Austrian grids have nodata value outside of the shape of Austria, so
  the edges of the grids are mostly nodata values.
- And we should avoid grid-based transformations too.

1 files changed, 24 insertions, 5 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 88a0ed14..69790dde 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1058,20 +1058,39 @@ static PJ* create_operation_to_geog_crs(PJ_CONTEXT* ctx, const PJ* crs) {
         ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
     proj_operation_factory_context_set_grid_availability_use(
         ctx, operation_ctx, PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
+    auto target_crs_2D = proj_crs_demote_to_2D(ctx, nullptr, crs);
     auto op_list_to_geodetic = proj_create_operations(
-        ctx, geodetic_crs, crs, operation_ctx);
+        ctx, geodetic_crs, target_crs_2D, operation_ctx);
+    proj_destroy(target_crs_2D);
     proj_operation_factory_context_destroy(operation_ctx);
     proj_destroy(geodetic_crs);
 
-    if( op_list_to_geodetic == nullptr ||
-        proj_list_get_count(op_list_to_geodetic) == 0 )
+    const int nOpCount = op_list_to_geodetic == nullptr  ? 0 :
+                         proj_list_get_count(op_list_to_geodetic);
+    if( nOpCount == 0 )
     {
         proj_context_log_debug(ctx, "Cannot compute transformation from geographic CRS to CRS");
         proj_list_destroy(op_list_to_geodetic);
         return nullptr;
     }
-    auto opGeogToCrs = proj_list_get(ctx, op_list_to_geodetic, 0);
-    assert(opGeogToCrs);
+    PJ* opGeogToCrs = nullptr;
+    // Use in priority operations *without* grids
+    for(int i = 0; i < nOpCount; i++ )
+    {
+        auto op = proj_list_get(ctx, op_list_to_geodetic, i);
+        assert(op);
+        if( proj_coordoperation_get_grid_used_count(ctx, op) == 0 )
+        {
+            opGeogToCrs = op;
+            break;
+        }
+        proj_destroy(op);
+    }
+    if( opGeogToCrs == nullptr )
+    {
+        opGeogToCrs = proj_list_get(ctx, op_list_to_geodetic, 0);
+        assert(opGeogToCrs);
+    }
     proj_list_destroy(op_list_to_geodetic);
     return opGeogToCrs;
 }