createFromProjString(): handle default parameters of '+krovak +type=crs', and handle +czech correctly (fixes #2199) (#2200)

* createFromProjString(): handle default parameters of '+krovak +type=crs', and handle +czech correctly (fixes #2199)

* Doc: add default values for Krovak

5 files changed, 88 insertions, 45 deletions

diff --git a/docs/source/operations/projections/krovak.rst b/docs/source/operations/projections/krovak.rst
index e8313e7b..59f66b04 100644
--- a/docs/source/operations/projections/krovak.rst
+++ b/docs/source/operations/projections/krovak.rst
@@ -33,20 +33,31 @@ Parameters
 
 .. note:: All parameters are optional for the Krovak projection.
 
+          The latitude of pseudo standard parallel is hardcoded to 78.5° and
+          the ellipsoid to Bessel.
+
 .. option:: +czech
 
     Reverse the sign of the output coordinates, as is tradition in the
     Czech Republic.
 
-.. include:: ../options/lon_0.rst
+.. option:: +lon_0=<value>
+
+    Longitude of projection center.
+
+    *Defaults to 24°50' (24.8333333333333)*
+
+.. option:: +lat_0=<value>
+
+    Latitude of projection center.
 
-.. include:: ../options/lat_0.rst
+    *Defaults to 49.5*
 
 .. option:: +k_0=<value>
 
     Scale factor. Determines scale factor used in the projection.
 
-    *Defaults to 0.9999.*
+    *Defaults to 0.9999*
 
 .. include:: ../options/x_0.rst
 
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 21754a44..a255b959 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -8089,7 +8089,7 @@ std::string PROJStringParser::Private::guessBodyName(double a) {
 GeodeticReferenceFrameNNPtr
 PROJStringParser::Private::buildDatum(Step &step, const std::string &title) {
 
-    const auto &ellpsStr = getParamValue(step, "ellps");
+    std::string ellpsStr = getParamValue(step, "ellps");
     const auto &datumStr = getParamValue(step, "datum");
     const auto &RStr = getParamValue(step, "R");
     const auto &aStr = getParamValue(step, "a");
@@ -8106,6 +8106,11 @@ PROJStringParser::Private::buildDatum(Step &step, const std::string &title) {
         !RStr.empty() || !aStr.empty() || !bStr.empty() || !rfStr.empty() ||
         !fStr.empty() || !esStr.empty() || !eStr.empty();
 
+    if (!numericParamPresent && ellpsStr.empty() && datumStr.empty() &&
+        step.name == "krovak") {
+        ellpsStr = "bessel";
+    }
+
     PrimeMeridianNNPtr pm(buildPrimeMeridian(step));
     PropertyMap grfMap;
 
@@ -8504,6 +8509,9 @@ PROJStringParser::Private::processAxisSwap(Step &step,
                 throw ParsingException("Unhandled order=" + orderStr);
             }
         }
+    } else if (step.name == "krovak" && hasParamValue(step, "czech")) {
+        axis[0] = west;
+        axis[1] = south;
     }
     return axis;
 }
@@ -8841,9 +8849,14 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                 Step::KeyValue("lonc", getParamValue(step, "lon_0")));
         }
     } else if (step.name == "krovak" &&
-               ((getParamValue(step, "axis") == "swu" && iAxisSwap < 0) ||
+               ((iAxisSwap < 0 && getParamValue(step, "axis") == "swu" &&
+                 !hasParamValue(step, "czech")) ||
                 (iAxisSwap > 0 &&
-                 getParamValue(steps_[iAxisSwap], "order") == "-2,-1"))) {
+                 getParamValue(steps_[iAxisSwap], "order") == "-2,-1" &&
+                 !hasParamValue(step, "czech")))) {
+        mapping = getMapping(EPSG_CODE_METHOD_KROVAK);
+    } else if (step.name == "krovak" && iAxisSwap < 0 &&
+               hasParamValue(step, "czech") && !hasParamValue(step, "axis")) {
         mapping = getMapping(EPSG_CODE_METHOD_KROVAK);
     } else if (step.name == "merc") {
         if (hasParamValue(step, "a") && hasParamValue(step, "b") &&
@@ -8992,9 +9005,6 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                 if (hasError) {
                     throw ParsingException("invalid value for " + proj_name);
                 }
-
-            } else if (param->unit_type == UnitOfMeasure::Type::SCALE) {
-                value = 1;
             }
             // For omerc, if gamma is missing, the default value is
             // alpha
@@ -9004,13 +9014,24 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                     value = getAngularValue(*paramValue);
                 }
             } else if (step.name == "krovak") {
+                // Keep it in sync with defaults of krovak.cpp
                 if (param->epsg_code ==
-                    EPSG_CODE_PARAMETER_COLATITUDE_CONE_AXIS) {
+                    EPSG_CODE_PARAMETER_LATITUDE_PROJECTION_CENTRE) {
+                    value = 49.5;
+                } else if (param->epsg_code ==
+                           EPSG_CODE_PARAMETER_LONGITUDE_OF_ORIGIN) {
+                    value = 24.833333333333333333;
+                } else if (param->epsg_code ==
+                           EPSG_CODE_PARAMETER_COLATITUDE_CONE_AXIS) {
                     value = 30.28813975277777776;
                 } else if (
                     param->epsg_code ==
                     EPSG_CODE_PARAMETER_LATITUDE_PSEUDO_STANDARD_PARALLEL) {
                     value = 78.5;
+                } else if (
+                    param->epsg_code ==
+                    EPSG_CODE_PARAMETER_SCALE_FACTOR_PSEUDO_STANDARD_PARALLEL) {
+                    value = 0.9999;
                 }
             } else if (step.name == "cea" && proj_name == "lat_ts") {
                 paramValue = &getParamValueK(step);
@@ -9034,6 +9055,8 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                         throw ParsingException("k/k_0 should be in [0,1]");
                     }
                 }
