Assorted set of fixes for PROJString to ISO19111 model:

- createFromPROJString(): take into account axisswap step for Krovak and Transverse Mercator (South Orientated)
- Geocentric export to PROJ4: use datum when possible, and add explicit units=m
- ESRI WKT parser: make it case insensitive to parameter and projection names, and more tolerant about possible parameter name aliases
- import from WKT1 for Polar_Stereographic: don't be case sensitive
- importFromPROJString: allow pm to override datum
- Equidistant cylindrical: add support for non-standard latitude of natural origin, used in a GDAL test case
- tmerc export to PROJString: use 'k' instead of 'k_0'
- pj_ellps: use official value from EPSG for reverse flattening of Airy ellipsoid
- GDAL compatibility: add support for importing odd formulations of Mercator as WKT1, but rejecting them when exporting to PROJ
- Add export of 'Geostationary Satellite (Sweep X)' to WKT1_GDAL via EXTENSION.PROJ4 node
- importFromPROJString: add support for +f
- WKT1 / PROJ4: add support for EXTENSION.PROJ4 nodes and +wktext
- exportToWKT: change way we deal with AXIS by default for WKT1_GDAL
- Improve etmerc handling
- Fix WKT import of peg_point_heading for Spherical_Cross_Track_Height
- International Map of the World Polyconic: change parameter mapping
- exportToPROJ: add alpha parameter
- Hotine_Oblique_Mercator_Two_Point_Natural_Origin: GDAL_WKT1 related fix
- GDAL compatibility improvements in import from PROJ4 / WKT1 for polar stereographic
- Add support for +towgs84 when importing a +proj=geocent
- import from WKT1: add support for an odd Mercator_1SP formulation handled by GDAL
- export to proj4 strings: add +units=m to projected CRS for better GDAL compatibility
- export to proj4 strings: add +no_defs to CRS for better GDAL compatibility

4 files changed, 852 insertions, 140 deletions

diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 207a8cd1..8c9f114b 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -263,6 +263,22 @@ TEST_F(CApi, proj_obj_as_wkt) {
         EXPECT_TRUE(std::string(wkt).find("AXIS") == std::string::npos) << wkt;
     }
 
