Merge pull request #2868 from rouault/proj_factors_with_projected_crs

proj_factors(): accept P to be a projected CRS (fixes #2854)

3 files changed, 78 insertions, 0 deletions

diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst
index 10bb07d0..bf5bb94b 100644
--- a/docs/source/development/reference/functions.rst
+++ b/docs/source/development/reference/functions.rst
@@ -756,6 +756,15 @@ Various
     distortion and meridian convergence. Depending on the underlying projection
     values will be calculated either numerically (default) or analytically.
 
+    Starting with PROJ 8.2, the P object can be a projected CRS, for example
+    instantiated from a EPSG CRS code. The factors computed will be those of the
+    map projection implied by the transformation from the base geographic CRS of
+    the projected CRS to the projected CRS.
+
+    The input geodetic coordinate lp should be such that lp.lam is the longitude
+    in radian, and lp.phi the latitude in radian (thus independently of the
+    definition of the base CRS, if P is a projected CRS).
+
     The function also calculates the partial derivatives of the given
     coordinate.
 
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 9c226647..132c21d3 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1896,6 +1896,48 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
     if (nullptr==P)
         return factors;
 
+    const auto type = proj_get_type(P);
+
+    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
+        // 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.
+
+        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 cs = proj_create_ellipsoidal_2D_cs(
+            ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0);
+        auto temp = proj_create_geographic_crs_from_datum(
+            ctx, "unnamed crs", datum ? datum : datum_ensemble,
+            cs);
+        proj_destroy(datum);
+        proj_destroy(datum_ensemble);
+        proj_destroy(cs);
+        proj_destroy(geodetic_crs);
+        auto newOp = proj_create_crs_to_crs_from_pj(ctx, temp, P, nullptr, nullptr);
+        proj_destroy(temp);
+        assert(newOp);
+        auto ret = proj_factors(newOp, lp);
+        proj_destroy(newOp);
+        return ret;
+    }
+
+    if( type != PJ_TYPE_CONVERSION &&
+        type != PJ_TYPE_TRANSFORMATION &&
+        type != PJ_TYPE_CONCATENATED_OPERATION &&
+        type != PJ_TYPE_OTHER_COORDINATE_OPERATION )
+    {
+        proj_log_error(P, _("Invalid type for P object"));
+        proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+        return factors;
+    }
+
     if (pj_factors(lp.lp, P, 0.0, &f))
         return factors;
 
diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
index 5fc8f0f1..5f96a93a 100644
--- a/test/unit/gie_self_tests.cpp
+++ b/test/unit/gie_self_tests.cpp
@@ -475,6 +475,33 @@ TEST(gie, info_functions) {
 
     proj_destroy(P);
 
+    // Test with a projected CRS
+    {
+        P = proj_create(PJ_DEFAULT_CTX, "EPSG:3395");
+
+        const auto factors2 = proj_factors(P, a);
+
+        EXPECT_EQ(factors.angular_distortion, factors2.angular_distortion);
+        EXPECT_EQ(factors.meridian_parallel_angle,
+                  factors2.meridian_parallel_angle);
+        EXPECT_EQ(factors.meridian_convergence, factors2.meridian_convergence);
+        EXPECT_EQ(factors.tissot_semimajor, factors2.tissot_semimajor);
+
+        proj_destroy(P);
+    }
+
+    // Test with a geographic CRS --> error
+    {
+        P = proj_create(PJ_DEFAULT_CTX, "EPSG:4326");
+
+        const auto factors2 = proj_factors(P, a);
+        EXPECT_NE(proj_errno(P), 0);
+        proj_errno_reset(P);
+        EXPECT_EQ(factors2.meridian_parallel_angle, 0);
+
+        proj_destroy(P);
+    }
+
     /* Check that proj_list_* functions work by looping through them */
     size_t n = 0;
     for (oper_list = proj_list_operations(); oper_list->id; ++oper_list)