Merge pull request #2985 from rouault/wkt1_hotine_without_rectified_grid_angle

WKT1 import: correctly deal with missing rectified_grid_angle parameter

3 files changed, 178 insertions, 7 deletions

diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index f8a4672a..454800dd 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3894,13 +3894,12 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard(
             if (mapping &&
                 mapping->epsg_code == EPSG_CODE_METHOD_MERCATOR_VARIANT_B &&
                 ci_equal(parameterName, "latitude_of_origin")) {
-                for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) {
-                    if (mapping->params[idx]->epsg_code ==
-                        EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN) {
-                        foundParameters[idx] = true;
-                        break;
-                    }
-                }
+                // Some illegal formulations of Mercator_2SP have a unexpected
+                // latitude_of_origin parameter. We accept it on import, but
+                // do not accept it when exporting to PROJ string, unless it is
+                // zero.
+                // No need to try to update foundParameters[] as this is a
+                // unexpected one.
                 parameterName = EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN;
                 propertiesParameter.set(
                     Identifier::CODE_KEY,
@@ -3969,6 +3968,48 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard(
         }
     }
 
+    if (mapping && (mapping->epsg_code ==
+                        EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_A ||
+                    mapping->epsg_code ==
+                        EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_B)) {
+        // Special case when importing some GDAL WKT of Hotine Oblique Mercator
+        // that have a Azimuth parameter but lacks the Rectified Grid Angle.
+        // We have code in the exportToPROJString() to deal with that situation,
+        // but also adds the rectified grid angle from the azimuth on import.
+        bool foundAngleRecifiedToSkewGrid = false;
+        bool foundAzimuth = false;
+        for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) {
+            if (foundParameters[idx] &&
+                mapping->params[idx]->epsg_code ==
+                    EPSG_CODE_PARAMETER_ANGLE_RECTIFIED_TO_SKEW_GRID) {
+                foundAngleRecifiedToSkewGrid = true;
+            } else if (foundParameters[idx] &&
+                       mapping->params[idx]->epsg_code ==
+                           EPSG_CODE_PARAMETER_AZIMUTH_INITIAL_LINE) {
+                foundAzimuth = true;
+            }
+        }
+        if (!foundAngleRecifiedToSkewGrid && foundAzimuth) {
+            for (size_t idx = 0; idx < parameters.size(); ++idx) {
+                if (parameters[idx]->getEPSGCode() ==
+                    EPSG_CODE_PARAMETER_AZIMUTH_INITIAL_LINE) {
+                    PropertyMap propertiesParameter;
+                    propertiesParameter.set(
+                        Identifier::CODE_KEY,
+                        EPSG_CODE_PARAMETER_ANGLE_RECTIFIED_TO_SKEW_GRID);
+                    propertiesParameter.set(Identifier::CODESPACE_KEY,
+                                            Identifier::EPSG);
+                    propertiesParameter.set(
+                        IdentifiedObject::NAME_KEY,
+                        EPSG_NAME_PARAMETER_ANGLE_RECTIFIED_TO_SKEW_GRID);
+                    parameters.push_back(
+                        OperationParameter::create(propertiesParameter));
+                    values.push_back(values[idx]);
+                }
+            }
+        }
+    }
+
     return Conversion::create(
                PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"),
                propertiesMethod, parameters, values)
diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp
index 9927892d..77e703f3 100644
--- a/src/iso19111/operation/conversion.cpp
+++ b/src/iso19111/operation/conversion.cpp
@@ -4005,6 +4005,16 @@ void Conversion::_exportToPROJString(
                         EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN) {
                         valueConverted = 1.0;
                     }
