proj_factor: fix when input is a compound CRS of a projected CRS, and…

… a projected CRS whose datum has a non-Greenwich prime meridian
rouault · Nov 10, 2023 · 6f9fd55 · 6f9fd55
1 parent e9c0e8c
commit 6f9fd55
Show file tree

Hide file tree

Showing 2 changed files with 64 additions and 6 deletions.
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
@@ -2886,6 +2886,16 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
 
     const auto type = proj_get_type(P);
 
+    if (type == PJ_TYPE_COMPOUND_CRS) {
+        auto ctx = P->ctx;
+        auto horiz = proj_crs_get_sub_crs(ctx, P, 0);
+        if (horiz) {
+            auto ret = proj_factors(horiz, lp);
+            proj_destroy(horiz);
+            return ret;
+        }
+    }
+
     if (type == PJ_TYPE_PROJECTED_CRS) {
         // If it is a projected CRS, then compute the factors on the conversion
         // associated to it. We need to start from a temporary geographic CRS
@@ -2899,14 +2909,34 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
         auto ctx = P->ctx;
         auto geodetic_crs = proj_get_source_crs(ctx, P);
         assert(geodetic_crs);
-        auto datum = proj_crs_get_datum(ctx, geodetic_crs);
-        auto datum_ensemble = proj_crs_get_datum_ensemble(ctx, geodetic_crs);
+        auto pm = proj_get_prime_meridian(ctx, geodetic_crs);
+        double pm_longitude = 0;
+        proj_prime_meridian_get_parameters(ctx, pm, &pm_longitude, nullptr,
+                                           nullptr);
+        proj_destroy(pm);
+        PJ *geogCRSNormalized;
         auto cs = proj_create_ellipsoidal_2D_cs(
             ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0);
-        auto geogCRSNormalized = proj_create_geographic_crs_from_datum(
-            ctx, "unnamed crs", datum ? datum : datum_ensemble, cs);
-        proj_destroy(datum);
-        proj_destroy(datum_ensemble);
+        if (pm_longitude != 0) {
+            auto ellipsoid = proj_get_ellipsoid(ctx, geodetic_crs);
+            double semi_major_metre = 0;
+            double inv_flattening = 0;
+            proj_ellipsoid_get_parameters(ctx, ellipsoid, &semi_major_metre,
+                                          nullptr, nullptr, &inv_flattening);
+            geogCRSNormalized = proj_create_geographic_crs(
+                ctx, "unname crs", "unnamed datum", proj_get_name(ellipsoid),
+                semi_major_metre, inv_flattening, "reference prime meridian", 0,
+                nullptr, 0, cs);
+            proj_destroy(ellipsoid);
+        } else {
+            auto datum = proj_crs_get_datum(ctx, geodetic_crs);
+            auto datum_ensemble =
+                proj_crs_get_datum_ensemble(ctx, geodetic_crs);
+            geogCRSNormalized = proj_create_geographic_crs_from_datum(
+                ctx, "unnamed crs", datum ? datum : datum_ensemble, cs);
+            proj_destroy(datum);
+            proj_destroy(datum_ensemble);
+        }
         proj_destroy(cs);
         auto conversion = proj_crs_get_coordoperation(ctx, P);
         auto projCS = proj_create_cartesian_2D_cs(

diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
@@ -590,6 +590,34 @@ TEST(gie, info_functions) {
         proj_destroy(P);
     }
 
+    // Test with a projected CRS whose datum has a non-Greenwich prime meridian
+    {
+        P = proj_create(PJ_DEFAULT_CTX, "EPSG:27571");
+
+        PJ_COORD c;
+        c.lp.lam = proj_torad(0.0689738);
+        c.lp.phi = proj_torad(49.508567);
+        const auto factors2 = proj_factors(P, c);
+
+        EXPECT_NEAR(factors2.meridional_scale, 1 - 12.26478760 * 1e-5, 1e-8);
+
+        proj_destroy(P);
+    }
+
+    // Test with a compound CRS of a projected CRS
+    {
+        P = proj_create(PJ_DEFAULT_CTX, "EPSG:5972");
+
+        PJ_COORD c;
+        c.lp.lam = proj_torad(10.729030600);
+        c.lp.phi = proj_torad(59.916494500);
+        const auto factors2 = proj_factors(P, c);
+
+        EXPECT_NEAR(factors2.meridional_scale, 1 - 28.54587730 * 1e-5, 1e-8);
+
+        proj_destroy(P);
+    }
+
     // Test with a geographic CRS --> error
     {
         P = proj_create(PJ_DEFAULT_CTX, "EPSG:4326");