From a92526a4ba60abc9c8909e5fa098b4102a44dc7e Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 6 Oct 2021 00:49:43 +0200 Subject: proj_create_crs_to_crs() + proj_trans(): fix when non-Greenwich prime meridian is involved This fixes a regression introduced in 7af1d5741da08d9546b907e0da2c21c54c61b27 / PROJ 7.2.0 where reprojection of area of use was broken when the source/target CRS did not use Greenwich as prime meridian. Fixes https://lists.osgeo.org/pipermail/gdal-dev/2021-October/054764.html Now with the fix: - using grid: $ echo 286415 431434 | PROJ_NETWORK=ON src/cs2cs -d 4 EPSG:20790 EPSG:3763 86412.4262 131434.1706 0.0000 - not using it: $ echo 286415 431434 | src/cs2cs -d 4 EPSG:20790 EPSG:3763 86412.5265 131433.8561 0.0000 --- src/4D_api.cpp | 21 +++++++++++++++------ test/unit/test_c_api.cpp | 43 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 6 deletions(-) diff --git a/src/4D_api.cpp b/src/4D_api.cpp index 6f8c3027..4e575f14 100644 --- a/src/4D_api.cpp +++ b/src/4D_api.cpp @@ -1612,15 +1612,24 @@ static PJ* create_operation_to_geog_crs(PJ_CONTEXT* ctx, const PJ* crs) { geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_2D_CRS || geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_3D_CRS ) { - auto datum = proj_crs_get_datum(ctx, geodetic_crs); - auto datum_ensemble = proj_crs_get_datum_ensemble(ctx, geodetic_crs); + auto datum = proj_crs_get_datum_forced(ctx, geodetic_crs); + assert( datum ); auto cs = proj_create_ellipsoidal_2D_cs( ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, nullptr, 0); - auto temp = proj_create_geographic_crs_from_datum( - ctx, "unnamed crs", datum ? datum : datum_ensemble, - cs); + auto ellps = proj_get_ellipsoid(ctx, datum); proj_destroy(datum); - proj_destroy(datum_ensemble); + double semi_major_metre = 0; + double inv_flattening = 0; + proj_ellipsoid_get_parameters(ctx, ellps, &semi_major_metre, + nullptr, nullptr, &inv_flattening); + // It is critical to set the prime meridian to 0 + auto temp = proj_create_geographic_crs( + ctx, "unnamed crs", "unnamed datum", + proj_get_name(ellps), + semi_major_metre, inv_flattening, + "Reference prime meridian", 0, nullptr, 0, + cs); + proj_destroy(ellps); proj_destroy(cs); proj_destroy(geodetic_crs); geodetic_crs = temp; diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index a9315fa7..fb8e79a0 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -1646,6 +1646,49 @@ TEST_F(CApi, proj_create_operations) { // --------------------------------------------------------------------------- +TEST_F(CApi, proj_create_operations_prime_meridian_non_greenwich) { + auto ctxt = proj_create_operation_factory_context(m_ctxt, nullptr); + ASSERT_NE(ctxt, nullptr); + ContextKeeper keeper_ctxt(ctxt); + + auto source_crs = proj_create_from_database( + m_ctxt, "EPSG", "27562", PJ_CATEGORY_CRS, false, + nullptr); // "NTF (Paris) / Lambert Centre France" + ASSERT_NE(source_crs, nullptr); + ObjectKeeper keeper_source_crs(source_crs); + + auto target_crs = proj_create_from_database( + m_ctxt, "EPSG", "4258", PJ_CATEGORY_CRS, false, nullptr); // ETRS89 + ASSERT_NE(target_crs, nullptr); + ObjectKeeper keeper_target_crs(target_crs); + + proj_operation_factory_context_set_spatial_criterion( + m_ctxt, ctxt, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION); + + proj_operation_factory_context_set_grid_availability_use( + m_ctxt, ctxt, PROJ_GRID_AVAILABILITY_IGNORED); + + auto res = proj_create_operations(m_ctxt, source_crs, target_crs, ctxt); + ASSERT_NE(res, nullptr); + ObjListKeeper keeper_res(res); + + { + PJ_COORD coord; + // lat,lon=49,-4 if using grid + coord.xy.x = 136555.58288992; + coord.xy.y = 463344.51894296; + int idx = proj_get_suggested_operation(m_ctxt, res, PJ_FWD, coord); + ASSERT_GE(idx, 0); + auto op = proj_list_get(m_ctxt, res, idx); + ASSERT_NE(op, nullptr); + ObjectKeeper keeper_op(op); + // Transformation using grid + EXPECT_EQ(proj_coordoperation_get_grid_used_count(m_ctxt, op), 1); + } +} + +// --------------------------------------------------------------------------- + TEST_F(CApi, proj_get_suggested_operation_with_operations_without_area_of_use) { auto ctxt = proj_create_operation_factory_context(m_ctxt, nullptr); ASSERT_NE(ctxt, nullptr); -- cgit v1.2.3