From aeb024f2df5fb1c4d83959b6b9edcb62e8385226 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 17 Sep 2020 18:34:06 +0200 Subject: Adjust createBoundCRSToWGS84IfPossible() and operation filtering (for POSGAR 2007 to WGS84 issues) (fixes #2356) - We make createBoundCRSToWGS84IfPossible() more restrictive. If there are more than one Helmert transformation from the CRS to WGS 84 covering the area of use of the CRS, we do not create a BoundCRS / +towgs84 - In createOperations() filtering, we are less aggressive in discarding operations that have the same area of use but worse accuracy. We do it only if they involve more transformation steps. We now get: ``` $ projinfo EPSG:5340 -o PROJ PROJ.4 string: +proj=longlat +ellps=GRS80 +no_defs +type=crs $ projinfo -s EPSG:5340 -t EPSG:4326 --spatial-test intersects --summary Candidate operations found: 2 EPSG:9264, POSGAR 2007 to WGS 84 (2), 0.5 m, Argentina EPSG:5351, POSGAR 2007 to WGS 84 (1), 1.0 m, Argentina ``` --- test/cli/proj_outIGNF.dist | 4 ++-- test/cli/testprojinfo_out.dist | 8 ++++++-- test/unit/test_c_api.cpp | 15 ++++++--------- test/unit/test_crs.cpp | 39 +++++++++++++-------------------------- test/unit/test_operation.cpp | 31 ++++++------------------------- 5 files changed, 33 insertions(+), 64 deletions(-) (limited to 'test') diff --git a/test/cli/proj_outIGNF.dist b/test/cli/proj_outIGNF.dist index c12b883b..35d053fa 100644 --- a/test/cli/proj_outIGNF.dist +++ b/test/cli/proj_outIGNF.dist @@ -8,7 +8,7 @@ 311552.5340 1906457.4840 0.0000 358799.172 6342652.486 0.000 960488.4138 1910172.8812 0.0000 1007068.686 6340907.237 0.000 600000.0000 1699510.8340 0.0000 645204.279 6133556.746 0.000 -1203792.5981 626873.17210 0.0000 * * inf +1203792.5981 626873.17210 0.0000 1238837.253 5057451.037 0.000 +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX 600000.0000 2600545.4523 0.0000 179040.148 5610495.275 0.000 135638.3592 2418760.4094 0.0000 -303729.363 5410118.356 0.000 @@ -17,7 +17,7 @@ 311552.5340 1906457.4840 0.0000 -96825.465 4909184.136 0.000 960488.4138 1910172.8812 0.0000 523880.019 4909191.141 0.000 600000.0000 1699510.8340 0.0000 179047.633 4708817.007 0.000 -1203792.5981 626873.17210 0.0000 * * inf +1203792.5981 626873.17210 0.0000 658259.467 3623786.764 0.000 +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX 2d20'11.4239243" 50d23'59.7718445" 0.0 179040.151 5610495.281 0.000 -3d57'49.4051448" 48d35'59.7121716" 0.0 -303729.365 5410118.352 0.000 diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist index 8ece2e62..5e1f1dd5 100644 --- a/test/cli/testprojinfo_out.dist +++ b/test/cli/testprojinfo_out.dist @@ -1044,8 +1044,9 @@ Testing -s +proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs -t EPSG:432 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=@foo.gtx +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 Testing -s "GDA94" -t "WGS 84 (G1762)" --spatial-test intersects --summary. Should include transformations through ITRF2008 and GDA2020 -Candidate operations found: 6 +Candidate operations found: 7 unknown id, GDA94 to GDA2020 (1) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.21 m, Australia - GDA +unknown id, GDA94 to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 5 m, Australia - GDA unknown id, Conversion from GDA94 (geog2D) to GDA94 (geocentric) + Inverse of ITRF2008 to GDA94 (1) + Inverse of WGS 84 (G1762) to ITRF2008 (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.04 m, Australia - onshore and EEZ unknown id, GDA94 to GDA2020 (2) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore unknown id, GDA94 to GDA2020 (3) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore @@ -1053,16 +1054,19 @@ unknown id, GDA94 to GDA2020 (5) + Conversion from GDA2020 (geog2D) to GDA2020 ( unknown id, GDA94 to GDA2020 (4) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Christmas Island - onshore Testing -s "AGD66" -t "WGS 84 (G1762)" --spatial-test intersects --summary. Should include a transformation through GDA2020 -Candidate operations found: 11 +Candidate operations found: 14 unknown id, AGD66 to WGS 84 (18) + WGS 84 to WGS 84 (G1762), 5 m, Australia - offshore unknown id, AGD66 to WGS 84 (16) + WGS 84 to WGS 84 (G1762), 7 m, Australia - onshore +unknown id, AGD66 to WGS 84 (20) + WGS 84 to WGS 84 (G1762), 11 m, Australia - onshore unknown id, AGD66 to WGS 84 (15) + WGS 84 to WGS 84 (G1762), 3 m, Australia - Northern Territory unknown id, AGD66 to WGS 84 (13) + WGS 84 to WGS 84 (G1762), 3 m, Australia - New South Wales and Victoria unknown id, AGD66 to WGS 84 (21) + WGS 84 to WGS 84 (G1762), 7 m, Papua New Guinea - mainland onshore unknown id, AGD66 to WGS 84 (14) + WGS 84 to WGS 84 (G1762), 3 m, Australia - Tasmania unknown id, AGD66 to WGS 84 (19) + WGS 84 to WGS 84 (G1762), 4 m, Papua New Guinea - PFTB +unknown id, AGD66 to WGS 84 (22) + WGS 84 to WGS 84 (G1762), 6 m, Papua New Guinea - PFTB unknown id, AGD66 to WGS 84 (23) + WGS 84 to WGS 84 (G1762), 6 m, Papua New Guinea - North Fly unknown id, AGD66 to GDA2020 (1) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - Australian Capital Territory +unknown id, AGD66 to WGS 84 (12) + WGS 84 to WGS 84 (G1762), 3 m, Australia - Australian Capital Territory unknown id, AGD66 to WGS 84 (17) + WGS 84 to WGS 84 (G1762), 3 m, Australia - onshore unknown id, Ballpark geographic offset from AGD66 to WGS 84 (G1762), unknown accuracy, World, has ballpark transformation diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index b9a6cdd5..564833db 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -540,7 +540,7 @@ TEST_F(CApi, proj_as_proj_string_approx_tmerc_option_yes) { // --------------------------------------------------------------------------- TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) { - auto crs = proj_create_from_database(m_ctxt, "EPSG", "3844", + auto crs = proj_create_from_database(m_ctxt, "EPSG", "4807", PJ_CATEGORY_CRS, false, nullptr); ObjectKeeper keeper(crs); ASSERT_NE(crs, nullptr); @@ -552,10 +552,8 @@ TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) { auto proj_4 = proj_as_proj_string(m_ctxt, res, PJ_PROJ_4, nullptr); ASSERT_NE(proj_4, nullptr); EXPECT_EQ(std::string(proj_4), - "+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 " - "+y_0=500000 +ellps=krass " - "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69 " - "+units=m +no_defs +type=crs"); + "+proj=longlat +ellps=clrk80ign +pm=paris " + "+towgs84=-168,-60,320,0,0,0,0 +no_defs +type=crs"); auto base_crs = proj_get_source_crs(m_ctxt, res); ObjectKeeper keeper_base_crs(base_crs); @@ -572,8 +570,7 @@ TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) { std::vector values(7, 0); EXPECT_TRUE(proj_coordoperation_get_towgs84_values(m_ctxt, transf, values.data(), 7, true)); - auto expected = std::vector{2.329, -147.042, -92.08, -0.309, - 0.325, 0.497, 5.69}; + auto expected = std::vector{-168, -60, 320, 0, 0, 0, 0}; EXPECT_EQ(values, expected); auto res2 = proj_crs_create_bound_crs(m_ctxt, base_crs, hub_crs, transf); @@ -1485,7 +1482,7 @@ TEST_F(CApi, proj_create_operations_discard_superseded) { ASSERT_NE(res, nullptr); ObjListKeeper keeper_res(res); - EXPECT_EQ(proj_list_get_count(res), 2); + EXPECT_EQ(proj_list_get_count(res), 4); } // --------------------------------------------------------------------------- @@ -1594,7 +1591,7 @@ TEST_F(CApi, proj_create_operations_with_pivot) { auto res = proj_create_operations(m_ctxt, source_crs, target_crs, ctxt); ASSERT_NE(res, nullptr); ObjListKeeper keeper_res(res); - EXPECT_EQ(proj_list_get_count(res), 7); + EXPECT_EQ(proj_list_get_count(res), 8); auto op = proj_list_get(m_ctxt, res, 0); ASSERT_NE(op, nullptr); ObjectKeeper keeper_op(op); diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp index caee0646..0e340560 100644 --- a/test/unit/test_crs.cpp +++ b/test/unit/test_crs.cpp @@ -5422,40 +5422,18 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) { { // Pulkovo 42 Romania auto crs_3844 = factory->createCoordinateReferenceSystem("3844"); - auto bound = crs_3844->createBoundCRSToWGS84IfPossible( - dbContext, CoordinateOperationContext::IntermediateCRSUse::NEVER); - EXPECT_NE(bound, crs_3844); - EXPECT_EQ(bound->createBoundCRSToWGS84IfPossible( + EXPECT_EQ(crs_3844->createBoundCRSToWGS84IfPossible( dbContext, CoordinateOperationContext::IntermediateCRSUse::NEVER), - bound); - auto boundCRS = nn_dynamic_pointer_cast(bound); - ASSERT_TRUE(boundCRS != nullptr); - EXPECT_EQ( - boundCRS->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 " - "+y_0=500000 +ellps=krass " - "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69 " - "+units=m +no_defs +type=crs"); + crs_3844); } { // Pulkovo 42 Poland auto crs_2171 = factory->createCoordinateReferenceSystem("2171"); - auto bound = crs_2171->createBoundCRSToWGS84IfPossible( - dbContext, CoordinateOperationContext::IntermediateCRSUse::NEVER); - EXPECT_NE(bound, crs_2171); - EXPECT_EQ(bound->createBoundCRSToWGS84IfPossible( + EXPECT_EQ(crs_2171->createBoundCRSToWGS84IfPossible( dbContext, CoordinateOperationContext::IntermediateCRSUse::NEVER), - bound); - auto boundCRS = nn_dynamic_pointer_cast(bound); - ASSERT_TRUE(boundCRS != nullptr); - EXPECT_EQ( - boundCRS->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=sterea +lat_0=50.625 +lon_0=21.0833333333333 " - "+k=0.9998 +x_0=4637000 +y_0=5647000 +ellps=krass " - "+towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 " - "+units=m +no_defs +type=crs"); + crs_2171); } { // NTF (Paris) @@ -5547,6 +5525,15 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) { ::testing::UnitTest::GetInstance()->elapsed_time(); EXPECT_LE(time_after - time_before, 500); } + { + // POSGAR 2007: it has 2 helmert shifts to WGS84 (#2356). Don't take + // an arbitrary one + auto crs_5340 = factory->createCoordinateReferenceSystem("5340"); + EXPECT_EQ(crs_5340->createBoundCRSToWGS84IfPossible( + dbContext, + CoordinateOperationContext::IntermediateCRSUse::NEVER), + crs_5340); + } } // --------------------------------------------------------------------------- diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 0238638d..eca5d59b 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -4599,7 +4599,7 @@ TEST(operation, geogCRS_to_geogCRS_context_NAD27_to_WGS84) { authFactory->createCoordinateReferenceSystem("4267"), // NAD27 authFactory->createCoordinateReferenceSystem("4326"), // WGS84 ctxt); - ASSERT_EQ(list.size(), 78U); + ASSERT_EQ(list.size(), 79U); EXPECT_EQ(list[0]->nameStr(), "NAD27 to WGS 84 (33)"); // 1.0 m, Canada - NAD27 EXPECT_EQ(list[1]->nameStr(), @@ -5050,7 +5050,7 @@ TEST(operation, geogCRS_to_geogCRS_context_concatenated_operation) { authFactory->createCoordinateReferenceSystem("4807"), // NTF(Paris) authFactory->createCoordinateReferenceSystem("4171"), // RGF93 ctxt); - ASSERT_EQ(list.size(), 5U); + ASSERT_EQ(list.size(), 4U); EXPECT_EQ(list[0]->nameStr(), "NTF (Paris) to RGF93 (1)"); EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), @@ -5158,7 +5158,7 @@ TEST(operation, geogCRS_to_geogCRS_CH1903_to_CH1903plus_context) { authFactory->createCoordinateReferenceSystem("4149"), // CH1903 authFactory->createCoordinateReferenceSystem("4150"), // CH1903+ ctxt); - ASSERT_TRUE(list.size() == 2U || list.size() == 3U); + ASSERT_TRUE(list.size() == 1U); EXPECT_EQ(list[0]->nameStr(), "CH1903 to CH1903+ (1)"); EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), @@ -5167,25 +5167,6 @@ TEST(operation, geogCRS_to_geogCRS_CH1903_to_CH1903plus_context) { "+step +proj=hgridshift +grids=ch_swisstopo_CHENyx06a.tif " "+step +proj=unitconvert +xy_in=rad +xy_out=deg " "+step +proj=axisswap +order=2,1"); - - if (list.size() == 2U) { - // Grids not there - EXPECT_EQ(list[1]->nameStr(), - "CH1903 to ETRS89 (1) + Inverse of CH1903+ to ETRS89 (1)"); - EXPECT_EQ( - list[1]->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=noop"); - } else { - // Grids available - EXPECT_EQ(list[1]->nameStr(), - "CH1903 to ETRS89 (2) + Inverse of CH1903+ to ETRS89 (1)"); - - EXPECT_EQ(list[2]->nameStr(), - "CH1903 to ETRS89 (1) + Inverse of CH1903+ to ETRS89 (1)"); - EXPECT_EQ( - list[2]->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=noop"); - } } // --------------------------------------------------------------------------- @@ -6069,7 +6050,7 @@ TEST(operation, projCRS_to_projCRS_context_compatible_area) { authFactory->createCoordinateReferenceSystem( "2171"), // Pulkovo 42 Poland I ctxt); - ASSERT_EQ(list.size(), 1U); + ASSERT_EQ(list.size(), 2U); EXPECT_EQ(list[0]->nameStr(), "Inverse of UTM zone 34N + Inverse of Pulkovo 1942(58) to WGS 84 " "(1) + Poland zone I"); @@ -6088,7 +6069,7 @@ TEST(operation, projCRS_to_projCRS_context_compatible_area_bis) { "3844"), // Pulkovo 42 Stereo 70 (Romania) authFactory->createCoordinateReferenceSystem("32634"), // UTM 34 ctxt); - ASSERT_EQ(list.size(), 1U); + ASSERT_EQ(list.size(), 3U); EXPECT_EQ(list[0]->nameStr(), "Inverse of Stereo 70 + " "Pulkovo 1942(58) to WGS 84 " "(19) + UTM zone 34N"); @@ -6107,7 +6088,7 @@ TEST(operation, projCRS_to_projCRS_context_one_incompatible_area) { authFactory->createCoordinateReferenceSystem( "2171"), // Pulkovo 42 Poland I ctxt); - ASSERT_EQ(list.size(), 1U); + ASSERT_EQ(list.size(), 2U); EXPECT_EQ(list[0]->nameStr(), "Inverse of UTM zone 31N + Inverse of Pulkovo 1942(58) to WGS 84 " "(1) + Poland zone I"); -- cgit v1.2.3 From 75074ce863e36299ec39b39decdd435eed15ad63 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 17 Sep 2020 19:46:07 +0200 Subject: createOperations(): tune sorting of transformations so that the ones with greater 'version numbers' are prefered over other ones, when all other comparison criteria are equal. Helps with Amersfoort RD New to EPSG:4326 --- test/cli/testprojinfo_out.dist | 2 +- test/unit/test_operation.cpp | 18 ++++++++++++++++++ 2 files changed, 19 insertions(+), 1 deletion(-) (limited to 'test') diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist index 5e1f1dd5..5d36f42e 100644 --- a/test/cli/testprojinfo_out.dist +++ b/test/cli/testprojinfo_out.dist @@ -1048,8 +1048,8 @@ Candidate operations found: 7 unknown id, GDA94 to GDA2020 (1) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.21 m, Australia - GDA unknown id, GDA94 to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 5 m, Australia - GDA unknown id, Conversion from GDA94 (geog2D) to GDA94 (geocentric) + Inverse of ITRF2008 to GDA94 (1) + Inverse of WGS 84 (G1762) to ITRF2008 (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.04 m, Australia - onshore and EEZ -unknown id, GDA94 to GDA2020 (2) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore unknown id, GDA94 to GDA2020 (3) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore +unknown id, GDA94 to GDA2020 (2) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore unknown id, GDA94 to GDA2020 (5) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Cocos (Keeling) Islands - onshore unknown id, GDA94 to GDA2020 (4) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Christmas Island - onshore diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index eca5d59b..23f3489b 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6214,6 +6214,24 @@ TEST(operation, projCRS_to_projCRS_through_geog3D) { // --------------------------------------------------------------------------- +TEST(operation, transform_from_amersfoort_rd_new_to_epsg_4326) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + auto list = CoordinateOperationFactory::create()->createOperations( + authFactory->createCoordinateReferenceSystem("28992"), + authFactory->createCoordinateReferenceSystem("4326"), ctxt); + ASSERT_EQ(list.size(), 2U); + // The order matters: "Amersfoort to WGS 84 (4)" replaces "Amersfoort to WGS + // 84 (3)" + EXPECT_EQ(list[0]->nameStr(), + "Inverse of RD New + Amersfoort to WGS 84 (4)"); + EXPECT_EQ(list[1]->nameStr(), + "Inverse of RD New + Amersfoort to WGS 84 (3)"); +} + +// --------------------------------------------------------------------------- + TEST(operation, boundCRS_of_geogCRS_to_geogCRS) { auto boundCRS = BoundCRS::createFromTOWGS84( GeographicCRS::EPSG_4807, std::vector{1, 2, 3, 4, 5, 6, 7}); -- cgit v1.2.3