Merge pull request #2847 from rouault/spherical_ocentric

Add support for GeodeticCRS using a Spherical ocentric coordinate system

2 files changed, 252 insertions, 8 deletions

diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 72ac1868..791a8b9e 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<GeodeticCRS>(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<GeographicCRS *>(crs.get()));
@@ -10032,6 +10063,29 @@ TEST(io, projparse_geocent_wktext) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_geoc) {
+    std::string input("+proj=longlat +geoc +datum=WGS84 +no_defs +type=crs");
+    auto obj = PROJStringParser().createFromPROJString(input);
+    auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(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);
@@ -10178,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<GeographicCRS>(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");
     }
 
     {
@@ -11476,6 +11530,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<GeodeticCRS>(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(
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index 56d9f743..b50542c6 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<GeographicCRS>(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<BoundCRS>(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<BoundCRS>(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,140 @@ 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<GeodeticCRS>(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<GeodeticCRS>(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_spherical_ocentric_to_geocentric) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=longlat +geoc +datum=WGS84 +type=crs");
+    auto src = nn_dynamic_pointer_cast<GeodeticCRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDest = PROJStringParser().createFromPROJString(
+        "+proj=geocent +datum=WGS84 +type=crs");
+    auto dest = nn_dynamic_pointer_cast<GeodeticCRS>(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");
+    auto src = nn_dynamic_pointer_cast<BoundCRS>(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<BoundCRS>(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<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDest = PROJStringParser().createFromPROJString(
+        "+proj=utm +zone=11 +datum=WGS84 +type=crs");
+    auto dest = nn_dynamic_pointer_cast<CRS>(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<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDest = PROJStringParser().createFromPROJString(
+        "+proj=utm +geoc +zone=11 +datum=WGS84 +type=crs");
+    auto dest = nn_dynamic_pointer_cast<CRS>(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");
+}