From e04c6a6c64809b432db77c701dae9b6014a69fe8 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 10 Nov 2023 18:27:07 +0100 Subject: [PATCH 1/2] 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 --- src/4D_api.cpp | 42 ++++++++++++++++++++++++++++++------ test/unit/gie_self_tests.cpp | 28 ++++++++++++++++++++++++ 2 files changed, 64 insertions(+), 6 deletions(-) diff --git a/src/4D_api.cpp b/src/4D_api.cpp index 6451a11971..d79fec8fd1 100644 --- 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 index 646e1cb718..6650126068 100644 --- 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"); From dc01f9c4ff85b9a64b1ea530140fc976dde501c0 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 10 Nov 2023 19:37:36 +0100 Subject: [PATCH 2/2] proj: projection factors: fix when input is a compound CRS of a projected CRS, and a projected CRS whose datum has a non-Greenwich prime meridian --- src/apps/proj.cpp | 51 +++++++++++++++++++++++++++++++------- test/cli/testproj | 15 +++++++++++ test/cli/testproj_out.dist | 2 ++ 3 files changed, 59 insertions(+), 9 deletions(-) diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp index 515461993f..099c4b5dce 100644 --- a/src/apps/proj.cpp +++ b/src/apps/proj.cpp @@ -537,7 +537,20 @@ int main(int argc, char **argv) { eargc--; // logic copied from proj_factors function if (PJ *P = proj_create(nullptr, ocrs.c_str())) { - const auto type = proj_get_type(P); + auto type = proj_get_type(P); + auto ctx = P->ctx; + if (type == PJ_TYPE_COMPOUND_CRS) { + auto horiz = proj_crs_get_sub_crs(ctx, P, 0); + if (horiz) { + if (proj_get_type(horiz) == PJ_TYPE_PROJECTED_CRS) { + proj_destroy(P); + P = horiz; + type = proj_get_type(P); + } else { + proj_destroy(horiz); + } + } + } if (type == PJ_TYPE_PROJECTED_CRS) { try { auto crs = dynamic_cast( @@ -548,18 +561,38 @@ int main(int argc, char **argv) { dir == NS_PROJ::cs::AxisDirection::SOUTH; } catch (...) { } - 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); Proj = proj_create_crs_to_crs_from_pj(ctx, geogCRSNormalized, P, nullptr, nullptr); diff --git a/test/cli/testproj b/test/cli/testproj index 375f53cdbf..67f937430a 100755 --- a/test/cli/testproj +++ b/test/cli/testproj @@ -41,6 +41,21 @@ echo "Test CRS option" $EXE EPSG:32620 -S >>${OUT} <>${OUT} <>${OUT} < ${OUT}.tmp +mv ${OUT}.tmp ${OUT} # # do 'diff' with distribution results diff --git a/test/cli/testproj_out.dist b/test/cli/testproj_out.dist index 73320d0a1e..f7b6036a8a 100644 --- a/test/cli/testproj_out.dist +++ b/test/cli/testproj_out.dist @@ -1,2 +1,4 @@ 2 49 2.000 49.000 500000.00 0.00 <0.9996 0.9996 0.9992 0 0.9996 0.9996> +600000.00 1200000.00 <0.999877 0.999877 0.999755 0 0.999877 0.999877> +500000.00 0.00 <0.9996 0.9996 0.9992 0 0.9996 0.9996>