+    // OUTPUT_AXIS=AUTO
+    {
+        const char *const options[] = {"OUTPUT_AXIS=AUTO", nullptr};
+        auto wkt = proj_obj_as_wkt(obj, PJ_WKT1_GDAL, options);
+        ASSERT_NE(wkt, nullptr);
+        EXPECT_TRUE(std::string(wkt).find("AXIS") == std::string::npos) << wkt;
+    }
+
+    // OUTPUT_AXIS=YES
+    {
+        const char *const options[] = {"OUTPUT_AXIS=YES", nullptr};
+        auto wkt = proj_obj_as_wkt(obj, PJ_WKT1_GDAL, options);
+        ASSERT_NE(wkt, nullptr);
+        EXPECT_TRUE(std::string(wkt).find("AXIS") != std::string::npos) << wkt;
+    }
+
     // unsupported option
     {
         const char *const options[] = {"unsupported=yes", nullptr};
@@ -307,7 +323,7 @@ TEST_F(CApi, proj_obj_as_proj_string) {
     {
         auto proj_4 = proj_obj_as_proj_string(obj, PJ_PROJ_4, nullptr);
         ASSERT_NE(proj_4, nullptr);
-        EXPECT_EQ(std::string(proj_4), "+proj=longlat +datum=WGS84");
+        EXPECT_EQ(std::string(proj_4), "+proj=longlat +datum=WGS84 +no_defs");
     }
 }
 
@@ -327,7 +343,7 @@ TEST_F(CApi, proj_obj_as_proj_string_incompatible_WKT1) {
 
 // ---------------------------------------------------------------------------
 
-TEST_F(CApi, proj_obj_as_proj_string_etmerc_option) {
+TEST_F(CApi, proj_obj_as_proj_string_etmerc_option_yes) {
     auto obj = proj_obj_create_from_proj_string(m_ctxt, "+proj=tmerc", nullptr);
     ObjectKeeper keeper(obj);
     ASSERT_NE(obj, nullptr);
@@ -335,8 +351,24 @@ TEST_F(CApi, proj_obj_as_proj_string_etmerc_option) {
     const char *options[] = {"USE_ETMERC=YES", nullptr};
     auto str = proj_obj_as_proj_string(obj, PJ_PROJ_4, options);
     ASSERT_NE(str, nullptr);
-    EXPECT_EQ(str, std::string("+proj=etmerc +lat_0=0 +lon_0=0 +k_0=1 +x_0=0 "
-                               "+y_0=0 +datum=WGS84"));
+    EXPECT_EQ(str, std::string("+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 "
+                               "+y_0=0 +datum=WGS84 +units=m +no_defs"));
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_obj_as_proj_string_etmerc_option_no) {
+    auto obj =
+        proj_obj_create_from_proj_string(m_ctxt, "+proj=utm +zone=31", nullptr);
+    ObjectKeeper keeper(obj);
+    ASSERT_NE(obj, nullptr);
+
+    const char *options[] = {"USE_ETMERC=NO", nullptr};
+    auto str = proj_obj_as_proj_string(obj, PJ_PROJ_4, options);
+    ASSERT_NE(str, nullptr);
+    EXPECT_EQ(str, std::string("+proj=tmerc +lat_0=0 +lon_0=3 +k=0.9996 "
+                               "+x_0=500000 +y_0=0 +datum=WGS84 +units=m "
+                               "+no_defs"));
 }
 
 // ---------------------------------------------------------------------------
@@ -356,7 +388,8 @@ TEST_F(CApi, proj_obj_crs_create_bound_crs_to_WGS84) {
     EXPECT_EQ(std::string(proj_4),
               "+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 "
               "+y_0=500000 +ellps=krass "
-              "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69");
+              "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69 "
+              "+units=m +no_defs");
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 43698252..84be1b22 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -302,22 +302,39 @@ TEST(crs, EPSG_4326_as_WKT2_2018_SIMPLIFIED) {
 
 TEST(crs, EPSG_4326_as_WKT1_GDAL) {
     auto crs = GeographicCRS::EPSG_4326;
-    WKTFormatterNNPtr f(
-        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL));
-    crs->exportToWKT(f.get());
-    EXPECT_EQ(f->toString(),
-              "GEOGCS[\"WGS 84\",\n"
-              "    DATUM[\"WGS_1984\",\n"
-              "        SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
-              "            AUTHORITY[\"EPSG\",\"7030\"]],\n"
-              "        AUTHORITY[\"EPSG\",\"6326\"]],\n"
-              "    PRIMEM[\"Greenwich\",0,\n"
-              "        AUTHORITY[\"EPSG\",\"8901\"]],\n"
-              "    UNIT[\"degree\",0.0174532925199433,\n"
-              "        AUTHORITY[\"EPSG\",\"9122\"]],\n"
-              "    AXIS[\"Latitude\",NORTH],\n"
-              "    AXIS[\"Longitude\",EAST],\n"
-              "    AUTHORITY[\"EPSG\",\"4326\"]]");
+    auto wkt = crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+    EXPECT_EQ(wkt, "GEOGCS[\"WGS 84\",\n"
+                   "    DATUM[\"WGS_1984\",\n"
+                   "        SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+                   "            AUTHORITY[\"EPSG\",\"7030\"]],\n"
+                   "        AUTHORITY[\"EPSG\",\"6326\"]],\n"
+                   "    PRIMEM[\"Greenwich\",0,\n"
+                   "        AUTHORITY[\"EPSG\",\"8901\"]],\n"
+                   "    UNIT[\"degree\",0.0174532925199433,\n"
+                   "        AUTHORITY[\"EPSG\",\"9122\"]],\n"
+                   "    AUTHORITY[\"EPSG\",\"4326\"]]");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(crs, EPSG_4326_as_WKT1_GDAL_with_axis) {
+    auto crs = GeographicCRS::EPSG_4326;
+    auto wkt = crs->exportToWKT(
+        &(WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL)
+              ->setOutputAxis(WKTFormatter::OutputAxisRule::YES)));
+    EXPECT_EQ(wkt, "GEOGCS[\"WGS 84\",\n"
+                   "    DATUM[\"WGS_1984\",\n"
+                   "        SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+                   "            AUTHORITY[\"EPSG\",\"7030\"]],\n"
+                   "        AUTHORITY[\"EPSG\",\"6326\"]],\n"
+                   "    PRIMEM[\"Greenwich\",0,\n"
+                   "        AUTHORITY[\"EPSG\",\"8901\"]],\n"
+                   "    UNIT[\"degree\",0.0174532925199433,\n"
+                   "        AUTHORITY[\"EPSG\",\"9122\"]],\n"
+                   "    AXIS[\"Latitude\",NORTH],\n"
+                   "    AXIS[\"Longitude\",EAST],\n"
+                   "    AUTHORITY[\"EPSG\",\"4326\"]]");
 }
 
 // ---------------------------------------------------------------------------
@@ -356,7 +373,7 @@ TEST(crs, EPSG_4326_as_PROJ_string) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +datum=WGS84");
+        "+proj=longlat +datum=WGS84 +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -403,28 +420,45 @@ TEST(crs, EPSG_4979_as_WKT2_2018_SIMPLIFIED) {
 
 // ---------------------------------------------------------------------------
 
-TEST(crs, EPSG_4979_as_WKT1_GDAL) {
+TEST(crs, EPSG_4979_as_WKT1_GDAL_with_axis) {
     auto crs = GeographicCRS::EPSG_4979;
-    WKTFormatterNNPtr f(
-        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL));
-    crs->exportToWKT(f.get());
-    // FIXME? WKT1 only supports 2 axis for GEOGCS. So this is an extension of
+    auto wkt = crs->exportToWKT(
+        &(WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL)
+              ->setOutputAxis(WKTFormatter::OutputAxisRule::YES)));
+    // WKT1 only supports 2 axis for GEOGCS. So this is an extension of
     // WKT1 as it
     // and GDAL doesn't really export such as beast, although it can import it
-    EXPECT_EQ(f->toString(),
-              "GEOGCS[\"WGS 84\",\n"
-              "    DATUM[\"WGS_1984\",\n"
-              "        SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
-              "            AUTHORITY[\"EPSG\",\"7030\"]],\n"
-              "        AUTHORITY[\"EPSG\",\"6326\"]],\n"
-              "    PRIMEM[\"Greenwich\",0,\n"
-              "        AUTHORITY[\"EPSG\",\"8901\"]],\n"
-              "    UNIT[\"degree\",0.0174532925199433,\n"
-              "        AUTHORITY[\"EPSG\",\"9122\"]],\n"
-              "    AXIS[\"Latitude\",NORTH],\n"
-              "    AXIS[\"Longitude\",EAST],\n"
-              "    AXIS[\"Ellipsoidal height\",UP],\n"
-              "    AUTHORITY[\"EPSG\",\"4979\"]]");
+    EXPECT_EQ(wkt, "GEOGCS[\"WGS 84\",\n"
+                   "    DATUM[\"WGS_1984\",\n"
+                   "        SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+                   "            AUTHORITY[\"EPSG\",\"7030\"]],\n"
+                   "        AUTHORITY[\"EPSG\",\"6326\"]],\n"
+                   "    PRIMEM[\"Greenwich\",0,\n"
+                   "        AUTHORITY[\"EPSG\",\"8901\"]],\n"
+                   "    UNIT[\"degree\",0.0174532925199433,\n"
+                   "        AUTHORITY[\"EPSG\",\"9122\"]],\n"
+                   "    AXIS[\"Latitude\",NORTH],\n"
+                   "    AXIS[\"Longitude\",EAST],\n"
+                   "    AXIS[\"Ellipsoidal height\",UP],\n"
+                   "    AUTHORITY[\"EPSG\",\"4979\"]]");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(crs, EPSG_4979_as_WKT1_GDAL) {
+    auto crs = GeographicCRS::EPSG_4979;
+    auto wkt = crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+    EXPECT_EQ(wkt, "GEOGCS[\"WGS 84\",\n"
+                   "    DATUM[\"WGS_1984\",\n"
+                   "        SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+                   "            AUTHORITY[\"EPSG\",\"7030\"]],\n"
+                   "        AUTHORITY[\"EPSG\",\"6326\"]],\n"
+                   "    PRIMEM[\"Greenwich\",0,\n"
+                   "        AUTHORITY[\"EPSG\",\"8901\"]],\n"
+                   "    UNIT[\"degree\",0.0174532925199433,\n"
+                   "        AUTHORITY[\"EPSG\",\"9122\"]],\n"
+                   "    AUTHORITY[\"EPSG\",\"4979\"]]");
 }
 
 // ---------------------------------------------------------------------------
@@ -512,8 +546,6 @@ TEST(crs, EPSG_4807_as_WKT1_GDAL) {
         "        AUTHORITY[\"EPSG\",\"8903\"]],\n"
         "    UNIT[\"grad\",0.015707963267949,\n"
         "        AUTHORITY[\"EPSG\",\"9105\"]],\n"
-        "    AXIS[\"Latitude\",NORTH],\n"
-        "    AXIS[\"Longitude\",EAST],\n"
         "    AUTHORITY[\"EPSG\",\"4807\"]]");
 }
 
@@ -554,7 +586,7 @@ TEST(crs, EPSG_4807_as_PROJ_string) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +ellps=clrk80ign +pm=paris");
+        "+proj=longlat +ellps=clrk80ign +pm=paris +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -580,7 +612,7 @@ TEST(crs, EPSG_4267) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +datum=NAD27");
+        "+proj=longlat +datum=NAD27 +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -619,7 +651,7 @@ TEST(crs, EPSG_4269) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +datum=NAD83");
+        "+proj=longlat +datum=NAD83 +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -687,7 +719,7 @@ TEST(crs, EPSG_27561_projected_with_geodetic_in_grad_as_PROJ_string_and_WKT1) {
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
         "+proj=lcc +lat_1=49.5 +lat_0=49.5 +lon_0=0 +k_0=0.999877341 "
-        "+x_0=600000 +y_0=200000 +ellps=clrk80ign +pm=paris");
+        "+x_0=600000 +y_0=200000 +ellps=clrk80ign +pm=paris +units=m +no_defs");
 
     auto nn_crs = NN_CHECK_ASSERT(crs);
     EXPECT_TRUE(nn_crs->isEquivalentTo(nn_crs.get()));
@@ -704,9 +736,7 @@ TEST(crs, EPSG_27561_projected_with_geodetic_in_grad_as_PROJ_string_and_WKT1) {
         "        DATUM[\"Nouvelle_Triangulation_Francaise_Paris\",\n"
         "            SPHEROID[\"Clarke 1880 (IGN)\",6378249.2,293.4660213]],\n"
         "        PRIMEM[\"Paris\",2.33722917000759],\n"
-        "        UNIT[\"grad\",0.015707963268],\n"
-        "        AXIS[\"Latitude\",NORTH],\n"
-        "        AXIS[\"Longitude\",EAST]],\n"
+        "        UNIT[\"grad\",0.015707963268]],\n"
         "    PROJECTION[\"Lambert_Conformal_Conic_1SP\"],\n"
         "    PARAMETER[\"latitude_of_origin\",55],\n"
         "    PARAMETER[\"central_meridian\",0],\n"
@@ -799,15 +829,15 @@ TEST(crs, EPSG_2222_projected_unit_foot_as_PROJ_string_and_WKT1) {
     EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
               "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=tmerc "
-              "+lat_0=31 +lon_0=-110.166666666667 +k_0=0.9999 +x_0=213360 "
+              "+lat_0=31 +lon_0=-110.166666666667 +k=0.9999 +x_0=213360 "
               "+y_0=0 +ellps=GRS80 +step +proj=unitconvert +xy_in=m +z_in=m "
               "+xy_out=ft +z_out=ft");
     EXPECT_EQ(
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=tmerc +lat_0=31 +lon_0=-110.166666666667 +k_0=0.9999 "
-        "+x_0=213360 +y_0=0 +datum=NAD83 +units=ft");
+        "+proj=tmerc +lat_0=31 +lon_0=-110.166666666667 +k=0.9999 "
+        "+x_0=213360 +y_0=0 +datum=NAD83 +units=ft +no_defs");
 
     auto wkt1 = crs->exportToWKT(
         WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
@@ -819,9 +849,7 @@ TEST(crs, EPSG_2222_projected_unit_foot_as_PROJ_string_and_WKT1) {
               "        PRIMEM[\"Greenwich\",0,\n"
               "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
               "        UNIT[\"degree\",0.0174532925199433,\n"
-              "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
-              "        AXIS[\"Latitude\",NORTH],\n"
-              "        AXIS[\"Longitude\",EAST]],\n"
+              "            AUTHORITY[\"EPSG\",\"9122\"]]],\n"
               "    PROJECTION[\"Transverse_Mercator\"],\n"
               "    PARAMETER[\"latitude_of_origin\",31],\n"
               "    PARAMETER[\"central_meridian\",-110.166666666667],\n"
@@ -867,9 +895,7 @@ TEST(crs, projected_with_parameter_unit_different_than_cs_unit_as_WKT1) {
               "            SPHEROID[\"GRS 1980\",6378137,298.257222101]],\n"
               "        PRIMEM[\"Greenwich\",0],\n"
               "        UNIT[\"degree\",0.0174532925199433,\n"
-              "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
-              "        AXIS[\"Latitude\",NORTH],\n"
-              "        AXIS[\"Longitude\",EAST]],\n"
+              "            AUTHORITY[\"EPSG\",\"9122\"]]],\n"
               "    PROJECTION[\"Transverse_Mercator\"],\n"
               "    PARAMETER[\"latitude_of_origin\",0],\n"
               "    PARAMETER[\"central_meridian\",9],\n"
@@ -1081,9 +1107,6 @@ TEST(crs, geocentricCRS_as_WKT1_GDAL) {
               "        AUTHORITY[\"EPSG\",\"8901\"]],\n"
               "    UNIT[\"metre\",1,\n"
               "        AUTHORITY[\"EPSG\",\"9001\"]],\n"
-              "    AXIS[\"Geocentric X\",OTHER],\n"
-              "    AXIS[\"Geocentric Y\",OTHER],\n"
-              "    AXIS[\"Geocentric Z\",NORTH],\n"
               "    AUTHORITY[\"EPSG\",\"4328\"]]");
 }
 
@@ -1097,7 +1120,7 @@ TEST(crs, geocentricCRS_as_PROJ_string) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=geocent +ellps=WGS84");
+        "+proj=geocent +datum=WGS84 +units=m +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -1597,8 +1620,6 @@ TEST(crs, projectedCRS_as_WKT1_GDAL) {
                     "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
                     "        UNIT[\"degree\",0.0174532925199433,\n"
                     "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
-                    "        AXIS[\"Latitude\",NORTH],\n"
-                    "        AXIS[\"Longitude\",EAST],\n"
                     "        AUTHORITY[\"EPSG\",\"4326\"]],\n"
                     "    PROJECTION[\"Transverse_Mercator\"],\n"
                     "    PARAMETER[\"latitude_of_origin\",0],\n"
@@ -1654,7 +1675,7 @@ TEST(crs, projectedCRS_as_PROJ_string) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=utm +zone=31 +datum=WGS84");
+        "+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -2870,8 +2891,6 @@ TEST(crs, compoundCRS_as_WKT1_GDAL) {
         "                AUTHORITY[\"EPSG\",\"8901\"]],\n"
         "            UNIT[\"degree\",0.0174532925199433,\n"
         "                AUTHORITY[\"EPSG\",\"9122\"]],\n"
-        "            AXIS[\"Latitude\",NORTH],\n"
-        "            AXIS[\"Longitude\",EAST],\n"
         "            AUTHORITY[\"EPSG\",\"4326\"]],\n"
         "        PROJECTION[\"Transverse_Mercator\"],\n"
         "        PARAMETER[\"latitude_of_origin\",0],\n"
@@ -2913,7 +2932,7 @@ TEST(crs, compoundCRS_as_PROJ_string) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=utm +zone=31 +datum=WGS84 +vunits=m");
+        "+proj=utm +zone=31 +datum=WGS84 +units=m +vunits=m +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -3301,9 +3320,7 @@ TEST(crs, boundCRS_to_WKT1) {
                     "        PRIMEM[\"Greenwich\",0,\n"
                     "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
                     "        UNIT[\"degree\",0.0174532925199433,\n"
-                    "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
-                    "        AXIS[\"Latitude\",NORTH],\n"
-                    "        AXIS[\"Longitude\",EAST]],\n"
+                    "            AUTHORITY[\"EPSG\",\"9122\"]]],\n"
                     "    PROJECTION[\"Transverse_Mercator\"],\n"
                     "    PARAMETER[\"latitude_of_origin\",0],\n"
                     "    PARAMETER[\"central_meridian\",3],\n"
@@ -3339,7 +3356,7 @@ TEST(crs, boundCRS_geographicCRS_to_PROJ_string) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +ellps=WGS84 +towgs84=1,2,3,4,5,6,7");
+        "+proj=longlat +ellps=WGS84 +towgs84=1,2,3,4,5,6,7 +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -3362,7 +3379,8 @@ TEST(crs, boundCRS_projectedCRS_to_PROJ_string) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=utm +zone=31 +ellps=WGS84 +towgs84=1,2,3,4,5,6,7");
+        "+proj=utm +zone=31 +ellps=WGS84 +towgs84=1,2,3,4,5,6,7 +units=m "
+        "+no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -3435,9 +3453,7 @@ TEST(crs, WKT1_DATUM_EXTENSION_to_WKT1_and_PROJ_string) {
         "            SPHEROID[\"intl\",6378388,297],\n"
         "            EXTENSION[\"PROJ4_GRIDS\",\"nzgd2kgrid0005.gsb\"]],\n"
         "        PRIMEM[\"Greenwich\",0],\n"
-        "        UNIT[\"degree\",0.0174532925199433],\n"
-        "        AXIS[\"Longitude\",EAST],\n"
-        "        AXIS[\"Latitude\",NORTH]],\n"
+        "        UNIT[\"degree\",0.0174532925199433]],\n"
         "    PROJECTION[\"New_Zealand_Map_Grid\"],\n"
         "    PARAMETER[\"latitude_of_origin\",-41],\n"
         "    PARAMETER[\"central_meridian\",173],\n"
@@ -3461,7 +3477,7 @@ TEST(crs, WKT1_DATUM_EXTENSION_to_WKT1_and_PROJ_string) {
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
         "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 "
-        "+ellps=intl +nadgrids=nzgd2kgrid0005.gsb +units=m");
+        "+ellps=intl +nadgrids=nzgd2kgrid0005.gsb +units=m +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -3473,7 +3489,7 @@ TEST(crs, WKT1_VERT_DATUM_EXTENSION_to_WKT1) {
                "        AUTHORITY[\"EPSG\",\"1027\"]],\n"
                "    UNIT[\"metre\",1,\n"
                "        AUTHORITY[\"EPSG\",\"9001\"]],\n"
-               "    AXIS[\"Up\",UP],\n"
+               "    AXIS[\"Gravity-related height\",UP],\n"
                "    AUTHORITY[\"EPSG\",\"3855\"]]";
 
     auto obj = WKTParser().createFromWKT(wkt);
@@ -3561,7 +3577,7 @@ TEST(crs, WKT1_VERT_DATUM_EXTENSION_to_PROJ_string) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+geoidgrids=egm08_25.gtx +vunits=m");
+        "+geoidgrids=egm08_25.gtx +vunits=m +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -4197,9 +4213,7 @@ TEST(crs, engineeringCRS_WKT2) {
 TEST(crs, engineeringCRS_WKT1) {
 
     auto expected = "LOCAL_CS[\"Engineering CRS\",\n"
-                    "    LOCAL_DATUM[\"Engineering datum\",32767],\n"
-                    "    AXIS[\"Easting\",EAST],\n"
-                    "    AXIS[\"Northing\",NORTH]]";
+                    "    LOCAL_DATUM[\"Engineering datum\",32767]]";
     EXPECT_EQ(
         createEngineeringCRS()->exportToWKT(
             WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
@@ -4522,7 +4536,8 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
                           .get()),
                   "+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 "
                   "+y_0=500000 +ellps=krass "
-                  "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69");
+                  "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69 "
+                  "+units=m +no_defs");
     }
     {
         // Pulkovo 42 Poland
@@ -4538,7 +4553,8 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
                           .get()),
                   "+proj=sterea +lat_0=50.625 +lon_0=21.0833333333333 "
                   "+k=0.9998 +x_0=4637000 +y_0=5647000 +ellps=krass "
-                  "+towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84");
+                  "+towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 "
+                  "+units=m +no_defs");
     }
     {
         // NTF (Paris)
@@ -4553,7 +4569,7 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
                           PROJStringFormatter::Convention::PROJ_4)
                           .get()),
                   "+proj=longlat +ellps=clrk80ign +pm=paris "
