Add +proj=mod_krovak projection method for Modified Krovak that appli…

…es to S-JTSK/05 in the Czech Republic Validated with the test point given in EPSG Guidance Note 7-2 Also validated with CUZK online calculator at https://geoportal.cuzk.cz/(S(g4ipaut0ckzb5lellu153lvf))/Default.aspx?head_tab=sekce-01-gp&mode=TextMeta&text=wcts&menu=19 Unfortunately this calculator does not offer 'elementary' transformations, and cannot just validate Modified Krovak alone. The best we can do is to use ETRS89 <--> S-JTSK/05 / Krovak modified, which involves a Helmert transformation as well for the datum transformation. * Validation of reverse modified Krovak: Online calculator: - Souřadnice/input: -5568990.91 -6050538.71 100 - Source CRS: S-JTSK/05 + Bpv (-Y-XH /east-north) - Destination CRS: ETRS89 (BLh / DEG) - Výsledek/result: 50.208297081 16.848326835 143.762 With PROJ: echo -5568990.91 -6050538.71 100 | bin/cs2cs -d 9 EPSG:5516 EPSG:4937 --3d 50.208297081 16.848326833 143.143233569 Note the different on longitude is totally neglectable: 0.1 mm, as computed by: $ echo 50.208297081 16.848326835 50.208297081 16.848326833 | bin/geod -F "%.4f" -I +ellps=GRS80 -90d 90d 0.0001 and that when outputing Krovak modified, the online calculator displays Y/X with a precision of 1 cm. The difference on Z is significant, but can be explained from the fact that the online calculator using Baltic Height for S-JTSK/05 and thus apply a Baltic Height<-->ETRS89 geoid that we don't apply here. * Validation of forward Krovak using same point: Online calculator: - Souřadnice/input: 50.208297081 16.848326835 143.762 - Source: ETRS89 (BLh / DEG) - Destination: S-JTSK/05 + Bpv (-Y-XH /east-north) - Výsledek/result: -5568990.91 -6050538.71 100 With PROJ: $ echo 50.208297081 16.848326835 143.762 | bin/cs2cs -d 2 EPSG:4937 EPSG:5516 --3d -5568990.91 -6050538.71 100.62
rouault · Feb 1, 2024 · 828a472 · 828a472
1 parent 3cc757f
commit 828a472
Show file tree

Hide file tree

Showing 11 changed files with 566 additions and 42 deletions.
diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst
@@ -93,6 +93,7 @@ Projections map the spherical 3D space to a flat 2D space.
    mil_os
    mill
    misrsom
+   mod_krovak
    moll
    murd1
    murd2

diff --git a/docs/source/operations/projections/krovak.rst b/docs/source/operations/projections/krovak.rst
@@ -5,11 +5,12 @@ Krovak
 ********************************************************************************
 
 +---------------------+----------------------------------------------------------+
-| **Classification**  | Conical                                                  |
+| **Classification**  | Conformal Conical                                        |
 +---------------------+----------------------------------------------------------+
 | **Available forms** | Forward and inverse, spherical and ellipsoidal           |
 +---------------------+----------------------------------------------------------+
-| **Defined area**    | Global, but more accurate around Czechoslovakia          |
+| **Defined area**    | Global, but more accurate around Czech Republic and      |
+|                     | Slovakia                                                 |
 +---------------------+----------------------------------------------------------+
 | **Alias**           | krovak                                                   |
 +---------------------+----------------------------------------------------------+
@@ -28,6 +29,16 @@ Krovak
 
    proj-string: ``+proj=krovak``
 
+By default, coordinates in the forward direction are output in easting, northing,
+and negative in the Czech Republic and Slovakia, with absolute value of
+easting/westing being smaller than absolute value of northing/southing.
+
+See also :ref:`mod_krovak` for a variation of Krovak used with the S-JTSK/05 datum
+in the Czech Republic.
+
+.. note:: Before PROJ 9.4, using other values for x_0 or y_0 than the default 0
+          would lead to incorrect results when not using the ``+czech`` switch.
+
 Parameters
 ################################################################################
 
@@ -39,7 +50,8 @@ Parameters
 .. option:: +czech
 
     Reverse the sign of the output coordinates, as is tradition in the
-    Czech Republic.
+    Czech Republic, to be westing, southing (positive values in Czech Republic
+    and Slovakia).
 
 .. option:: +lon_0=<value>
 

