diff options
| -rw-r--r-- | docs/source/development/reference/functions.rst | 9 | ||||
| -rw-r--r-- | src/4D_api.cpp | 42 | ||||
| -rw-r--r-- | test/unit/gie_self_tests.cpp | 27 |
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) |