-                  "+towgs84=-168,-60,320,0,0,0,0");
+                  "+towgs84=-168,-60,320,0,0,0,0 +no_defs");
     }
     {
         // NTF (Paris) / Lambert zone II + NGF-IGN69 height
@@ -4569,7 +4585,8 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
                           .get()),
                   "+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 "
                   "+x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris "
-                  "+towgs84=-168,-60,320,0,0,0,0 +vunits=m");
+                  "+towgs84=-168,-60,320,0,0,0,0 +units=m "
+                  "+vunits=m +no_defs");
     }
     {
         auto crs = createVerticalCRS();
@@ -4589,7 +4606,7 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
                           .get()),
                   "+proj=stere +lat_0=-90 +lat_ts=-67 +lon_0=140 +x_0=300000 "
                   "+y_0=-2299363.482 +ellps=intl "
-                  "+towgs84=324.912,153.282,172.026,0,0,0,0");
+                  "+towgs84=324.912,153.282,172.026,0,0,0,0 +units=m +no_defs");
     }
     {
         auto factoryIGNF =
@@ -4604,7 +4621,8 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
                           PROJStringFormatter::Convention::PROJ_4)
                           .get()),
                   "+proj=geocent +ellps=intl "
-                  "+towgs84=109.753,-528.133,-362.244,0,0,0,0");
+                  "+towgs84=109.753,-528.133,-362.244,0,0,0,0 +units=m "
+                  "+no_defs");
     }
 }
 
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 7d26d82a..26597f09 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -453,6 +453,69 @@ static std::string contentWKT2_EPSG_4326(
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt1_geographic_with_PROJ4_extension) {
+    auto wkt = "GEOGCS[\"WGS 84\",\n"
+               "    DATUM[\"unknown\",\n"
+               "        SPHEROID[\"WGS84\",6378137,298.257223563]],\n"
+               "    PRIMEM[\"Greenwich\",0],\n"
+               "    UNIT[\"degree\",0.0174532925199433],\n"
+               "    EXTENSION[\"PROJ4\",\"+proj=longlat +foo=bar +wktext\"]]";
+    auto obj = WKTParser().createFromWKT(wkt);
+
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
+        wkt);
+
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=longlat +foo=bar +wktext");
+
+    EXPECT_TRUE(
+        crs->exportToWKT(WKTFormatter::create().get()).find("EXTENSION") ==
+        std::string::npos);
+
+    EXPECT_TRUE(
+        crs->exportToWKT(
+               WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI).get())
+            .find("EXTENSION") == std::string::npos);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_geocentric_with_PROJ4_extension) {
+    auto wkt = "GEOCCS[\"WGS 84\",\n"
+               "    DATUM[\"unknown\",\n"
+               "        SPHEROID[\"WGS84\",6378137,298.257223563]],\n"
+               "    PRIMEM[\"Greenwich\",0],\n"
+               "    UNIT[\"Meter\",1],\n"
+               "    EXTENSION[\"PROJ4\",\"+proj=geocent +foo=bar +wktext\"]]";
+    auto obj = WKTParser().createFromWKT(wkt);
+
+    auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
+        wkt);
+
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=geocent +foo=bar +wktext");
+
+    EXPECT_TRUE(
+        crs->exportToWKT(WKTFormatter::create().get()).find("EXTENSION") ==
+        std::string::npos);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, wkt2_GEODCRS_EPSG_4326) {
     auto obj = WKTParser().createFromWKT("GEODCRS" + contentWKT2_EPSG_4326);
     auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
@@ -860,6 +923,49 @@ TEST(wkt_parse, wkt1_projected_no_axis) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt1_projected_with_PROJ4_extension) {
+    auto wkt = "PROJCS[\"unnamed\",\n"
+               "    GEOGCS[\"WGS 84\",\n"
+               "        DATUM[\"unknown\",\n"
+               "            SPHEROID[\"WGS84\",6378137,298.257223563]],\n"
+               "        PRIMEM[\"Greenwich\",0],\n"
+               "        UNIT[\"degree\",0.0174532925199433]],\n"
+               "    PROJECTION[\"Mercator_1SP\"],\n"
+               "    PARAMETER[\"central_meridian\",0],\n"
+               "    PARAMETER[\"scale_factor\",1],\n"
+               "    PARAMETER[\"false_easting\",0],\n"
+               "    PARAMETER[\"false_northing\",0],\n"
+               "    UNIT[\"Meter\",1],\n"
+               "    AXIS[\"Easting\",EAST],\n"
+               "    AXIS[\"Northing\",NORTH],\n"
+               "    EXTENSION[\"PROJ4\",\"+proj=merc +wktext\"]]";
+    auto obj = WKTParser().createFromWKT(wkt);
+
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
+        wkt);
+
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=merc +wktext");
+
+    EXPECT_TRUE(
+        crs->exportToWKT(WKTFormatter::create().get()).find("EXTENSION") ==
+        std::string::npos);
+
+    EXPECT_TRUE(
+        crs->exportToWKT(
+               WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI).get())
+            .find("EXTENSION") == std::string::npos);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, wkt1_krovak_south_west) {
     auto wkt =
         "PROJCS[\"S-JTSK / Krovak\","
@@ -942,13 +1048,15 @@ TEST(wkt_parse, wkt1_krovak_south_west) {
     auto expectedPROJString =
         "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
         "+step +proj=krovak +axis=swu +lat_0=49.5 "
-        "+lon_0=24.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel";
+        "+lon_0=24.8333333333333 +alpha=30.2881397222222 "
+        "+k=0.9999 +x_0=0 +y_0=0 +ellps=bessel";
     EXPECT_EQ(projString, expectedPROJString);
 
     obj = PROJStringParser().createFromPROJString(projString);
     crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
     ASSERT_TRUE(crs != nullptr);
     auto wkt2 = crs->exportToWKT(WKTFormatter::create().get());
+    EXPECT_TRUE(wkt2.find("METHOD[\"Krovak\"") != std::string::npos) << wkt2;
     EXPECT_TRUE(
         wkt2.find("PARAMETER[\"Latitude of pseudo standard parallel\",78.5,") !=
         std::string::npos)
@@ -959,6 +1067,17 @@ TEST(wkt_parse, wkt1_krovak_south_west) {
         << wkt2;
     EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
               expectedPROJString);
+
+    obj = PROJStringParser().createFromPROJString(
+        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step +proj=krovak +lat_0=49.5 "
+        "+lon_0=24.8333333333333 +alpha=30.2881397222222 "
+        "+k=0.9999 +x_0=0 +y_0=0 +ellps=bessel "
+        "+step +proj=axisswap +order=-2,-1");
+    crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    wkt2 = crs->exportToWKT(WKTFormatter::create().get());
+    EXPECT_TRUE(wkt2.find("METHOD[\"Krovak\"") != std::string::npos) << wkt2;
 }
 
 // ---------------------------------------------------------------------------
@@ -1043,7 +1162,174 @@ TEST(wkt_parse, wkt1_krovak_north_oriented) {
     EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
               "+step +proj=krovak +lat_0=49.5 +lon_0=24.8333333333333 "
-              "+k=0.9999 +x_0=0 +y_0=0 +ellps=bessel");
+              "+alpha=30.2881397222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_polar_stereographic_latitude_of_origin_70) {
+    auto wkt = "PROJCS[\"unknown\",\n"
+               "    GEOGCS[\"unknown\",\n"
+               "        DATUM[\"WGS_1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+               "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+               "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+               "        UNIT[\"degree\",0.0174532925199433,\n"
+               "            AUTHORITY[\"EPSG\",\"9122\"]]],\n"
+               "    PROJECTION[\"Polar_Stereographic\"],\n"
+               "    PARAMETER[\"latitude_of_origin\",70],\n"
+               "    PARAMETER[\"central_meridian\",2],\n"
+               "    PARAMETER[\"false_easting\",3],\n"
+               "    PARAMETER[\"false_northing\",4],\n"
+               "    UNIT[\"metre\",1,\n"
+               "        AUTHORITY[\"EPSG\",\"9001\"]]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto projString = crs->exportToPROJString(
+        PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+            .get());
+    auto expectedPROJString = "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=2 "
+                              "+x_0=3 +y_0=4 +datum=WGS84 +units=m +no_defs";
+    EXPECT_EQ(projString, expectedPROJString);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_polar_stereographic_latitude_of_origin_90) {
+    auto wkt = "PROJCS[\"unknown\",\n"
+               "    GEOGCS[\"unknown\",\n"
+               "        DATUM[\"WGS_1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+               "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+               "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+               "        UNIT[\"degree\",0.0174532925199433,\n"
+               "            AUTHORITY[\"EPSG\",\"9122\"]]],\n"
+               "    PROJECTION[\"Polar_Stereographic\"],\n"
+               "    PARAMETER[\"latitude_of_origin\",90],\n"
+               "    PARAMETER[\"central_meridian\",2],\n"
+               "    PARAMETER[\"false_easting\",3],\n"
+               "    PARAMETER[\"false_northing\",4],\n"
+               "    UNIT[\"metre\",1,\n"
+               "        AUTHORITY[\"EPSG\",\"9001\"]]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto projString = crs->exportToPROJString(
+        PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+            .get());
+    auto expectedPROJString = "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=2 "
+                              "+x_0=3 +y_0=4 +datum=WGS84 +units=m +no_defs";
+    EXPECT_EQ(projString, expectedPROJString);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_polar_stereographic_latitude_of_origin_90_scale_factor_1) {
+    auto wkt = "PROJCS[\"unknown\",\n"
+               "    GEOGCS[\"unknown\",\n"
+               "        DATUM[\"WGS_1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+               "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+               "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+               "        UNIT[\"degree\",0.0174532925199433,\n"
+               "            AUTHORITY[\"EPSG\",\"9122\"]]],\n"
+               "    PROJECTION[\"Polar_Stereographic\"],\n"
+               "    PARAMETER[\"latitude_of_origin\",90],\n"
+               "    PARAMETER[\"central_meridian\",2],\n"
+               "    PARAMETER[\"scale_factor\",1],\n"
+               "    PARAMETER[\"false_easting\",3],\n"
+               "    PARAMETER[\"false_northing\",4],\n"
+               "    UNIT[\"metre\",1,\n"
+               "        AUTHORITY[\"EPSG\",\"9001\"]]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto projString = crs->exportToPROJString(
+        PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+            .get());
+    auto expectedPROJString = "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=2 "
+                              "+x_0=3 +y_0=4 +datum=WGS84 +units=m +no_defs";
+    EXPECT_EQ(projString, expectedPROJString);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_polar_stereographic_scale_factor) {
+    auto wkt = "PROJCS[\"unknown\",\n"
+               "    GEOGCS[\"unknown\",\n"
+               "        DATUM[\"WGS_1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+               "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+               "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+               "        UNIT[\"degree\",0.0174532925199433,\n"
+               "            AUTHORITY[\"EPSG\",\"9122\"]]],\n"
+               "    PROJECTION[\"Polar_Stereographic\"],\n"
+               "    PARAMETER[\"latitude_of_origin\",90],\n"
+               "    PARAMETER[\"central_meridian\",2],\n"
+               "    PARAMETER[\"scale_factor\",0.99],\n"
+               "    PARAMETER[\"false_easting\",3],\n"
+               "    PARAMETER[\"false_northing\",4],\n"
+               "    UNIT[\"metre\",1,\n"
+               "        AUTHORITY[\"EPSG\",\"9001\"]]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto projString = crs->exportToPROJString(
+        PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+            .get());
+    auto expectedPROJString = "+proj=stere +lat_0=90 +lon_0=2 +k=0.99 +x_0=3 "
+                              "+y_0=4 +datum=WGS84 +units=m +no_defs";
+    EXPECT_EQ(projString, expectedPROJString);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_Spherical_Cross_Track_Height) {
+    auto wkt = "PROJCS[\"unknown\",\n"
+               "    GEOGCS[\"unknown\",\n"
+               "        DATUM[\"WGS_1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+               "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+               "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+               "        UNIT[\"degree\",0.0174532925199433,\n"
+               "            AUTHORITY[\"EPSG\",\"9122\"]]],\n"
+               "    PROJECTION[\"Spherical_Cross_Track_Height\"],\n"
+               "    PARAMETER[\"peg_point_latitude\",1],\n"
+               "    PARAMETER[\"peg_point_longitude\",2],\n"
+               "    PARAMETER[\"peg_point_heading\",3],\n"
+               "    PARAMETER[\"peg_point_height\",4],\n"
+               "    UNIT[\"metre\",1,\n"
+               "        AUTHORITY[\"EPSG\",\"9001\"]]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto projString = crs->exportToPROJString(
+        PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+            .get());
+    auto expectedPROJString = "+proj=sch +plat_0=1 +plon_0=2 +phdg_0=3 +h_0=4 "
+                              "+datum=WGS84 +units=m +no_defs";
+    EXPECT_EQ(projString, expectedPROJString);
 }
 
 // ---------------------------------------------------------------------------
