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

1 files changed, 15 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;