+            } else if (param->unit_type == UnitOfMeasure::Type::SCALE) {
+                value = 1;
             }
 
             PropertyMap propertiesParameter;
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index ad32d98b..270e7f95 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -2894,17 +2894,18 @@ TEST(crs, Krovak_North_Orientated_as_WKT1_ESRI) {
     ASSERT_TRUE(crs != nullptr);
 
     auto expected = "PROJCS[\"unknown\",GEOGCS[\"GCS_unknown\","
-                    "DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\","
-                    "6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],"
+                    "DATUM[\"D_Unknown_based_on_Bessel_1841_ellipsoid\","
+                    "SPHEROID[\"Bessel_1841\",6377397.155,299.1528128]],"
+                    "PRIMEM[\"Greenwich\",0.0],"
                     "UNIT[\"Degree\",0.0174532925199433]],"
                     "PROJECTION[\"Krovak\"],"
                     "PARAMETER[\"False_Easting\",0.0],"
                     "PARAMETER[\"False_Northing\",0.0],"
                     "PARAMETER[\"Pseudo_Standard_Parallel_1\",78.5],"
-                    "PARAMETER[\"Scale_Factor\",1.0],"
+                    "PARAMETER[\"Scale_Factor\",0.9999],"
                     "PARAMETER[\"Azimuth\",30.2881397527778],"
-                    "PARAMETER[\"Longitude_Of_Center\",0.0],"
-                    "PARAMETER[\"Latitude_Of_Center\",0.0],"
+                    "PARAMETER[\"Longitude_Of_Center\",24.8333333333333],"
+                    "PARAMETER[\"Latitude_Of_Center\",49.5],"
                     "PARAMETER[\"X_Scale\",-1.0],"
                     "PARAMETER[\"Y_Scale\",1.0],"
                     "PARAMETER[\"XY_Plane_Rotation\",90.0],"
@@ -2926,17 +2927,18 @@ TEST(crs, Krovak_as_WKT1_ESRI) {
     ASSERT_TRUE(crs != nullptr);
 
     auto expected = "PROJCS[\"unknown\",GEOGCS[\"GCS_unknown\","
-                    "DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\","
-                    "6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],"
+                    "DATUM[\"D_Unknown_based_on_Bessel_1841_ellipsoid\","
+                    "SPHEROID[\"Bessel_1841\",6377397.155,299.1528128]],"
+                    "PRIMEM[\"Greenwich\",0.0],"
                     "UNIT[\"Degree\",0.0174532925199433]],"
                     "PROJECTION[\"Krovak\"],"
                     "PARAMETER[\"False_Easting\",0.0],"
                     "PARAMETER[\"False_Northing\",0.0],"
                     "PARAMETER[\"Pseudo_Standard_Parallel_1\",78.5],"
-                    "PARAMETER[\"Scale_Factor\",1.0],"
+                    "PARAMETER[\"Scale_Factor\",0.9999],"
                     "PARAMETER[\"Azimuth\",30.2881397527778],"
-                    "PARAMETER[\"Longitude_Of_Center\",0.0],"
-                    "PARAMETER[\"Latitude_Of_Center\",0.0],"
+                    "PARAMETER[\"Longitude_Of_Center\",24.8333333333333],"
+                    "PARAMETER[\"Latitude_Of_Center\",49.5],"
                     "PARAMETER[\"X_Scale\",1.0],"
                     "PARAMETER[\"Y_Scale\",1.0],"
                     "PARAMETER[\"XY_Plane_Rotation\",0.0],"
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index f829bf91..0570bb7e 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -1220,7 +1220,7 @@ TEST(wkt_parse, wkt1_krovak_south_west) {
         "    PROJECTION[\"Krovak\"],"
         "    PARAMETER[\"latitude_of_center\",49.5],"
         "    PARAMETER[\"longitude_of_center\",24.83333333333333],"
-        "    PARAMETER[\"azimuth\",30.28813972222222],"
+        "    PARAMETER[\"azimuth\",30.2881397527778],"
         "    PARAMETER[\"pseudo_standard_parallel_1\",78.5],"
         "    PARAMETER[\"scale_factor\",0.9999],"
         "    PARAMETER[\"false_easting\",0],"
@@ -1254,7 +1254,7 @@ TEST(wkt_parse, wkt1_krovak_south_west) {
         "        PARAMETER[\"Longitude of origin\",24.8333333333333,\n"
         "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
         "            ID[\"EPSG\",8833]],\n"
-        "        PARAMETER[\"Co-latitude of cone axis\",30.2881397222222,\n"
+        "        PARAMETER[\"Co-latitude of cone axis\",30.2881397527778,\n"
         "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
         "            ID[\"EPSG\",1036]],\n"
         "        PARAMETER[\"Latitude of pseudo standard parallel\",78.5,\n"
@@ -1284,37 +1284,44 @@ TEST(wkt_parse, wkt1_krovak_south_west) {
     auto projString =
         crs->exportToPROJString(PROJStringFormatter::create().get());
     auto expectedPROJString = "+proj=krovak +axis=swu +lat_0=49.5 "
-                              "+lon_0=24.8333333333333 +alpha=30.2881397222222 "
+                              "+lon_0=24.8333333333333 +alpha=30.2881397527778 "
                               "+k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m "
                               "+no_defs +type=crs";
     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());
+    auto crs2 = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs2 != nullptr);
+    auto wkt2 = crs2->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)
         << wkt2;
     EXPECT_TRUE(
-        wkt2.find("PARAMETER[\"Co-latitude of cone axis\",30.2881397222222,") !=
+        wkt2.find("PARAMETER[\"Co-latitude of cone axis\",30.2881397527778,") !=
         std::string::npos)
         << wkt2;
-    EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
+    EXPECT_EQ(crs2->exportToPROJString(PROJStringFormatter::create().get()),
+              expectedPROJString);
+
+    obj = PROJStringParser().createFromPROJString(
+        "+proj=krovak +czech +type=crs");
+    crs2 = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs2 != nullptr);
+    EXPECT_EQ(crs2->exportToPROJString(PROJStringFormatter::create().get()),
               expectedPROJString);
 
     obj = PROJStringParser().createFromPROJString(
         "+type=crs +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 "
+        "+lon_0=24.8333333333333 +alpha=30.2881397527778 "
         "+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());
+    crs2 = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs2 != nullptr);
+    wkt2 = crs2->exportToWKT(WKTFormatter::create().get());
     EXPECT_TRUE(wkt2.find("METHOD[\"Krovak\"") != std::string::npos) << wkt2;
 }
 
