Merge pull request #2847 from rouault/spherical_ocentric

Add support for GeodeticCRS using a Spherical ocentric coordinate system

4 files changed, 260 insertions, 16 deletions

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} <<EOF
 0d00'00.000"W 0d00'00.000"N 0.0
 79d00'00.000"W 45d00'00.000"N 0.0
@@ -156,10 +156,10 @@ $EXE +proj=latlong +datum=WGS84 +geoc \
 EOF
 #
 echo "#############################################################" >> ${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} <<EOF
 0d00'00.000"W 0d00'00.000"N 0.0
 79d00'00.000"W 44d48'27.276"N 0.000
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index 9bb85933..06aa1a2d 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -43,13 +43,13 @@ Test geocentric x/y/z consumption.
 861996.98       -4434590.01 4487348.41	79dW	45dN 0.001
 0.00    -0.00 6356752.31	0dE	90dN -0.004
 #############################################################
-Test conversion from geocentric latlong to geodetic latlong
+Test conversion from geodetic latlong to geocentric latlong
 0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 0.000
 79d00'00.000"W 45d00'00.000"N 0.0	79dW	44d48'27.276"N 0.000
 12d00'00.000"W 45d00'00.000"N 0.0	12dW	44d48'27.276"N 0.000
 0d00'00.000"W 90d00'00.000"N 0.0	0dE	90dN 0.000
 #############################################################
-Test conversion from geodetic latlong to geocentric latlong
+Test conversion from geocentric latlong to geodetic latlong
 0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 0.000
 79d00'00.000"W 44d48'27.276"N 0.000	79dW	45dN 0.000
 12d00'00.000"W 44d48'27.276"N 0.0	12dW	45dN 0.000
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");
+}