proj_create_crs_to_crs() + proj_trans(): fix when non-Greenwich prime meridian is involved

This fixes a regression introduced in 7af1d5741da08d9546b907e0da2c21c54c61b27 / PROJ 7.2.0
where reprojection of area of use was broken when the source/target CRS
did not use Greenwich as prime meridian.
Fixes https://lists.osgeo.org/pipermail/gdal-dev/2021-October/054764.html

Now with the fix:
- using grid:
$ echo 286415 431434 | PROJ_NETWORK=ON src/cs2cs -d 4 EPSG:20790  EPSG:3763
86412.4262	131434.1706 0.0000

- not using it:
$ echo 286415 431434 | src/cs2cs -d 4 EPSG:20790  EPSG:3763
86412.5265	131433.8561 0.0000

2 files changed, 58 insertions, 6 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 6f8c3027..4e575f14 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1612,15 +1612,24 @@ static PJ* create_operation_to_geog_crs(PJ_CONTEXT* ctx, const PJ* crs) {
         geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
         geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_3D_CRS )
     {
-        auto datum = proj_crs_get_datum(ctx, geodetic_crs);
-        auto datum_ensemble = proj_crs_get_datum_ensemble(ctx, geodetic_crs);
+        auto datum = proj_crs_get_datum_forced(ctx, geodetic_crs);
+        assert( datum );
         auto cs = proj_create_ellipsoidal_2D_cs(
             ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, nullptr, 0);
-        auto temp = proj_create_geographic_crs_from_datum(
-            ctx, "unnamed crs", datum ? datum : datum_ensemble,
-            cs);
+        auto ellps = proj_get_ellipsoid(ctx, datum);
         proj_destroy(datum);
-        proj_destroy(datum_ensemble);
+        double semi_major_metre = 0;
+        double inv_flattening = 0;
+        proj_ellipsoid_get_parameters(ctx, ellps, &semi_major_metre,
+                                      nullptr, nullptr, &inv_flattening);
+        // It is critical to set the prime meridian to 0
+        auto temp = proj_create_geographic_crs(
+            ctx, "unnamed crs", "unnamed datum",
+            proj_get_name(ellps),
+            semi_major_metre, inv_flattening,
+            "Reference prime meridian", 0, nullptr, 0,
+            cs);
+        proj_destroy(ellps);
         proj_destroy(cs);
         proj_destroy(geodetic_crs);
         geodetic_crs = temp;
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index a9315fa7..fb8e79a0 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -1646,6 +1646,49 @@ TEST_F(CApi, proj_create_operations) {
 
 // ---------------------------------------------------------------------------
 
+TEST_F(CApi, proj_create_operations_prime_meridian_non_greenwich) {
+    auto ctxt = proj_create_operation_factory_context(m_ctxt, nullptr);
+    ASSERT_NE(ctxt, nullptr);
+    ContextKeeper keeper_ctxt(ctxt);
+
+    auto source_crs = proj_create_from_database(
+        m_ctxt, "EPSG", "27562", PJ_CATEGORY_CRS, false,
+        nullptr); // "NTF (Paris) / Lambert Centre France"
+    ASSERT_NE(source_crs, nullptr);
+    ObjectKeeper keeper_source_crs(source_crs);
+
+    auto target_crs = proj_create_from_database(
+        m_ctxt, "EPSG", "4258", PJ_CATEGORY_CRS, false, nullptr); // ETRS89
+    ASSERT_NE(target_crs, nullptr);
+    ObjectKeeper keeper_target_crs(target_crs);
+
+    proj_operation_factory_context_set_spatial_criterion(
+        m_ctxt, ctxt, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
+
+    proj_operation_factory_context_set_grid_availability_use(
+        m_ctxt, ctxt, PROJ_GRID_AVAILABILITY_IGNORED);
+
+    auto res = proj_create_operations(m_ctxt, source_crs, target_crs, ctxt);
+    ASSERT_NE(res, nullptr);
+    ObjListKeeper keeper_res(res);
+
+    {
+        PJ_COORD coord;
+        // lat,lon=49,-4 if using grid
+        coord.xy.x = 136555.58288992;
+        coord.xy.y = 463344.51894296;
+        int idx = proj_get_suggested_operation(m_ctxt, res, PJ_FWD, coord);
+        ASSERT_GE(idx, 0);
+        auto op = proj_list_get(m_ctxt, res, idx);
+        ASSERT_NE(op, nullptr);
+        ObjectKeeper keeper_op(op);
+        // Transformation using grid
+        EXPECT_EQ(proj_coordoperation_get_grid_used_count(m_ctxt, op), 1);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 TEST_F(CApi, proj_get_suggested_operation_with_operations_without_area_of_use) {
     auto ctxt = proj_create_operation_factory_context(m_ctxt, nullptr);
     ASSERT_NE(ctxt, nullptr);