@@ -3508,6 +3794,56 @@ TEST(wkt_parse, esri_projcs) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt1_esri_case_insensitive_names) {
+    auto wkt = "PROJCS[\"WGS_1984_UTM_Zone_31N\",GEOGCS[\"GCS_WGS_1984\","
+               "DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,"
+               "298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\","
+               "0.0174532925199433]],PROJECTION[\"transverse_mercator\"],"
+               "PARAMETER[\"false_easting\",500000.0],"
+               "PARAMETER[\"false_northing\",0.0],"
+               "PARAMETER[\"central_meridian\",3.0],"
+               "PARAMETER[\"scale_factor\",0.9996],"
+               "PARAMETER[\"latitude_of_origin\",0.0],UNIT[\"Meter\",1.0]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    int zone = 0;
+    bool north = false;
+    EXPECT_TRUE(crs->derivingConversion()->isUTM(zone, north));
+    EXPECT_EQ(zone, 31);
+    EXPECT_TRUE(north);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_esri_non_expected_param_name) {
+    // We try to be lax on parameter names.
+    auto wkt =
+        "PROJCS[\"WGS_1984_UTM_Zone_31N\",GEOGCS[\"GCS_WGS_1984\","
+        "DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,"
+        "298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\","
+        "0.0174532925199433]],PROJECTION[\"transverse_mercator\"],"
+        "PARAMETER[\"false_easting\",500000.0],"
+        "PARAMETER[\"false_northing\",0.0],"
+        "PARAMETER[\"longitude_of_center\",3.0]," // should be Central_Meridian
+        "PARAMETER[\"scale_factor\",0.9996],"
+        "PARAMETER[\"latitude_of_origin\",0.0],UNIT[\"Meter\",1.0]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    int zone = 0;
+    bool north = false;
+    EXPECT_TRUE(crs->derivingConversion()->isUTM(zone, north));
+    EXPECT_EQ(zone, 31);
+    EXPECT_TRUE(north);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, wkt1_esri_krovak_south_west) {
     auto wkt = "PROJCS[\"S-JTSK_Krovak\",GEOGCS[\"GCS_S_JTSK\","
                "DATUM[\"D_S_JTSK\","
@@ -5305,6 +5641,58 @@ TEST(io, projparse_longlat_a_rf) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_longlat_a_f_zero) {
+    auto obj =
+        PROJStringParser().createFromPROJString("+proj=longlat +a=2 +f=0");
+    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);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_longlat_a_f_non_zero) {
+    auto obj =
+        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);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_longlat_R) {
     auto obj = PROJStringParser().createFromPROJString("+proj=longlat +R=2");
     auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
@@ -5339,7 +5727,7 @@ TEST(io, projparse_longlat_pm_paris) {
     f->simulCurNodeHasId();
     crs->exportToWKT(f.get());
     auto expected = "GEODCRS[\"unknown\",\n"
-                    "    DATUM[\"Unknown based on WGS84 ellipsoid\",\n"
+                    "    DATUM[\"Unknown based on WGS 84 ellipsoid\",\n"
                     "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
                     "            LENGTHUNIT[\"metre\",1]]],\n"
                     "    PRIMEM[\"Paris\",2.5969213,\n"
@@ -5393,7 +5781,7 @@ TEST(io, projparse_longlat_pm_numeric) {
     f->simulCurNodeHasId();
     crs->exportToWKT(f.get());
     auto expected = "GEODCRS[\"unknown\",\n"
-                    "    DATUM[\"Unknown based on WGS84 ellipsoid\",\n"
+                    "    DATUM[\"Unknown based on WGS 84 ellipsoid\",\n"
                     "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
                     "            LENGTHUNIT[\"metre\",1]]],\n"
                     "    PRIMEM[\"unknown\",2.5,\n"
@@ -5411,6 +5799,21 @@ TEST(io, projparse_longlat_pm_numeric) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_longlat_pm_overriding_datum) {
+    // It is arguable that we allow the prime meridian of a datum defined by
+    // its name to be overriden, but this is found at least in a regression test
+    // of GDAL. So let's keep the ellipsoid part of the datum in that case and
+    // use the specified prime meridian.
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=WGS84 +pm=ferro");
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->datum()->nameStr(), "Unknown based on WGS 84 ellipsoid");
+    EXPECT_EQ(crs->datum()->primeMeridian()->nameStr(), "Ferro");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_longlat_complex) {
     std::string input =
         "+proj=pipeline +step +proj=longlat +ellps=clrk80ign "
@@ -5455,7 +5858,7 @@ TEST(io, projparse_longlat_towgs84_3_terms) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +ellps=GRS80 +towgs84=1.2,2,3,0,0,0,0");
+        "+proj=longlat +ellps=GRS80 +towgs84=1.2,2,3,0,0,0,0 +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -5490,7 +5893,7 @@ TEST(io, projparse_longlat_towgs84_7_terms) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +ellps=GRS80 +towgs84=1.2,2,3,4,5,6,7");
+        "+proj=longlat +ellps=GRS80 +towgs84=1.2,2,3,4,5,6,7 +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -5514,7 +5917,7 @@ TEST(io, projparse_longlat_nadgrids) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +ellps=GRS80 +nadgrids=foo.gsb");
+        "+proj=longlat +ellps=GRS80 +nadgrids=foo.gsb +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -5542,7 +5945,7 @@ TEST(io, projparse_longlat_geoidgrids) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +ellps=GRS80 +geoidgrids=foo.gtx +vunits=m");
+        "+proj=longlat +ellps=GRS80 +geoidgrids=foo.gtx +vunits=m +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -5839,6 +6242,13 @@ TEST(io, projparse_tmerc_south_oriented) {
         "            LENGTHUNIT[\"metre\",1]]]";
 
     EXPECT_EQ(f->toString(), expected);
