proj_factors(): make it work with projected CRS with non-metre unit and/or northing/easting axis order (fixes OSGeo#3824)

Author: rouault

src/4D_api.cpp

@@ -2894,6 +2894,9 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
         // using the same datum as the one of the projected CRS, and with
         // input coordinates being in longitude, latitude order in radian,
         // to be consistent with the expectations of the lp input parameter.
+        // We also need to create a modified projected CRS with a normalized
+        // easting/northing axis order in metre, so the resulting operation is
+        // just a single step pipeline with no axisswap or unitconvert steps.
 
         auto ctx = P->ctx;
         auto geodetic_crs = proj_get_source_crs(ctx, P);
@@ -2902,16 +2905,27 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
         auto datum_ensemble = proj_crs_get_datum_ensemble(ctx, geodetic_crs);
         auto cs = proj_create_ellipsoidal_2D_cs(
             ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0);
-        auto temp = proj_create_geographic_crs_from_datum(
+        auto 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(
+            ctx, PJ_CART2D_EASTING_NORTHING, "metre", 1.0);
+        auto projCRSNormalized = proj_create_projected_crs(
+            ctx, nullptr, geodetic_crs, conversion, projCS);
+        assert(projCRSNormalized);
         proj_destroy(geodetic_crs);
-        auto newOp =
-            proj_create_crs_to_crs_from_pj(ctx, temp, P, nullptr, nullptr);
-        proj_destroy(temp);
+        proj_destroy(conversion);
+        proj_destroy(projCS);
+        auto newOp = proj_create_crs_to_crs_from_pj(
+            ctx, geogCRSNormalized, projCRSNormalized, nullptr, nullptr);
+        proj_destroy(geogCRSNormalized);
+        proj_destroy(projCRSNormalized);
         assert(newOp);
+        // For debugging:
+        // printf("%s\n", proj_as_proj_string(ctx, newOp, PJ_PROJ_5, nullptr));
         auto ret = proj_factors(newOp, lp);
         proj_destroy(newOp);
         return ret;

test/unit/gie_self_tests.cpp

@@ -522,6 +522,72 @@ TEST(gie, info_functions) {
         proj_destroy(P);
     }
 
+    // Test with a projected CRS with feet unit
+    {
+        PJ_COORD c;
+        c.lp.lam = proj_torad(-110);
+        c.lp.phi = proj_torad(30);
+
+        P = proj_create(PJ_DEFAULT_CTX,
+                        "+proj=tmerc +lat_0=31 +lon_0=-110.166666666667 "
+                        "+k=0.9999 +x_0=213360 +y_0=0 +ellps=GRS80 +units=ft");
+        factors = proj_factors(P, c);
+        EXPECT_NEAR(factors.meridional_scale, 0.99990319, 1e-8)
+            << factors.meridional_scale;
+        EXPECT_NEAR(factors.parallel_scale, 0.99990319, 1e-8)
+            << factors.parallel_scale;
+        EXPECT_NEAR(factors.angular_distortion, 0, 1e-7)
+            << factors.angular_distortion;
+        EXPECT_NEAR(factors.meridian_parallel_angle, M_PI_2, 1e-7)
+            << factors.meridian_parallel_angle;
+        proj_destroy(P);
+
+        P = proj_create(PJ_DEFAULT_CTX, "EPSG:2222");
+
+        const auto factors2 = proj_factors(P, c);
+
+        EXPECT_NEAR(factors.meridional_scale, factors2.meridional_scale, 1e-10);
+        EXPECT_NEAR(factors.parallel_scale, factors2.parallel_scale, 1e-10);
+        EXPECT_NEAR(factors.angular_distortion, factors2.angular_distortion,
+                    1e-10);
+        EXPECT_NEAR(factors.meridian_parallel_angle,
+                    factors2.meridian_parallel_angle, 1e-109);
+
+        proj_destroy(P);
+    }
+
+    // Test with a projected CRS with northing, easting axis order
+    {
+        PJ_COORD c;
+        c.lp.lam = proj_torad(9);
+        c.lp.phi = proj_torad(0);
+
+        P = proj_create(PJ_DEFAULT_CTX, "+proj=utm +zone=32 +ellps=GRS80");
+        factors = proj_factors(P, c);
+        EXPECT_NEAR(factors.meridional_scale, 0.9996, 1e-8)
+            << factors.meridional_scale;
+        EXPECT_NEAR(factors.parallel_scale, 0.9996, 1e-8)
+            << factors.parallel_scale;
+        EXPECT_NEAR(factors.angular_distortion, 0, 1e-7)
+            << factors.angular_distortion;
+        EXPECT_NEAR(factors.meridian_parallel_angle, M_PI_2, 1e-7)
+            << factors.meridian_parallel_angle;
+        proj_destroy(P);
+
+        P = proj_create(PJ_DEFAULT_CTX, "EPSG:3044");
+
+        const auto factors2 = proj_factors(P, c);
+
+        EXPECT_NEAR(factors.meridional_scale, factors2.meridional_scale, 1e-10);
+        EXPECT_NEAR(factors.parallel_scale, factors2.parallel_scale, 1e-10);
+        EXPECT_NEAR(factors.angular_distortion, factors2.angular_distortion,
+                    1e-10);
+        EXPECT_NEAR(factors.meridian_parallel_angle,
+                    factors2.meridian_parallel_angle, 1e-109);
+
+        proj_destroy(P);
+    }
+
     // Test with a geographic CRS --> error
     {
         P = proj_create(PJ_DEFAULT_CTX, "EPSG:4326");