+                    if ((mapping->epsg_code ==
+                             EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_A ||
+                         mapping->epsg_code ==
+                             EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_B) &&
+                        param->epsg_code ==
+                            EPSG_CODE_PARAMETER_ANGLE_RECTIFIED_TO_SKEW_GRID) {
+                        // Do not use 0 as the default value for +gamma of
+                        // proj=omerc
+                        continue;
+                    }
                 } else if (param->unit_type ==
                            common::UnitOfMeasure::Type::ANGULAR) {
                     valueConverted =
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 196552d4..5ffa461b 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -1770,6 +1770,126 @@ TEST(wkt_parse, wkt1_Spherical_Cross_Track_Height) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt1_hotine_oblique_mercator_without_rectified_grid_angle) {
+    auto wkt = "PROJCS[\"NAD_1983_Michigan_GeoRef_Meters\","
+               "GEOGCS[\"NAD83(1986)\","
+               "DATUM[\"North_American_Datum_1983\","
+               "SPHEROID[\"GRS_1980\",6378137,298.257222101]],"
+               "PRIMEM[\"Greenwich\",0],"
+               "UNIT[\"Degree\",0.017453292519943295]],"
+               "PROJECTION[\"Hotine_Oblique_Mercator\"],"
+               "PARAMETER[\"false_easting\",2546731.496],"
+               "PARAMETER[\"false_northing\",-4354009.816],"
+               "PARAMETER[\"latitude_of_center\",45.30916666666666],"
+               "PARAMETER[\"longitude_of_center\",-86],"
+               "PARAMETER[\"azimuth\",-22.74444],"
+               "PARAMETER[\"scale_factor\",0.9996],"
+               "UNIT[\"Meter\",1]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    // Check that we have added automatically rectified_grid_angle
+    auto got_wkt = crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+    EXPECT_TRUE(got_wkt.find("PARAMETER[\"rectified_grid_angle\",-22.74444]") !=
+                std::string::npos)
+        << got_wkt;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_hotine_oblique_mercator_with_rectified_grid_angle) {
+    auto wkt = "PROJCS[\"NAD_1983_Michigan_GeoRef_Meters\","
+               "GEOGCS[\"NAD83(1986)\","
+               "DATUM[\"North_American_Datum_1983\","
+               "SPHEROID[\"GRS_1980\",6378137,298.257222101]],"
+               "PRIMEM[\"Greenwich\",0],"
+               "UNIT[\"Degree\",0.017453292519943295]],"
+               "PROJECTION[\"Hotine_Oblique_Mercator\"],"
+               "PARAMETER[\"false_easting\",2546731.496],"
+               "PARAMETER[\"false_northing\",-4354009.816],"
+               "PARAMETER[\"latitude_of_center\",45.30916666666666],"
+               "PARAMETER[\"longitude_of_center\",-86],"
+               "PARAMETER[\"azimuth\",-22.74444],"
+               "PARAMETER[\"rectified_grid_angle\",-23],"
+               "PARAMETER[\"scale_factor\",0.9996],"
+               "UNIT[\"Meter\",1]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    // Check that we have not overriden rectified_grid_angle
+    auto got_wkt = crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+    EXPECT_TRUE(got_wkt.find("PARAMETER[\"rectified_grid_angle\",-23]") !=
+                std::string::npos)
+        << got_wkt;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(proj_export, wkt2_hotine_oblique_mercator_without_rectified_grid_angle) {
+    auto wkt = "PROJCRS[\"NAD_1983_Michigan_GeoRef_Meters\",\n"
+               "    BASEGEOGCRS[\"NAD83(1986)\",\n"
+               "        DATUM[\"North American Datum 1983\",\n"
+               "            ELLIPSOID[\"GRS_1980\",6378137,298.257222101,\n"
+               "                LENGTHUNIT[\"metre\",1]],\n"
+               "            ID[\"EPSG\",6269]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n"
+               "    CONVERSION[\"unnamed\",\n"
+               "        METHOD[\"Hotine Oblique Mercator (variant A)\",\n"
+               "            ID[\"EPSG\",9812]],\n"
+               "        PARAMETER[\"False easting\",2546731.496,\n"
+               "            LENGTHUNIT[\"Meter\",1],\n"
+               "            ID[\"EPSG\",8806]],\n"
+               "        PARAMETER[\"False northing\",-4354009.816,\n"
+               "            LENGTHUNIT[\"Meter\",1],\n"
+               "            ID[\"EPSG\",8807]],\n"
+               "        PARAMETER[\"Latitude of projection centre\","
+               "                  45.3091666666667,\n"
+               "            ANGLEUNIT[\"Degree\",0.0174532925199433],\n"
+               "            ID[\"EPSG\",8811]],\n"
+               "        PARAMETER[\"Longitude of projection centre\",-86,\n"
+               "            ANGLEUNIT[\"Degree\",0.0174532925199433],\n"
+               "            ID[\"EPSG\",8812]],\n"
+               "        PARAMETER[\"Azimuth of initial line\",-22.74444,\n"
+               "            ANGLEUNIT[\"Degree\",0.0174532925199433],\n"
+               "            ID[\"EPSG\",8813]],\n"
+               "        PARAMETER[\"Scale factor on initial line\",0.9996,\n"
+               "            SCALEUNIT[\"unity\",1],\n"
+               "            ID[\"EPSG\",8815]]],\n"
+               "    CS[Cartesian,2],\n"
+               "        AXIS[\"(E)\",east,\n"
+               "            ORDER[1],\n"
+               "            LENGTHUNIT[\"Meter\",1]],\n"
+               "        AXIS[\"(N)\",north,\n"
+               "            ORDER[2],\n"
+               "            LENGTHUNIT[\"Meter\",1]]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    // We don't do any particular handling of missing Angle from Rectified
+    // to Skew Grid on import, but on export to PROJ string,
+    // check that we don't add a dummy gamma value.
+    auto expectedPROJString = "+proj=omerc +no_uoff +lat_0=45.3091666666667 "
+                              "+lonc=-86 +alpha=-22.74444 "
+                              "+k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 "
+                              "+datum=NAD83 +units=m +no_defs +type=crs";
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        expectedPROJString);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, wkt2_projected) {
     auto wkt = "PROJCRS[\"WGS 84 / UTM zone 31N\",\n"
                "    BASEGEODCRS[\"WGS 84\",\n"