+
+    obj = PROJStringParser().createFromPROJString(
+        "+proj=pipeline +step +proj=tmerc +x_0=1 +lat_0=1 +k_0=2 +step "
+        "+proj=axisswap +order=-1,-2");
+    crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    EXPECT_EQ(crs->derivingConversion()->method()->nameStr(),
+              "Transverse Mercator (South Orientated)");
 }
 
 // ---------------------------------------------------------------------------
@@ -6093,14 +6503,26 @@ TEST(io, projparse_etmerc) {
     auto obj = PROJStringParser().createFromPROJString("+proj=etmerc");
     auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
     ASSERT_TRUE(crs != nullptr);
-    WKTFormatterNNPtr f(WKTFormatter::create());
-    f->simulCurNodeHasId();
-    f->setMultiLine(false);
-    crs->exportToWKT(f.get());
-    auto wkt = f->toString();
-    EXPECT_TRUE(wkt.find("METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]]") !=
-                std::string::npos)
-        << wkt;
+    auto wkt2 = crs->exportToWKT(
+        &WKTFormatter::create()->simulCurNodeHasId().setMultiLine(false));
+    EXPECT_TRUE(
+        wkt2.find("METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]]") !=
+        std::string::npos)
+        << wkt2;
+
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 "
+        "+datum=WGS84 +units=m +no_defs");
+
+    auto wkt1 = crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=etmerc +lat_0=0 "
+                          "+lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m "
+                          "+no_defs\"]") != std::string::npos)
+        << wkt1;
 }
 
 // ---------------------------------------------------------------------------