diff --git a/docs/source/operations/projections/mod_krovak.rst b/docs/source/operations/projections/mod_krovak.rst
@@ -0,0 +1,77 @@
+.. _mod_krovak:
+
+********************************************************************************
+Modified Krovak
+********************************************************************************
+
+.. versionadded:: 9.4.0
+
++---------------------+----------------------------------------------------------+
+| **Classification**  | Conical                                                  |
++---------------------+----------------------------------------------------------+
+| **Available forms** | Forward and inverse, spherical and ellipsoidal           |
++---------------------+----------------------------------------------------------+
+| **Defined area**    | Czech Republic                                           |
++---------------------+----------------------------------------------------------+
+| **Alias**           | mod_krovak                                               |
++---------------------+----------------------------------------------------------+
+| **Domain**          | 2D                                                       |
++---------------------+----------------------------------------------------------+
+| **Input type**      | Geodetic coordinates                                     |
++---------------------+----------------------------------------------------------+
+| **Output type**     | Projected coordinates                                    |
++---------------------+----------------------------------------------------------+
+
+
+.. figure:: ./images/krovak.png
+   :width: 500 px
+   :align: center
+   :alt:   Modified Krovak
+
+   proj-string: ``+proj=mod_krovak``
+
+Modified Krovak builts upon traditional :ref:`krovak`, with corrective terms that
+are better suited when using it with the S-JTSK/05 datum. This method is specific
+to the Czech Republic. Due to the corrective terms, this projection method is
+no longer strictly conformal.
+
+By default, coordinates in the forward direction are output in easting, northing,
+and negative in the Czech Republic, with absolute value of easting/westing
+being smaller than absolute value of northing/southing.
+To distinguish it from regular Krovak, the usual value for ``+x_0`` and ``+y_0``
+in Modified Krovak is typically 5,000,000.
+
+Parameters
+################################################################################
+
+.. note:: All parameters are optional for the Modified 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, to be westing, southing (positive values in Czech Republic)
+
+.. option:: +lon_0=<value>
+
+    Longitude of projection center.
+
+    *Defaults to 24°50' (24.8333333333333)*
+
+.. option:: +lat_0=<value>
+
+    Latitude of projection center.
+
+    *Defaults to 49.5*
+
+.. option:: +k_0=<value>
+
+    Scale factor. Determines scale factor used in the projection.
+
+    *Defaults to 0.9999*
+
+.. include:: ../options/x_0.rst
+
+.. include:: ../options/y_0.rst
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
@@ -10635,7 +10635,7 @@ PROJStringParser::Private::buildDatum(Step &step, const std::string &title) {
         !fStr.empty() || !esStr.empty() || !eStr.empty();
 
     if (!numericParamPresent && ellpsStr.empty() && datumStr.empty() &&
-        step.name == "krovak") {
+        (step.name == "krovak" || step.name == "mod_krovak")) {
         ellpsStr = "bessel";
     }
 
@@ -11059,7 +11059,8 @@ PROJStringParser::Private::processAxisSwap(Step &step,
                 throw ParsingException("Unhandled order=" + orderStr);
             }
         }
-    } else if (step.name == "krovak" && hasParamValue(step, "czech")) {
+    } else if ((step.name == "krovak" || step.name == "mod_krovak") &&
+               hasParamValue(step, "czech")) {
         axis[0] = west;
         axis[1] = south;
     }
