Fix ingestion of +proj=cea with +k_0

Fixes #1881

Digging into the implementation of proj=cea, it appears that
k_0 and lat_ts are intended to be exclusive ways of specifying the
same concept. EPSG only models the variant using lat_s.
So if k_0 is found and lat_ts is absent, compute the equivalent
value of lat_ts from k_0.

Note: k_0 should normally be in the [0,1] range. In case creative users
would use something outside, we raise an exception, even if the cea
implementation could potentially deal with any k_0 value. Hopefully
this is a (reasonable) limitation that will address nominal use cases.

3 files changed, 60 insertions, 15 deletions

diff --git a/docs/source/operations/projections/cea.rst b/docs/source/operations/projections/cea.rst
index 4542fe71..80fe1df3 100644
--- a/docs/source/operations/projections/cea.rst
+++ b/docs/source/operations/projections/cea.rst
@@ -29,3 +29,9 @@ Parameters
 .. include:: ../options/x_0.rst
 
 .. include:: ../options/y_0.rst
+
+.. note::
+
+    ``lat_ts`` and ``k_0`` are mutually exclusive. If ``lat_ts``
+    is specified, it is equivalent to setting ``k_0`` to
+    :math:`\frac{\cos \phi_{ts}}{\sqrt{1 - e^2 \sin^2 \phi_{ts}}}`
\ No newline at end of file
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index d88581b4..0fc11f6f 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -8736,22 +8736,43 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
 
             } else if (param->unit_type == UnitOfMeasure::Type::SCALE) {
                 value = 1;
-            } else {
-                // For omerc, if gamma is missing, the default value is
-                // alpha
-                if (step.name == "omerc" && proj_name == "gamma") {
-                    paramValue = &getParamValue(step, "alpha");
-                    if (!paramValue->empty()) {
-                        value = getAngularValue(*paramValue);
+            }
+            // For omerc, if gamma is missing, the default value is
+            // alpha
+            else if (step.name == "omerc" && proj_name == "gamma") {
+                paramValue = &getParamValue(step, "alpha");
+                if (!paramValue->empty()) {
+                    value = getAngularValue(*paramValue);
+                }
+            } else if (step.name == "krovak") {
+                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 (step.name == "cea" && proj_name == "lat_ts") {
+                paramValue = &getParamValueK(step);
+                if (!paramValue->empty()) {
+                    bool hasError = false;
+                    const double k = getNumericValue(*paramValue, &hasError);
+                    if (hasError) {
+                        throw ParsingException("invalid value for k/k_0");
                     }
-                } else if (step.name == "krovak") {
-                    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;
+                    if (k >= 0 && k <= 1) {
+                        const double es =
+                            geogCRS->ellipsoid()->squaredEccentricity();
+                        if (es < 0) {
+                            throw ParsingException("Invalid flattening");
+                        }
+                        value =
+                            Angle(acos(k * sqrt((1 - es) / (1 - k * k * es))),
+                                  UnitOfMeasure::RADIAN)
+                                .convertToUnit(UnitOfMeasure::DEGREE);
+                    } else {
+                        throw ParsingException("k/k_0 should be in [0,1]");
                     }
                 }
             }
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index fd38847c..cc780be1 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -7968,6 +7968,24 @@ TEST(io, projparse_cea_ellipsoidal) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_cea_ellipsoidal_with_k_0) {
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=cea +ellps=GRS80 +k_0=0.99 +type=crs");
+    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("PARAMETER[\"Latitude of 1st standard parallel\",8.1365") !=
+        std::string::npos)
+        << wkt;
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_geos_sweep_x) {
     auto obj = PROJStringParser().createFromPROJString(
         "+proj=geos +sweep=x +h=1 +type=crs");