@@ -6181,6 +6603,39 @@ TEST(io, projparse_merc_stere_polar_variant_A) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_merc_stere_polar_k_and_lat_ts) {
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=stere +lat_0=90 +lat_ts=90 +k=1");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto wkt = crs->exportToWKT(
+        &(WKTFormatter::create()->simulCurNodeHasId().setMultiLine(false)));
+    EXPECT_TRUE(
+        wkt.find(
+            "METHOD[\"Polar Stereographic (variant B)\",ID[\"EPSG\",9829]]") !=
+        std::string::npos)
+        << wkt;
+    EXPECT_TRUE(wkt.find("PARAMETER[\"Latitude of standard parallel\",90") !=
+                std::string::npos)
+        << wkt;
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 "
+        "+units=m +no_defs");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_merc_stere_polar_k_and_lat_ts_incompatible) {
+    EXPECT_THROW(PROJStringParser().createFromPROJString(
+                     "+proj=stere +lat_0=90 +lat_ts=70 +k=0.994"),
+                 ParsingException);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_merc_stere) {
     auto obj = PROJStringParser().createFromPROJString("+proj=stere +lat_0=30");
     auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
@@ -6239,7 +6694,8 @@ TEST(io, projparse_utm_south) {
 // ---------------------------------------------------------------------------
 
 TEST(io, projparse_non_earth_ellipsoid) {
-    std::string input("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=1");
+    std::string input("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=1 +units=m "
+                      "+no_defs");
     auto obj = PROJStringParser().createFromPROJString(input);
     auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
     ASSERT_TRUE(crs != nullptr);
@@ -6435,7 +6891,7 @@ TEST(io, projparse_projected_unknown) {
         "84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY["
         "\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\","
         "\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\","
-        "\"9122\"]],AXIS[\"Longitude\",EAST],AXIS[\"Latitude\",NORTH]],"
+        "\"9122\"]]],"
         "PROJECTION[\"custom_proj4\"],UNIT[\"metre\",1,AUTHORITY[\"EPSG\","
         "\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],EXTENSION["
         "\"PROJ4\",\"+proj=mbt_s +datum=WGS84 +unused_flag +lat_0=45 "
@@ -6456,7 +6912,7 @@ TEST(io, projparse_projected_unknown) {
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
         "+proj=mbt_s +unused_flag +lat_0=45 +lon_0=0 +k=1 +x_0=10 "
-        "+y_0=0 +datum=WGS84");
+        "+y_0=0 +datum=WGS84 +units=m +no_defs");
 
     EXPECT_EQ(
         crs->exportToPROJString(
@@ -6512,6 +6968,39 @@ TEST(io, projparse_geocent) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_geocent_towgs84) {
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=geocent +ellps=WGS84 +towgs84=1,2,3");
+    auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    WKTFormatterNNPtr f(WKTFormatter::create());
+    f->simulCurNodeHasId();
+    f->setMultiLine(false);
+    crs->exportToWKT(f.get());
+    auto wkt = f->toString();
+    EXPECT_TRUE(
+        wkt.find("METHOD[\"Geocentric translations (geocentric domain)") !=
+        std::string::npos)
+        << wkt;
+    EXPECT_TRUE(wkt.find("PARAMETER[\"X-axis translation\",1") !=
+                std::string::npos)
+        << wkt;
+    EXPECT_TRUE(wkt.find("PARAMETER[\"Y-axis translation\",2") !=
+                std::string::npos)
+        << wkt;
+    EXPECT_TRUE(wkt.find("PARAMETER[\"Z-axis translation\",3") !=
+                std::string::npos)
+        << wkt;
+
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=geocent +ellps=WGS84 +towgs84=1,2,3,0,0,0,0 +units=m +no_defs");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_cart_unit) {
     std::string input(
         "+proj=pipeline +step +proj=cart +ellps=WGS84 +step "
@@ -6544,6 +7033,48 @@ TEST(io, projparse_cart_unit_numeric) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_longlat_wktext) {
+    std::string input("+proj=longlat +foo +wktext");
+    auto obj = PROJStringParser().createFromPROJString(input);
+    auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        input);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_geocent_wktext) {
+    std::string input("+proj=geocent +foo +wktext");
+    auto obj = PROJStringParser().createFromPROJString(input);
+    auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        input);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_projected_wktext) {
+    std::string input("+proj=merc +foo +wktext");
+    auto obj = PROJStringParser().createFromPROJString(input);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        input);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_ob_tran_longlat) {
     std::string input(
         "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
@@ -6872,7 +7403,8 @@ TEST(io, projparse_longlat_title) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=longlat +ellps=intl +towgs84=109.753,-528.133,-362.244,0,0,0,0");
+        "+proj=longlat +ellps=intl +towgs84=109.753,-528.133,-362.244,0,0,0,0 "
+        "+no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -6895,7 +7427,7 @@ TEST(io, projparse_projected_title) {
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
         "+proj=utm +zone=43 +south +ellps=intl "
-        "+towgs84=109.753,-528.133,-362.244,0,0,0,0");
+        "+towgs84=109.753,-528.133,-362.244,0,0,0,0 +units=m +no_defs");
 }
 
 // ---------------------------------------------------------------------------
@@ -6992,6 +7524,10 @@ TEST(io, projparse_longlat_errors) {
                      "+proj=longlat +a=1 +rf=invalid"),
                  ParsingException);
 
