Add proj_create_derived_geographic_crs() and proj_create_conversion_pole_rotation_grib_convention() to address GRIB datasets using a pole rotation method

1 files changed, 91 insertions, 0 deletions

diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index d8816acf..f87f6589 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -2411,6 +2411,13 @@ TEST_F(CApi, proj_create_projections) {
         ObjectKeeper keeper_projCRS(projCRS);
         ASSERT_NE(projCRS, nullptr);
     }
+    {
+        auto projCRS = proj_create_conversion_pole_rotation_grib_convention(
+            m_ctxt, 0, 0, 0, "Degree", 0.0174532925199433);
+        ObjectKeeper keeper_projCRS(projCRS);
+        ASSERT_NE(projCRS, nullptr);
+    }
+
     /* END: Generated by scripts/create_c_api_projections.py*/
 }
 
@@ -4343,4 +4350,88 @@ TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) {
     EXPECT_EQ(std::string(proj_5_bis), std::string(proj_5));
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_create_derived_geographic_crs) {
+
+    PJ *crs_4326 = proj_create(m_ctxt, "EPSG:4326");
+    ObjectKeeper keeper_crs_4326(crs_4326);
+    ASSERT_NE(crs_4326, nullptr);
+
+    PJ *conversion = proj_create_conversion_pole_rotation_grib_convention(
+        m_ctxt, 2, 3, 4, "Degree", 0.0174532925199433);
+    ObjectKeeper keeper_conversion(conversion);
+    ASSERT_NE(conversion, nullptr);
+
+    PJ *cs = proj_crs_get_coordinate_system(m_ctxt, crs_4326);
+    ObjectKeeper keeper_cs(cs);
+    ASSERT_NE(cs, nullptr);
+
+    ASSERT_EQ(
+        proj_create_derived_geographic_crs(m_ctxt, "my rotated CRS",
+                                           conversion, // wrong type of object
+                                           conversion, cs),
+        nullptr);
+
+    ASSERT_EQ(
+        proj_create_derived_geographic_crs(m_ctxt, "my rotated CRS", crs_4326,
+                                           crs_4326, // wrong type of object
+                                           cs),
+        nullptr);
+
+    ASSERT_EQ(proj_create_derived_geographic_crs(
+                  m_ctxt, "my rotated CRS", crs_4326, conversion,
+                  conversion // wrong type of object
+                  ),
+              nullptr);
+
+    PJ *derived_crs = proj_create_derived_geographic_crs(
+        m_ctxt, "my rotated CRS", crs_4326, conversion, cs);
+    ObjectKeeper keeper_derived_crs(derived_crs);
+    ASSERT_NE(derived_crs, nullptr);
+
+    EXPECT_FALSE(proj_is_derived_crs(m_ctxt, crs_4326));
+    EXPECT_TRUE(proj_is_derived_crs(m_ctxt, derived_crs));
+
+    auto wkt = proj_as_wkt(m_ctxt, derived_crs, PJ_WKT2_2019, nullptr);
+    const char *expected_wkt =
+        "GEOGCRS[\"my rotated CRS\",\n"
+        "    BASEGEOGCRS[\"WGS 84\",\n"
+        "        DATUM[\"World Geodetic System 1984\",\n"
+        "            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                LENGTHUNIT[\"metre\",1]]],\n"
+        "        PRIMEM[\"Greenwich\",0,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
+        "    DERIVINGCONVERSION[\"Pole rotation (GRIB convention)\",\n"
+        "        METHOD[\"Pole rotation (GRIB convention)\"],\n"
+        "        PARAMETER[\"Latitude of the southern pole (GRIB "
+        "convention)\",2,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433,\n"
+        "                ID[\"EPSG\",9122]]],\n"
+        "        PARAMETER[\"Longitude of the southern pole (GRIB "
+        "convention)\",3,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433,\n"
+        "                ID[\"EPSG\",9122]]],\n"
+        "        PARAMETER[\"Axis rotation (GRIB convention)\",4,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433,\n"
+        "                ID[\"EPSG\",9122]]]],\n"
+        "    CS[ellipsoidal,2],\n"
+        "        AXIS[\"geodetic latitude (Lat)\",north,\n"
+        "            ORDER[1],\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433,\n"
+        "                ID[\"EPSG\",9122]]],\n"
+        "        AXIS[\"geodetic longitude (Lon)\",east,\n"
+        "            ORDER[2],\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433,\n"
+        "                ID[\"EPSG\",9122]]]]";
+
+    ASSERT_NE(wkt, nullptr);
+    EXPECT_EQ(wkt, std::string(expected_wkt));
+
+    auto proj_5 = proj_as_proj_string(m_ctxt, derived_crs, PJ_PROJ_5, nullptr);
+    ASSERT_NE(proj_5, nullptr);
+    EXPECT_EQ(proj_5, std::string("+proj=ob_tran +o_proj=longlat +o_lon_p=-4 "
+                                  "+o_lat_p=-2 +lon_0=3 +datum=WGS84 +no_defs "
+                                  "+type=crs"));
+}
 } // namespace