cs2cs: upgrade to use proj_create_crs_to_crs()

3 files changed, 271 insertions, 41 deletions

diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 8c9f114b..6c22cade 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -809,6 +809,22 @@ TEST_F(CApi, proj_obj_get_source_target_crs_transformation) {
 
 // ---------------------------------------------------------------------------
 
+TEST_F(CApi, proj_obj_get_source_crs_of_projected_crs) {
+    auto crs = proj_obj_create_from_wkt(
+        m_ctxt,
+        createProjectedCRS()->exportToWKT(WKTFormatter::create().get()).c_str(),
+        nullptr);
+    ASSERT_NE(crs, nullptr);
+    ObjectKeeper keeper(crs);
+
+    auto sourceCRS = proj_obj_get_source_crs(crs);
+    ASSERT_NE(sourceCRS, nullptr);
+    ObjectKeeper keeper_sourceCRS(sourceCRS);
+    EXPECT_EQ(std::string(proj_obj_get_name(sourceCRS)), "WGS 84");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST_F(CApi, proj_obj_get_source_target_crs_conversion_without_crs) {
     auto obj = proj_obj_create_from_database(
         m_ctxt, "EPSG", "16031", PJ_OBJ_CATEGORY_COORDINATE_OPERATION, false,
@@ -1955,4 +1971,69 @@ TEST_F(CApi, proj_obj_create_projections) {
     /* END: Generated by scripts/create_c_api_projections.py*/
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_obj_cs_get_axis_info) {
+    {
+        auto crs = proj_obj_create_from_database(
+            m_ctxt, "EPSG", "4326", PJ_OBJ_CATEGORY_CRS, false, nullptr);
+        ASSERT_NE(crs, nullptr);
+        ObjectKeeper keeper(crs);
+
+        auto cs = proj_obj_crs_get_coordinate_system(crs);
+        ASSERT_NE(cs, nullptr);
+        ObjectKeeper keeperCs(cs);
+
+        const char *type = proj_obj_cs_get_type(cs);
+        ASSERT_NE(type, nullptr);
+        EXPECT_EQ(std::string(type), "Ellipsoidal");
+
+        EXPECT_EQ(proj_obj_cs_get_axis_count(cs), 2);
+
+        EXPECT_FALSE(proj_obj_cs_get_axis_info(cs, -1, nullptr, nullptr,
+                                               nullptr, nullptr, nullptr));
+
+        EXPECT_FALSE(proj_obj_cs_get_axis_info(cs, 2, nullptr, nullptr, nullptr,
+                                               nullptr, nullptr));
+
+        EXPECT_TRUE(proj_obj_cs_get_axis_info(cs, 0, nullptr, nullptr, nullptr,
+                                              nullptr, nullptr));
+
+        const char *name = nullptr;
+        const char *abbrev = nullptr;
+        const char *direction = nullptr;
+        double unitConvFactor = 0.0;
+        const char *unitName = nullptr;
+
+        EXPECT_TRUE(proj_obj_cs_get_axis_info(cs, 0, &name, &abbrev, &direction,
+                                              &unitConvFactor, &unitName));
+        ASSERT_NE(name, nullptr);
+        ASSERT_NE(abbrev, nullptr);
+        ASSERT_NE(direction, nullptr);
+        ASSERT_NE(unitName, nullptr);
+        EXPECT_EQ(std::string(name), "Geodetic latitude");
+        EXPECT_EQ(std::string(abbrev), "Lat");
+        EXPECT_EQ(std::string(direction), "north");
+        EXPECT_EQ(unitConvFactor, 0.017453292519943295) << unitConvFactor;
+        EXPECT_EQ(std::string(unitName), "degree");
+    }
+
+    // Non CRS object
+    {
+        auto obj = proj_obj_create_from_database(
+            m_ctxt, "EPSG", "1170", PJ_OBJ_CATEGORY_COORDINATE_OPERATION, false,
+            nullptr);
+        ASSERT_NE(obj, nullptr);
+        ObjectKeeper keeper(obj);
+        EXPECT_EQ(proj_obj_crs_get_coordinate_system(obj), nullptr);
+
+        EXPECT_EQ(proj_obj_cs_get_type(obj), nullptr);
+
+        EXPECT_EQ(proj_obj_cs_get_axis_count(obj), -1);
+
+        EXPECT_FALSE(proj_obj_cs_get_axis_info(obj, 0, nullptr, nullptr,
+                                               nullptr, nullptr, nullptr));
+    }
+}
+
 } // namespace
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 99f58739..a4974234 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -5686,23 +5686,33 @@ TEST(io, projparse_longlat_a_f_non_zero) {
         PROJStringParser().createFromPROJString("+proj=longlat +a=2 +f=0.5");
     auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
     ASSERT_TRUE(crs != nullptr);
-    WKTFormatterNNPtr f(WKTFormatter::create());
-    f->simulCurNodeHasId();
-    crs->exportToWKT(f.get());
-    auto expected = "GEODCRS[\"unknown\",\n"
-                    "    DATUM[\"unknown\",\n"
-                    "        ELLIPSOID[\"unknown\",2,2,\n"
-                    "            LENGTHUNIT[\"metre\",1]]],\n"
-                    "    PRIMEM[\"Reference meridian\",0,\n"
-                    "        ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
-                    "    CS[ellipsoidal,2],\n"
-                    "        AXIS[\"longitude\",east,\n"
-                    "            ORDER[1],\n"
-                    "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
-                    "        AXIS[\"latitude\",north,\n"
-                    "            ORDER[2],\n"
-                    "            ANGLEUNIT[\"degree\",0.0174532925199433]]]";
-    EXPECT_EQ(f->toString(), expected);
+    EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+    auto rf = crs->ellipsoid()->computedInverseFlattening();
+    EXPECT_EQ(rf, 2) << rf;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_longlat_a_e) {
+    auto obj =
+        PROJStringParser().createFromPROJString("+proj=longlat +a=2 +e=0.5");
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+    auto rf = crs->ellipsoid()->computedInverseFlattening();
+    EXPECT_NEAR(rf, 7.46410161513775, 1e-14) << rf;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_longlat_a_es) {
+    auto obj =
+        PROJStringParser().createFromPROJString("+proj=longlat +a=2 +es=0.5");
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+    auto rf = crs->ellipsoid()->computedInverseFlattening();
+    EXPECT_EQ(rf, 3.4142135623730958) << rf;
 }
 
 // ---------------------------------------------------------------------------
@@ -5711,23 +5721,31 @@ TEST(io, projparse_longlat_R) {
     auto obj = PROJStringParser().createFromPROJString("+proj=longlat +R=2");
     auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
     ASSERT_TRUE(crs != nullptr);
-    WKTFormatterNNPtr f(WKTFormatter::create());
-    f->simulCurNodeHasId();
-    crs->exportToWKT(f.get());
-    auto expected = "GEODCRS[\"unknown\",\n"
-                    "    DATUM[\"unknown\",\n"
-                    "        ELLIPSOID[\"unknown\",2,0,\n"
-                    "            LENGTHUNIT[\"metre\",1]]],\n"
-                    "    PRIMEM[\"Reference meridian\",0,\n"
-                    "        ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
-                    "    CS[ellipsoidal,2],\n"
-                    "        AXIS[\"longitude\",east,\n"
-                    "            ORDER[1],\n"
-                    "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
-                    "        AXIS[\"latitude\",north,\n"
-                    "            ORDER[2],\n"
-                    "            ANGLEUNIT[\"degree\",0.0174532925199433]]]";
-    EXPECT_EQ(f->toString(), expected);
+    EXPECT_TRUE(crs->ellipsoid()->isSphere());
+    EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_longlat_a) {
+    auto obj = PROJStringParser().createFromPROJString("+proj=longlat +a=2");
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_TRUE(crs->ellipsoid()->isSphere());
+    EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_longlat_a_override_ellps) {
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=longlat +a=2 +ellps=WGS84");
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_TRUE(!crs->ellipsoid()->isSphere());
+    EXPECT_EQ(crs->ellipsoid()->semiMajorAxis().getSIValue(), 2);
+    EXPECT_EQ(crs->ellipsoid()->computedInverseFlattening(), 298.25722356300003)
+        << crs->ellipsoid()->computedInverseFlattening();
 }
 
 // ---------------------------------------------------------------------------
@@ -7505,6 +7523,16 @@ TEST(io, projparse_init) {
         EXPECT_EQ(co->exportToPROJString(PROJStringFormatter::create().get()),
                   "+proj=pipeline +step +init=epsg:4326 +step +proj=longlat");
     }
+
+    {
+        auto obj = PROJStringParser().createFromPROJString(
+            "init=epsg:4326 proj=longlat ellps=GRS80");
+        auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
+                  "+proj=pipeline +step +proj=longlat +ellps=GRS80 "
+                  "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -7536,10 +7564,6 @@ TEST(io, projparse_errors) {
     EXPECT_THROW(PROJStringParser().createFromPROJString(
                      "proj=pipeline step init=epsg:4326 init=epsg:4326"),
                  ParsingException);
-
-    EXPECT_THROW(PROJStringParser().createFromPROJString(
-                     "proj=pipeline step init=epsg:4326 proj=longlat"),
-                 ParsingException);
 }
 
 // ---------------------------------------------------------------------------
@@ -7577,9 +7601,6 @@ TEST(io, projparse_longlat_errors) {
         PROJStringParser().createFromPROJString("+proj=longlat +R=invalid"),
         ParsingException);
 
-    EXPECT_THROW(PROJStringParser().createFromPROJString("+proj=longlat +a=1"),
-                 ParsingException);
-
     EXPECT_THROW(PROJStringParser().createFromPROJString("+proj=longlat +b=1"),
                  ParsingException);
 
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 528aa6b8..818c10ec 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4477,6 +4477,38 @@ TEST(operation, geogCRS_to_geogCRS_CH1903_to_CH1903plus_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, geogCRS_to_geogCRS_init_IGNF_to_init_IGNF_context) {
+
+    auto dbContext = DatabaseContext::create();
+
+    auto sourceCRS_obj = PROJStringParser()
+                             .attachDatabaseContext(dbContext)
+                             .setUsePROJ4InitRules(true)
+                             .createFromPROJString("+init=IGNF:NTFG");
+    auto sourceCRS = nn_dynamic_pointer_cast<CRS>(sourceCRS_obj);
+    ASSERT_TRUE(sourceCRS != nullptr);
+
+    auto targetCRS_obj = PROJStringParser()
+                             .attachDatabaseContext(dbContext)
+                             .setUsePROJ4InitRules(true)
+                             .createFromPROJString("+init=IGNF:RGF93G");
+    auto targetCRS = nn_dynamic_pointer_cast<CRS>(targetCRS_obj);
+    ASSERT_TRUE(targetCRS != nullptr);
+
+    auto authFactory = AuthorityFactory::create(dbContext, std::string());
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        NN_CHECK_ASSERT(sourceCRS), NN_CHECK_ASSERT(targetCRS), ctxt);
+    ASSERT_EQ(list.size(), 1);
+
+    EXPECT_EQ(list[0]->nameStr(),
+              "NOUVELLE TRIANGULATION DE LA FRANCE (NTF) vers RGF93 (ETRS89)");
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=hgridshift +grids=ntf_r93.gsb");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, geogCRS_to_geogCRS_3D) {
 
     auto geogcrs_m_obj =
@@ -5110,6 +5142,83 @@ TEST(operation, boundCRS_of_geogCRS_to_unrelated_geogCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, createOperation_boundCRS_identified_by_datum) {
+    auto objSrc =
+        PROJStringParser().createFromPROJString("+proj=longlat +datum=WGS84");
+    auto src = nn_dynamic_pointer_cast<GeographicCRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDest = PROJStringParser().createFromPROJString(
+        "+proj=utm +zone=32 +a=6378249.2 +b=6356515 "
+        "+towgs84=-263.0,6.0,431.0 +no_defs");
+    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 +proj=cart +ellps=WGS84 "
+              "+step +proj=helmert +x=263 +y=-6 +z=-431 "
+              "+step +inv +proj=cart +ellps=clrk80ign "
+              "+step +proj=utm +zone=32 +ellps=clrk80ign");
+
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), std::string());
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest), ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    EXPECT_TRUE(list[0]->isEquivalentTo(op.get()));
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, boundCRS_of_clrk_66_geogCRS_to_nad83_geogCRS) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=latlong +ellps=clrk66 +nadgrids=ntv1_can.dat,conus");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDest =
+        PROJStringParser().createFromPROJString("+proj=latlong +datum=NAD83");
+    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 +proj=hgridshift +grids=ntv1_can.dat,conus "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, boundCRS_of_clrk_66_projCRS_to_nad83_geogCRS) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=utm +zone=17 +ellps=clrk66 +nadgrids=ntv1_can.dat,conus");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDest =
+        PROJStringParser().createFromPROJString("+proj=latlong +datum=NAD83");
+    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 +inv +proj=utm +zone=17 +ellps=clrk66 "
+              "+step +proj=hgridshift +grids=ntv1_can.dat,conus "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, boundCRS_of_projCRS_to_geogCRS) {
     auto utm31 = ProjectedCRS::create(
         PropertyMap(), GeographicCRS::EPSG_4807,
@@ -5864,6 +5973,25 @@ TEST(operation, createOperation_on_crs_with_canonical_bound_crs) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, createOperation_fallback_to_proj4_strings) {
+    auto objDest = PROJStringParser().createFromPROJString(
+        "+proj=longlat +geoc +ellps=WGS84");
+    auto dest = nn_dynamic_pointer_cast<GeographicCRS>(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 +inv +proj=longlat +datum=WGS84 "
+              "+step +proj=longlat +geoc +ellps=WGS84 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, mercator_variant_A_to_variant_B) {
     auto projCRS = ProjectedCRS::create(
         PropertyMap(), GeographicCRS::EPSG_4326,