+    EXPECT_THROW(PROJStringParser().createFromPROJString(
+                     "+proj=longlat +a=1 +f=invalid"),
+                 ParsingException);
+
     EXPECT_THROW(
         PROJStringParser().createFromPROJString("+proj=longlat +R=invalid"),
         ParsingException);
@@ -7005,6 +7541,9 @@ TEST(io, projparse_longlat_errors) {
     EXPECT_THROW(PROJStringParser().createFromPROJString("+proj=longlat +rf=1"),
                  ParsingException);
 
+    EXPECT_THROW(PROJStringParser().createFromPROJString("+proj=longlat +f=0"),
+                 ParsingException);
+
     EXPECT_THROW(
         PROJStringParser().createFromPROJString("+proj=longlat +pm=unknown"),
         ParsingException);
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 2d2688a8..e62eddf9 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -1185,13 +1185,13 @@ TEST(operation, tmerc_export) {
     auto conv = Conversion::createTransverseMercator(
         PropertyMap(), Angle(1), Angle(2), Scale(3), Length(4), Length(5));
     EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=tmerc +lat_0=1 +lon_0=2 +k_0=3 +x_0=4 +y_0=5");
+              "+proj=tmerc +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5");
 
     {
         auto formatter = PROJStringFormatter::create();
         formatter->setUseETMercForTMerc(true);
         EXPECT_EQ(conv->exportToPROJString(formatter.get()),
-                  "+proj=etmerc +lat_0=1 +lon_0=2 +k_0=3 +x_0=4 +y_0=5");
+                  "+proj=etmerc +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5");
     }
 
     EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()),
@@ -1270,7 +1270,7 @@ TEST(operation, tmerc_south_oriented_export) {
         PropertyMap(), Angle(1), Angle(2), Scale(3), Length(4), Length(5));
 
     EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=tmerc +axis=wsu +lat_0=1 +lon_0=2 +k_0=3 +x_0=4 +y_0=5");
+              "+proj=tmerc +axis=wsu +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5");
 
     EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()),
               "CONVERSION[\"Transverse Mercator (South Orientated)\",\n"
@@ -1329,7 +1329,7 @@ TEST(operation, tmerc_south_oriented_export) {
     EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
               "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=tmerc "
-              "+axis=wsu +lat_0=0 +lon_0=29 +k_0=1 +x_0=0 +y_0=0 +ellps=WGS84");
+              "+axis=wsu +lat_0=0 +lon_0=29 +k=1 +x_0=0 +y_0=0 +ellps=WGS84");
 }
 
 // ---------------------------------------------------------------------------
@@ -1986,7 +1986,7 @@ TEST(operation, createEquidistantCylindrical) {
         PropertyMap(), Angle(1), Angle(2), Length(3), Length(4));
 
     EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=eqc +lat_ts=1 +lon_0=2 +x_0=3 +y_0=4");
+              "+proj=eqc +lat_ts=1 +lat_0=0 +lon_0=2 +x_0=3 +y_0=4");
 
     EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()),
               "CONVERSION[\"Equidistant Cylindrical\",\n"
@@ -2022,7 +2022,7 @@ TEST(operation, createEquidistantCylindricalSpherical) {
         PropertyMap(), Angle(1), Angle(2), Length(3), Length(4));
 
     EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=eqc +lat_ts=1 +lon_0=2 +x_0=3 +y_0=4");
+              "+proj=eqc +lat_ts=1 +lat_0=0 +lon_0=2 +x_0=3 +y_0=4");
 
     EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()),
               "CONVERSION[\"Equidistant Cylindrical (Spherical)\",\n"
@@ -2053,6 +2053,28 @@ TEST(operation, createEquidistantCylindricalSpherical) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, equidistant_cylindrical_lat_0) {
+
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=eqc +ellps=sphere +lat_0=-2 +lat_ts=1 +lon_0=-10");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto wkt = crs->exportToWKT(WKTFormatter::create().get());
+    EXPECT_TRUE(wkt.find("PARAMETER[\"Latitude of natural origin\",-2") !=
+                std::string::npos)
+        << wkt;
+
+    auto projString = crs->exportToPROJString(
+        PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+            .get());
+    EXPECT_EQ(projString,
+              "+proj=eqc +lat_ts=1 +lat_0=-2 +lon_0=-10 +x_0=0 +y_0=0 "
+              "+ellps=sphere +units=m +no_defs");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, gall_export) {
 
     auto conv =
@@ -2173,10 +2195,18 @@ TEST(operation, geostationary_satellite_sweep_x_export) {
               "        LENGTHUNIT[\"metre\",1],\n"
               "        ID[\"EPSG\",8807]]]");
 
-    EXPECT_THROW(
-        conv->exportToWKT(
-            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
-        FormattingException);
+    auto crs = ProjectedCRS::create(
+        PropertyMap(), GeographicCRS::EPSG_4326, conv,
+        CartesianCS::createEastingNorthing(UnitOfMeasure::METRE));
+    auto wkt1 = crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+    EXPECT_TRUE(wkt1.find("PROJECTION[\"Geostationary_Satellite\"]") !=
+                std::string::npos)
+        << wkt1;
+    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=geos +sweep=x +lon_0=1 "
+                          "+h=2 +x_0=3 +y_0=4 +datum=WGS84 +units=m "
+                          "+no_defs\"]]") != std::string::npos)
+        << wkt1;
 }
 
 // ---------------------------------------------------------------------------