@@ -11519,6 +11520,16 @@ PROJStringParser::Private::buildProjectedCRS(int iStep,
     } else if (step.name == "krovak" && iAxisSwap < 0 &&
                hasParamValue(step, "czech") && !hasParamValue(step, "axis")) {
         mapping = getMapping(EPSG_CODE_METHOD_KROVAK);
+    } else if (step.name == "mod_krovak" &&
+               ((iAxisSwap < 0 && getParamValue(step, "axis") == "swu" &&
+                 !hasParamValue(step, "czech")) ||
+                (iAxisSwap > 0 &&
+                 getParamValue(steps_[iAxisSwap], "order") == "-2,-1" &&
+                 !hasParamValue(step, "czech")))) {
+        mapping = getMapping(EPSG_CODE_METHOD_KROVAK_MODIFIED);
+    } else if (step.name == "mod_krovak" && iAxisSwap < 0 &&
+               hasParamValue(step, "czech") && !hasParamValue(step, "axis")) {
+        mapping = getMapping(EPSG_CODE_METHOD_KROVAK_MODIFIED);
     } else if (step.name == "merc") {
         if (hasParamValue(step, "a") && hasParamValue(step, "b") &&
             getParamValue(step, "a") == getParamValue(step, "b") &&
@@ -11692,7 +11703,7 @@ PROJStringParser::Private::buildProjectedCRS(int iStep,
                 if (!paramValue->empty()) {
                     value = getAngularValue(*paramValue);
                 }
-            } else if (step.name == "krovak") {
+            } else if (step.name == "krovak" || step.name == "mod_krovak") {
                 // Keep it in sync with defaults of krovak.cpp
                 if (param->epsg_code ==
                     EPSG_CODE_PARAMETER_LATITUDE_PROJECTION_CENTRE) {
@@ -12278,7 +12289,8 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
                 continue;
             }
             foundKeys.insert(kv.key);
-            if (step.name == "krovak" && kv.key == "alpha") {
+            if ((step.name == "krovak" || step.name == "mod_krovak") &&
+                kv.key == "alpha") {
                 // We recognize it in our CRS parsing code
                 recognizedByPROJ = true;
             } else {

diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp
@@ -3925,7 +3925,9 @@ void Conversion::_exportToPROJString(
                 "y_0", parameterValueNumericAsSI(
                            EPSG_CODE_PARAMETER_NORTHING_PROJECTION_CENTRE));
         }
-    } else if (methodEPSGCode == EPSG_CODE_METHOD_KROVAK_NORTH_ORIENTED) {
+    } else if (methodEPSGCode == EPSG_CODE_METHOD_KROVAK_NORTH_ORIENTED ||
+               methodEPSGCode ==
+                   EPSG_CODE_METHOD_KROVAK_MODIFIED_NORTH_ORIENTED) {
         double colatitude =
             parameterValueNumeric(EPSG_CODE_PARAMETER_COLATITUDE_CONE_AXIS,
                                   common::UnitOfMeasure::DEGREE);
@@ -4184,7 +4186,9 @@ void Conversion::_exportToPROJString(
             if (mapping->proj_name_aux) {
                 bool addAux = true;
                 if (internal::starts_with(mapping->proj_name_aux, "axis=")) {
-                    if (mapping->epsg_code == EPSG_CODE_METHOD_KROVAK) {
+                    if (mapping->epsg_code == EPSG_CODE_METHOD_KROVAK ||
+                        mapping->epsg_code ==
+                            EPSG_CODE_METHOD_KROVAK_MODIFIED) {
                         auto projCRS = dynamic_cast<const crs::ProjectedCRS *>(
                             l_targetCRS.get());
                         if (projCRS) {

diff --git a/src/iso19111/operation/parammappings.cpp b/src/iso19111/operation/parammappings.cpp
@@ -760,6 +760,13 @@ static const MethodMapping projectionMethodMappings[] = {
     {EPSG_NAME_METHOD_KROVAK, EPSG_CODE_METHOD_KROVAK, "Krovak", "krovak",
      "axis=swu", krovakParameters},
 
+    {EPSG_NAME_METHOD_KROVAK_MODIFIED_NORTH_ORIENTED,
+     EPSG_CODE_METHOD_KROVAK_MODIFIED_NORTH_ORIENTED, nullptr, "mod_krovak",
+     nullptr, krovakParameters},
+
+    {EPSG_NAME_METHOD_KROVAK_MODIFIED, EPSG_CODE_METHOD_KROVAK_MODIFIED,
+     nullptr, "mod_krovak", "axis=swu", krovakParameters},
+
     {EPSG_NAME_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA,
      EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA,
      "Lambert_Azimuthal_Equal_Area", "laea", nullptr, paramsLaea},

diff --git a/src/pj_list.h b/src/pj_list.h
@@ -104,6 +104,7 @@ PROJ_HEAD(merc, "Mercator")
 PROJ_HEAD(mil_os, "Miller Oblated Stereographic")
 PROJ_HEAD(mill, "Miller Cylindrical")
 PROJ_HEAD(misrsom, "Space oblique for MISR")
+PROJ_HEAD(mod_krovak, "Modified Krovak")
 PROJ_HEAD(moll, "Mollweide")
 PROJ_HEAD(molobadekas,
           "Molodensky-Badekas transform") /* implemented in PJ_helmert.c */

diff --git a/src/proj_constants.h b/src/proj_constants.h
@@ -173,6 +173,13 @@
 #define EPSG_NAME_METHOD_KROVAK "Krovak"
 #define EPSG_CODE_METHOD_KROVAK 9819
 
+#define EPSG_NAME_METHOD_KROVAK_MODIFIED "Krovak Modified"
+#define EPSG_CODE_METHOD_KROVAK_MODIFIED 1042
+
+#define EPSG_NAME_METHOD_KROVAK_MODIFIED_NORTH_ORIENTED                        \
+    "Krovak Modified (North Orientated)"
+#define EPSG_CODE_METHOD_KROVAK_MODIFIED_NORTH_ORIENTED 1043
+
 #define EPSG_NAME_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA                          \
     "Lambert Azimuthal Equal Area"
 #define EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA 9820
-Original file line number
+Diff line change
@@ Expand Up / @@ -93,6 +93,7 @@ Projections map the spherical 3D space to a flat 2D space. @@
        mil_os
        mill
        misrsom
+       mod_krovak
        moll
        murd1
        murd2
@@ Expand Down @@