@@ -1336,7 +1343,7 @@ TEST(wkt_parse, wkt1_krovak_north_oriented) {
         "    PROJECTION[\"Krovak\"],"
         "    PARAMETER[\"latitude_of_center\",49.5],"
         "    PARAMETER[\"longitude_of_center\",24.83333333333333],"
-        "    PARAMETER[\"azimuth\",30.28813972222222],"
+        "    PARAMETER[\"azimuth\",30.2881397527778],"
         "    PARAMETER[\"pseudo_standard_parallel_1\",78.5],"
         "    PARAMETER[\"scale_factor\",0.9999],"
         "    PARAMETER[\"false_easting\",0],"
@@ -1372,7 +1379,7 @@ TEST(wkt_parse, wkt1_krovak_north_oriented) {
         "        PARAMETER[\"Longitude of origin\",24.8333333333333,\n"
         "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
         "            ID[\"EPSG\",8833]],\n"
-        "        PARAMETER[\"Co-latitude of cone axis\",30.2881397222222,\n"
+        "        PARAMETER[\"Co-latitude of cone axis\",30.2881397527778,\n"
         "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
         "            ID[\"EPSG\",1036]],\n"
         "        PARAMETER[\"Latitude of pseudo standard parallel\",78.5,\n"
@@ -1399,7 +1406,7 @@ TEST(wkt_parse, wkt1_krovak_north_oriented) {
 
     EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=krovak +lat_0=49.5 +lon_0=24.8333333333333 "
-              "+alpha=30.2881397222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel "
+              "+alpha=30.2881397527778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel "
               "+units=m +no_defs +type=crs");
 }
 
@@ -5487,7 +5494,7 @@ TEST(wkt_parse,
                "UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Krovak\"],"
                "PARAMETER[\"latitude_of_center\",49.5],"
                "PARAMETER[\"longitude_of_center\",24.83333333333333],"
-               "PARAMETER[\"azimuth\",30.28813972222222],"
+               "PARAMETER[\"azimuth\",30.2881397527778],"
                "PARAMETER[\"pseudo_standard_parallel_1\",78.5],"
                "PARAMETER[\"scale_factor\",0.9999],"
                "PARAMETER[\"false_easting\",0],"
@@ -5521,7 +5528,7 @@ TEST(wkt_parse,
         "        PARAMETER[\"Longitude of origin\",24.8333333333333,\n"
         "            ANGLEUNIT[\"Degree\",0.0174532925199433],\n"
         "            ID[\"EPSG\",8833]],\n"
-        "        PARAMETER[\"Co-latitude of cone axis\",30.2881397222222,\n"
+        "        PARAMETER[\"Co-latitude of cone axis\",30.2881397527778,\n"
         "            ANGLEUNIT[\"Degree\",0.0174532925199433],\n"
         "            ID[\"EPSG\",1036]],\n"
         "        PARAMETER[\"Latitude of pseudo standard parallel\",78.5,\n"
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 4aa9d01e..a655aff0 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -2653,12 +2653,12 @@ TEST(operation, imw_polyconic_export) {
 
 TEST(operation, krovak_north_oriented_export) {
     auto conv = Conversion::createKrovakNorthOriented(
-        PropertyMap(), Angle(49.5), Angle(42.5), Angle(30.28813972222222),
+        PropertyMap(), Angle(49.5), Angle(42.5), Angle(30.2881397527778),
         Angle(78.5), Scale(0.9999), Length(5), Length(6));
     EXPECT_TRUE(conv->validateParameters().empty());
 
     EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.2881397222222 "
+              "+proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.2881397527778 "
               "+k=0.9999 +x_0=5 +y_0=6");
 
     EXPECT_EQ(
@@ -2672,7 +2672,7 @@ TEST(operation, krovak_north_oriented_export) {
         "    PARAMETER[\"Longitude of origin\",42.5,\n"
         "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
         "        ID[\"EPSG\",8833]],\n"
-        "    PARAMETER[\"Co-latitude of cone axis\",30.2881397222222,\n"
+        "    PARAMETER[\"Co-latitude of cone axis\",30.2881397527778,\n"
         "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
         "        ID[\"EPSG\",1036]],\n"
         "    PARAMETER[\"Latitude of pseudo standard parallel\",78.5,\n"
@@ -2694,7 +2694,7 @@ TEST(operation, krovak_north_oriented_export) {
         "PROJECTION[\"Krovak\"],\n"
         "PARAMETER[\"latitude_of_center\",49.5],\n"
         "PARAMETER[\"longitude_of_center\",42.5],\n"
-        "PARAMETER[\"azimuth\",30.2881397222222],\n"
+        "PARAMETER[\"azimuth\",30.2881397527778],\n"
         "PARAMETER[\"pseudo_standard_parallel_1\",78.5],\n"
         "PARAMETER[\"scale_factor\",0.9999],\n"
         "PARAMETER[\"false_easting\",5],\n"
@@ -2705,13 +2705,13 @@ TEST(operation, krovak_north_oriented_export) {
 
 TEST(operation, krovak_export) {
     auto conv = Conversion::createKrovak(
-        PropertyMap(), Angle(49.5), Angle(42.5), Angle(30.28813972222222),
+        PropertyMap(), Angle(49.5), Angle(42.5), Angle(30.2881397527778),
         Angle(78.5), Scale(0.9999), Length(5), Length(6));
     EXPECT_TRUE(conv->validateParameters().empty());
 
     EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=krovak +axis=swu +lat_0=49.5 +lon_0=42.5 "
-              "+alpha=30.2881397222222 +k=0.9999 +x_0=5 "
+              "+alpha=30.2881397527778 +k=0.9999 +x_0=5 "
               "+y_0=6");
 
     EXPECT_EQ(
@@ -2725,7 +2725,7 @@ TEST(operation, krovak_export) {
         "    PARAMETER[\"Longitude of origin\",42.5,\n"
         "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
         "        ID[\"EPSG\",8833]],\n"
-        "    PARAMETER[\"Co-latitude of cone axis\",30.2881397222222,\n"
+        "    PARAMETER[\"Co-latitude of cone axis\",30.2881397527778,\n"
         "        ANGLEUNIT[\"degree\",0.0174532925199433],\n"
         "        ID[\"EPSG\",1036]],\n"
         "    PARAMETER[\"Latitude of pseudo standard parallel\",78.5,\n"
@@ -2747,7 +2747,7 @@ TEST(operation, krovak_export) {
         "PROJECTION[\"Krovak\"],\n"
         "PARAMETER[\"latitude_of_center\",49.5],\n"
         "PARAMETER[\"longitude_of_center\",42.5],\n"
-        "PARAMETER[\"azimuth\",30.2881397222222],\n"
+        "PARAMETER[\"azimuth\",30.2881397527778],\n"
         "PARAMETER[\"pseudo_standard_parallel_1\",78.5],\n"
         "PARAMETER[\"scale_factor\",0.9999],\n"
         "PARAMETER[\"false_easting\",5],\n"