@@ -2410,7 +2440,7 @@ TEST(operation, hotine_oblique_mercator_two_point_natural_origin_export) {
         conv->exportToWKT(
             WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
         "PROJECTION[\"Hotine_Oblique_Mercator_Two_Point_Natural_Origin\"],\n"
-        "PARAMETER[\"latitude_of_origin\",1],\n"
+        "PARAMETER[\"latitude_of_center\",1],\n"
         "PARAMETER[\"latitude_of_point_1\",2],\n"
         "PARAMETER[\"longitude_of_point_1\",3],\n"
         "PARAMETER[\"latitude_of_point_2\",4],\n"
@@ -2435,12 +2465,12 @@ TEST(operation, imw_polyconic_export) {
               "    PARAMETER[\"Longitude of natural origin\",1,\n"
               "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
               "        ID[\"EPSG\",8802]],\n"
-              "    PARAMETER[\"Latitude of 1st standard parallel\",3,\n"
-              "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
-              "        ID[\"EPSG\",8823]],\n"
-              "    PARAMETER[\"Latitude of 2nd standard parallel\",4,\n"
-              "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
-              "        ID[\"EPSG\",8824]],\n"
+              "    PARAMETER[\"Latitude of 1st point\",3,\n"
+              "        ANGLEUNIT[\"degree\",0.0174532925199433,\n"
+              "            ID[\"EPSG\",9122]]],\n"
+              "    PARAMETER[\"Latitude of 2nd point\",4,\n"
+              "        ANGLEUNIT[\"degree\",0.0174532925199433,\n"
+              "            ID[\"EPSG\",9122]]],\n"
               "    PARAMETER[\"False easting\",5,\n"
               "        LENGTHUNIT[\"metre\",1],\n"
               "        ID[\"EPSG\",8806]],\n"
@@ -2453,8 +2483,8 @@ TEST(operation, imw_polyconic_export) {
             WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()),
         "PROJECTION[\"International_Map_of_the_World_Polyconic\"],\n"
         "PARAMETER[\"central_meridian\",1],\n"
-        "PARAMETER[\"standard_parallel_1\",3],\n"
-        "PARAMETER[\"standard_parallel_2\",4],\n"
+        "PARAMETER[\"Latitude_Of_1st_Point\",3],\n"
+        "PARAMETER[\"Latitude_Of_2nd_Point\",4],\n"
         "PARAMETER[\"false_easting\",5],\n"
         "PARAMETER[\"false_northing\",6]");
 }
@@ -2467,7 +2497,8 @@ TEST(operation, krovak_north_oriented_export) {
         Angle(78.5), Scale(0.9999), Length(5), Length(6));
 
     EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=krovak +lat_0=49.5 +lon_0=42.5 +k=0.9999 +x_0=5 +y_0=6");
+              "+proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.2881397222222 "
+              "+k=0.9999 +x_0=5 +y_0=6");
 
     EXPECT_EQ(
         conv->exportToWKT(WKTFormatter::create().get()),
@@ -2517,7 +2548,8 @@ TEST(operation, krovak_export) {
         Angle(78.5), Scale(0.9999), Length(5), Length(6));
 
     EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=krovak +axis=swu +lat_0=49.5 +lon_0=42.5 +k=0.9999 +x_0=5 "
+              "+proj=krovak +axis=swu +lat_0=49.5 +lon_0=42.5 "
+              "+alpha=30.2881397222222 +k=0.9999 +x_0=5 "
               "+y_0=6");
 
     EXPECT_EQ(
@@ -2706,6 +2738,36 @@ TEST(operation, wkt1_import_mercator_variant_A) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, wkt1_import_mercator_variant_A_that_is_variant_B) {
+    // Adresses https://trac.osgeo.org/gdal/ticket/3026
+    auto wkt = "PROJCS[\"test\",\n"
+               "    GEOGCS[\"WGS 84\",\n"
+               "        DATUM[\"WGS 1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563]],\n"
+               "        PRIMEM[\"Greenwich\",0],\n"
+               "        UNIT[\"degree\",0.0174532925199433]],\n"
+               "    PROJECTION[\"Mercator_1SP\"],\n"
+               "    PARAMETER[\"latitude_of_origin\",-1],\n"
+               "    PARAMETER[\"central_meridian\",2],\n"
+               "    PARAMETER[\"scale_factor\",1],\n"
+               "    PARAMETER[\"false_easting\",3],\n"
+               "    PARAMETER[\"false_northing\",4],\n"
+               "    UNIT[\"metre\",1]]";
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto conversion = crs->derivingConversion();
+    auto convRef = Conversion::createMercatorVariantB(
+        PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"), Angle(-1),
+        Angle(2), Length(3), Length(4));
+
+    EXPECT_TRUE(conversion->isEquivalentTo(convRef.get(),
+                                           IComparable::Criterion::EQUIVALENT));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, mercator_variant_B_export) {
     auto conv = Conversion::createMercatorVariantB(
         PropertyMap(), Angle(1), Angle(2), Length(3), Length(4));
@@ -2742,6 +2804,70 @@ TEST(operation, mercator_variant_B_export) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, odd_mercator_1sp_with_non_null_latitude) {
+    auto obj = WKTParser().createFromWKT(
+        "PROJCS[\"unnamed\",\n"
+        "    GEOGCS[\"WGS 84\",\n"
+        "        DATUM[\"WGS_1984\",\n"
+        "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+        "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+        "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+        "        PRIMEM[\"Greenwich\",0,\n"
+        "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+        "        UNIT[\"degree\",0.0174532925199433,\n"
+        "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+        "        AUTHORITY[\"EPSG\",\"4326\"]],\n"
+        "    PROJECTION[\"Mercator_1SP\"],\n"
+        "    PARAMETER[\"latitude_of_origin\",30],\n"
+        "    PARAMETER[\"central_meridian\",0],\n"
+        "    PARAMETER[\"scale_factor\",0.99],\n"
+        "    PARAMETER[\"false_easting\",0],\n"
+        "    PARAMETER[\"false_northing\",0],\n"
+        "    UNIT[\"metre\",1],\n"
+        "    AXIS[\"Easting\",EAST],\n"
+        "    AXIS[\"Northing\",NORTH]]");
+
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    EXPECT_THROW(crs->exportToPROJString(PROJStringFormatter::create().get()),
+                 FormattingException);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, odd_mercator_2sp_with_latitude_of_origin) {
+    auto obj = WKTParser().createFromWKT(
+        "PROJCS[\"unnamed\",\n"
+        "    GEOGCS[\"WGS 84\",\n"
+        "        DATUM[\"WGS_1984\",\n"
+        "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+        "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+        "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+        "        PRIMEM[\"Greenwich\",0,\n"
+        "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+        "        UNIT[\"degree\",0.0174532925199433,\n"
+        "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+        "        AUTHORITY[\"EPSG\",\"4326\"]],\n"
+        "    PROJECTION[\"Mercator_2SP\"],\n"
+        "    PARAMETER[\"standard_parallel_1\",30],\n"
+        "    PARAMETER[\"latitude_of_origin\",40],\n"
+        "    PARAMETER[\"central_meridian\",0],\n"
+        "    PARAMETER[\"false_easting\",0],\n"
+        "    PARAMETER[\"false_northing\",0],\n"
+        "    UNIT[\"metre\",1],\n"
+        "    AXIS[\"Easting\",EAST],\n"
+        "    AXIS[\"Northing\",NORTH]]");
+
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    EXPECT_THROW(crs->exportToPROJString(PROJStringFormatter::create().get()),
+                 FormattingException);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, webmerc_export) {
     auto conv = Conversion::createPopularVisualisationPseudoMercator(
         PropertyMap(), Angle(0), Angle(2), Length(3), Length(4));
@@ -2790,8 +2916,6 @@ TEST(operation, webmerc_export) {
         "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
         "        UNIT[\"degree\",0.0174532925199433,\n"
         "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
-        "        AXIS[\"Latitude\",NORTH],\n"
-        "        AXIS[\"Longitude\",EAST],\n"
         "        AUTHORITY[\"EPSG\",\"4326\"]],\n"
         "    PROJECTION[\"Mercator_1SP\"],\n"
         "    PARAMETER[\"central_meridian\",2],\n"
@@ -2994,9 +3118,7 @@ TEST(operation, webmerc_import_from_WKT2_EPSG_3785_deprecated) {
         "            SPHEROID[\"Popular Visualisation Sphere\",6378137,0],\n"
         "            TOWGS84[0,0,0,0,0,0,0]],\n"
         "        PRIMEM[\"Greenwich\",0],\n"
-        "        UNIT[\"degree\",0.0174532925199433],\n"
-        "        AXIS[\"Latitude\",NORTH],\n"
-        "        AXIS[\"Longitude\",EAST]],\n"
+        "        UNIT[\"degree\",0.0174532925199433]],\n"
         "    PROJECTION[\"Mercator_1SP\"],\n"
         "    PARAMETER[\"central_meridian\",0],\n"
         "    PARAMETER[\"scale_factor\",1],\n"
@@ -3730,7 +3852,7 @@ TEST(operation, conversion_inverse) {
               "        ID[\"EPSG\",8807]]]");
 
     EXPECT_EQ(inv->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +inv +proj=tmerc +lat_0=1 +lon_0=2 +k_0=3 "
+              "+proj=pipeline +step +inv +proj=tmerc +lat_0=1 +lon_0=2 +k=3 "
               "+x_0=4 +y_0=5");
 
     EXPECT_TRUE(inv->isEquivalentTo(inv.get()));
@@ -5463,7 +5585,7 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) {
                                                                     compound2);
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +inv +proj=tmerc +lat_0=1 +lon_0=2 +k_0=3 "
+              "+proj=pipeline +step +inv +proj=tmerc +lat_0=1 +lon_0=2 +k=3 "
               "+x_0=4 +y_0=5 +ellps=WGS84 +step "
               "+proj=vgridshift +grids=bla.gtx +multiplier=0.001 +step "
               "+proj=utm +zone=32 "
@@ -5473,7 +5595,7 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) {
         formatter->setUseETMercForTMerc(true);
         EXPECT_EQ(op->exportToPROJString(formatter.get()),
                   "+proj=pipeline +step +inv +proj=etmerc +lat_0=1 +lon_0=2 "
-                  "+k_0=3 +x_0=4 +y_0=5 +ellps=WGS84 +step "
+                  "+k=3 +x_0=4 +y_0=5 +ellps=WGS84 +step "
                   "+proj=vgridshift +grids=bla.gtx +multiplier=0.001 +step "
                   "+proj=utm +zone=32 "
                   "+ellps=WGS84");
@@ -5485,7 +5607,7 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) {
                   "+proj=pipeline +step +inv +proj=utm +zone=32 +ellps=WGS84 "
                   "+step +inv +proj=vgridshift +grids=bla.gtx "
                   "+multiplier=0.001 +step +proj=etmerc +lat_0=1 +lon_0=2 "
-                  "+k_0=3 +x_0=4 +y_0=5 +ellps=WGS84");
+                  "+k=3 +x_0=4 +y_0=5 +ellps=WGS84");
     }
 
     auto opInverse = CoordinateOperationFactory::create()->createOperation(