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

Updated doc:

    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).

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)