From 5527b10ed140e20fac8e183317514fd59e4c8b99 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 7 Sep 2021 19:57:36 +0200 Subject: Support importing/exporting WKT & PROJJSON of 2D axis spherical planetocentric geodetic CRS --- test/unit/test_io.cpp | 84 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) (limited to 'test') diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 72ac1868..96e3a849 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -894,6 +894,37 @@ TEST(wkt_parse, wkt2_EPSG_4979) { // --------------------------------------------------------------------------- +TEST(wkt_parse, wkt2_spherical_planetocentric) { + const auto wkt = + "GEODCRS[\"Mercury (2015) / Ocentric\",\n" + " DATUM[\"Mercury (2015)\",\n" + " ELLIPSOID[\"Mercury (2015)\",2440530,1075.12334801762,\n" + " LENGTHUNIT[\"metre\",1]],\n" + " ANCHOR[\"Hun Kal: 20.0\"]],\n" + " PRIMEM[\"Reference Meridian\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " CS[spherical,2],\n" + " AXIS[\"planetocentric latitude (U)\",north,\n" + " ORDER[1],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " AXIS[\"planetocentric longitude (V)\",east,\n" + " ORDER[2],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " ID[\"IAU\",19902,2015],\n" + " REMARK[\"Source of IAU Coordinate systems: " + "doi://10.1007/s10569-017-9805-5\"]]"; + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_TRUE(crs->isSphericalPlanetocentric()); + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT2_2019).get()), + wkt); +} + +// --------------------------------------------------------------------------- + static void checkGeocentric(GeodeticCRSPtr crs) { // Explicitly check this is NOT a GeographicCRS EXPECT_TRUE(!dynamic_cast(crs.get())); @@ -11476,6 +11507,59 @@ TEST(json_import, geographic_crs) { // --------------------------------------------------------------------------- +TEST(json_import, spherical_planetocentric) { + const auto json = "{\n" + " \"$schema\": \"foo\",\n" + " \"type\": \"GeodeticCRS\",\n" + " \"name\": \"Mercury (2015) / Ocentric\",\n" + " \"datum\": {\n" + " \"type\": \"GeodeticReferenceFrame\",\n" + " \"name\": \"Mercury (2015)\",\n" + " \"anchor\": \"Hun Kal: 20.0\",\n" + " \"ellipsoid\": {\n" + " \"name\": \"Mercury (2015)\",\n" + " \"semi_major_axis\": 2440530,\n" + " \"inverse_flattening\": 1075.12334801762\n" + " },\n" + " \"prime_meridian\": {\n" + " \"name\": \"Reference Meridian\",\n" + " \"longitude\": 0\n" + " }\n" + " },\n" + " \"coordinate_system\": {\n" + " \"subtype\": \"spherical\",\n" + " \"axis\": [\n" + " {\n" + " \"name\": \"Planetocentric latitude\",\n" + " \"abbreviation\": \"U\",\n" + " \"direction\": \"north\",\n" + " \"unit\": \"degree\"\n" + " },\n" + " {\n" + " \"name\": \"Planetocentric longitude\",\n" + " \"abbreviation\": \"V\",\n" + " \"direction\": \"east\",\n" + " \"unit\": \"degree\"\n" + " }\n" + " ]\n" + " },\n" + " \"id\": {\n" + " \"authority\": \"IAU\",\n" + " \"code\": 19902\n" + " },\n" + " \"remarks\": \"Source of IAU Coordinate systems: " + "doi://10.1007/s10569-017-9805-5\"\n" + "}"; + auto obj = createFromUserInput(json, nullptr); + auto gcrs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(gcrs != nullptr); + EXPECT_TRUE(gcrs->isSphericalPlanetocentric()); + EXPECT_EQ(gcrs->exportToJSON(&(JSONFormatter::create()->setSchema("foo"))), + json); +} + +// --------------------------------------------------------------------------- + TEST(json_import, geographic_crs_errors) { EXPECT_THROW( createFromUserInput( -- cgit v1.2.3 From 85733181ee7c2777139f5d1db94f2beabb737e96 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 8 Sep 2021 12:50:37 +0200 Subject: createOperations(): deal with spherical planetocentric geodetic CRS This also fixes conversion between geocentric latlong and geodetic latlong with cs2cs. This was dealt with in PR 1093, but in the wrong direction (the geocentric latitude must be <= in absolute value to the geodetic one) The issue here was linked to the semantics of the +geoc specifier, which affects the semantics of the input coordinates in the forward direction (+geoc means that the input coordinate is is a geocentric latitude), which defeats the logic of doing A to B by using the inverse path of A and the forward path of B. --- test/cli/testvarious | 12 ++-- test/cli/tv_out.dist | 4 +- test/unit/test_io.cpp | 27 +++++++- test/unit/test_operationfactory.cpp | 126 ++++++++++++++++++++++++++++++++++-- 4 files changed, 153 insertions(+), 16 deletions(-) (limited to 'test') diff --git a/test/cli/testvarious b/test/cli/testvarious index 6e9cd43c..80340369 100755 --- a/test/cli/testvarious +++ b/test/cli/testvarious @@ -144,10 +144,10 @@ $EXE +proj=geocent +datum=WGS84 \ EOF # echo "#############################################################" >> ${OUT} -echo Test conversion from geocentric latlong to geodetic latlong >> ${OUT} +echo Test conversion from geodetic latlong to geocentric latlong >> ${OUT} # -$EXE +proj=latlong +datum=WGS84 +geoc \ - +to +proj=latlong +datum=WGS84 \ +$EXE +proj=latlong +datum=WGS84 \ + +to +proj=latlong +datum=WGS84 +geoc \ -E >>${OUT} <> ${OUT} -echo Test conversion from geodetic latlong to geocentric latlong >> ${OUT} +echo Test conversion from geocentric latlong to geodetic latlong >> ${OUT} # -$EXE +proj=latlong +datum=WGS84 \ - +to +proj=latlong +datum=WGS84 +geoc \ +$EXE +proj=latlong +datum=WGS84 +geoc \ + +to +proj=latlong +datum=WGS84 \ -E >>${OUT} <(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_TRUE(crs->isSphericalPlanetocentric()); +#if 1 + EXPECT_THROW( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + FormattingException); +#else + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + input); +#endif +} + +// --------------------------------------------------------------------------- + TEST(io, projparse_projected_wktext) { std::string input("+proj=merc +foo +wktext +type=crs"); auto obj = PROJStringParser().createFromPROJString(input); @@ -10209,13 +10232,13 @@ TEST(io, projparse_init) { } { - auto obj = createFromUserInput("+title=mytitle +geoc +init=epsg:4326", + auto obj = createFromUserInput("+title=mytitle +over +init=epsg:4326", dbContext, true); auto crs = nn_dynamic_pointer_cast(obj); ASSERT_TRUE(crs != nullptr); EXPECT_EQ(crs->nameStr(), "mytitle"); EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=longlat +geoc +datum=WGS84 +no_defs +type=crs"); + "+proj=longlat +over +datum=WGS84 +no_defs +type=crs"); } { diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 56d9f743..ec7d8700 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -6092,7 +6092,7 @@ TEST(operation, createOperation_on_crs_with_canonical_bound_crs) { TEST(operation, createOperation_fallback_to_proj4_strings) { auto objDest = PROJStringParser().createFromPROJString( - "+proj=longlat +geoc +datum=WGS84 +type=crs"); + "+proj=longlat +over +datum=WGS84 +type=crs"); auto dest = nn_dynamic_pointer_cast(objDest); ASSERT_TRUE(dest != nullptr); @@ -6102,7 +6102,7 @@ TEST(operation, createOperation_fallback_to_proj4_strings) { EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline +step +proj=axisswap +order=2,1 " "+step +proj=unitconvert +xy_in=deg +xy_out=rad " - "+step +proj=longlat +geoc +datum=WGS84 " + "+step +proj=longlat +over +datum=WGS84 " "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); } @@ -6110,12 +6110,12 @@ TEST(operation, createOperation_fallback_to_proj4_strings) { TEST(operation, createOperation_fallback_to_proj4_strings_bound_of_geog) { auto objSrc = PROJStringParser().createFromPROJString( - "+proj=longlat +geoc +ellps=GRS80 +towgs84=0,0,0 +type=crs"); + "+proj=longlat +over +ellps=GRS80 +towgs84=0,0,0 +type=crs"); auto src = nn_dynamic_pointer_cast(objSrc); ASSERT_TRUE(src != nullptr); auto objDest = PROJStringParser().createFromPROJString( - "+proj=longlat +geoc +ellps=clrk66 +towgs84=0,0,0 +type=crs"); + "+proj=longlat +over +ellps=clrk66 +towgs84=0,0,0 +type=crs"); auto dest = nn_dynamic_pointer_cast(objDest); ASSERT_TRUE(dest != nullptr); @@ -6125,8 +6125,8 @@ TEST(operation, createOperation_fallback_to_proj4_strings_bound_of_geog) { EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline " "+step +proj=unitconvert +xy_in=deg +xy_out=rad " - "+step +inv +proj=longlat +geoc +ellps=GRS80 +towgs84=0,0,0 " - "+step +proj=longlat +geoc +ellps=clrk66 +towgs84=0,0,0 " + "+step +inv +proj=longlat +over +ellps=GRS80 +towgs84=0,0,0 " + "+step +proj=longlat +over +ellps=clrk66 +towgs84=0,0,0 " "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); } @@ -6832,3 +6832,117 @@ TEST(operation, derivedGeographicCRS_with_to_wgs84_to_geographicCRS) { "+proj=pipeline +step +proj=axisswap +order=2,1 " + pipeline); } } + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_spherical_ocentric_to_geographic) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), GeographicCRS::EPSG_4326); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=WGS84 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_geographic_to_spherical_ocentric) { + auto objDest = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + GeographicCRS::EPSG_4326, NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=geoc +ellps=WGS84 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_bound_of_spherical_ocentric_to_same_type) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +ellps=GRS80 +towgs84=1,2,3 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDest = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +ellps=clrk66 +towgs84=4,5,6 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=GRS80 " + "+step +proj=push +v_3 " + "+step +proj=cart +ellps=GRS80 " + "+step +proj=helmert +x=-3 +y=-3 +z=-3 " + "+step +inv +proj=cart +ellps=clrk66 " + "+step +proj=pop +v_3 " + "+step +proj=geoc +ellps=clrk66 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_spherical_ocentric_to_projected) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDest = PROJStringParser().createFromPROJString( + "+proj=utm +zone=11 +datum=WGS84 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=WGS84 " + "+step +proj=utm +zone=11 +ellps=WGS84"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, + createOperation_spherical_ocentric_to_projected_of_spherical_ocentric) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDest = PROJStringParser().createFromPROJString( + "+proj=utm +geoc +zone=11 +datum=WGS84 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=WGS84 " + "+step +proj=utm +zone=11 +ellps=WGS84"); +} -- cgit v1.2.3 From 078952e7f078e029d66ab6ca1ed594dfecadd1fc Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 8 Sep 2021 14:34:50 +0200 Subject: createOperations(): use an explicit conversion operation for geodetic <--> geocentric latitude --- test/unit/test_operationfactory.cpp | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) (limited to 'test') diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index ec7d8700..b50542c6 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -6873,6 +6873,29 @@ TEST(operation, createOperation_geographic_to_spherical_ocentric) { // --------------------------------------------------------------------------- +TEST(operation, createOperation_spherical_ocentric_to_geocentric) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDest = PROJStringParser().createFromPROJString( + "+proj=geocent +datum=WGS84 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=WGS84 " + "+step +proj=cart +ellps=WGS84"); +} + +// --------------------------------------------------------------------------- + TEST(operation, createOperation_bound_of_spherical_ocentric_to_same_type) { auto objSrc = PROJStringParser().createFromPROJString( "+proj=longlat +geoc +ellps=GRS80 +towgs84=1,2,3 +type=crs"); -- cgit v1.2.3