From a89b135465ba62201b853d468c7c98422b2cd38c Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 15 Jun 2024 17:44:35 +0200 Subject: [PATCH] Add toWGS84AutocorrectWrongValues() method and use it in PROJ.4 and WKT1 CRS import This adds an heuristics to detect wrong sign in rotation terms of TOWGS84 clause wrongly imported from a Coordinate Frame rotation, such as the ones currently emitted by epsg dot io. Refs #4170 --- include/proj/io.hpp | 5 ++ scripts/reference_exported_symbols.txt | 1 + src/iso19111/factory.cpp | 92 ++++++++++++++++++++++ src/iso19111/io.cpp | 29 +++++++ test/unit/test_factory.cpp | 102 +++++++++++++++++++++++++ test/unit/test_io.cpp | 69 +++++++++++++++++ 6 files changed, 298 insertions(+) diff --git a/include/proj/io.hpp b/include/proj/io.hpp index 68d9b1c8c7..aa03b7803b 100644 --- a/include/proj/io.hpp +++ b/include/proj/io.hpp @@ -960,6 +960,11 @@ class PROJ_GCC_DLL DatabaseContext { PROJ_DLL std::vector getVersionedAuthoritiesFromName(const std::string &authName); + PROJ_FOR_TEST bool + toWGS84AutocorrectWrongValues(double &tx, double &ty, double &tz, + double &rx, double &ry, double &rz, + double &scale_difference) const; + //! @endcond protected: diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt index f30edbbefa..9035283b5b 100644 --- a/scripts/reference_exported_symbols.txt +++ b/scripts/reference_exported_symbols.txt @@ -400,6 +400,7 @@ osgeo::proj::io::DatabaseContext::lookForGridInfo(std::string const&, bool, std: osgeo::proj::io::DatabaseContext::startInsertStatementsSession() osgeo::proj::io::DatabaseContext::stopInsertStatementsSession() osgeo::proj::io::DatabaseContext::suggestsCodeFor(dropbox::oxygen::nn > const&, std::string const&, bool) +osgeo::proj::io::DatabaseContext::toWGS84AutocorrectWrongValues(double&, double&, double&, double&, double&, double&, double&) const osgeo::proj::io::FactoryException::~FactoryException() osgeo::proj::io::FactoryException::FactoryException(char const*) osgeo::proj::io::FactoryException::FactoryException(osgeo::proj::io::FactoryException const&) diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index f20c5f363b..bd4acf2b65 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -3851,6 +3851,98 @@ DatabaseContext::getTransformationsForGridName( return res; } +// --------------------------------------------------------------------------- + +// Fixes wrong towgs84 values returned by epsg.io when using a Coordinate Frame +// transformation, where they neglect to reverse the sign of the rotation terms. +// Cf https://github.com/OSGeo/PROJ/issues/4170 and +// https://github.com/maptiler/epsg.io/issues/194 +// We do that only when we found a valid Coordinate Frame rotation that +// has the same numeric values (and no corresponding Position Vector +// transformation with same values, or Coordinate Frame transformation with +// opposite sign for rotation terms, both are highly unlikely) +bool DatabaseContext::toWGS84AutocorrectWrongValues( + double &tx, double &ty, double &tz, double &rx, double &ry, double &rz, + double &scale_difference) const { + if (rx == 0 && ry == 0 && rz == 0) + return false; + // 9606: Coordinate Frame rotation (geog2D domain) + // 9607: Position Vector transformation (geog2D domain) + std::string sql( + "SELECT DISTINCT method_code " + "FROM helmert_transformation_table WHERE " + "abs(tx - ?) <= 1e-8 * abs(tx) AND " + "abs(ty - ?) <= 1e-8 * abs(ty) AND " + "abs(tz - ?) <= 1e-8 * abs(tz) AND " + "abs(rx - ?) <= 1e-8 * abs(rx) AND " + "abs(ry - ?) <= 1e-8 * abs(ry) AND " + "abs(rz - ?) <= 1e-8 * abs(rz) AND " + "abs(scale_difference - ?) <= 1e-8 * abs(scale_difference) AND " + "method_auth_name = 'EPSG' AND " + "method_code IN (9606, 9607) AND " + "translation_uom_auth_name = 'EPSG' AND " + "translation_uom_code = 9001 AND " // metre + "rotation_uom_auth_name = 'EPSG' AND " + "rotation_uom_code = 9104 AND " // arc-second + "scale_difference_uom_auth_name = 'EPSG' AND " + "scale_difference_uom_code = 9202 AND " // parts per million + "deprecated = 0"); + ListOfParams params; + params.emplace_back(tx); + params.emplace_back(ty); + params.emplace_back(tz); + params.emplace_back(rx); + params.emplace_back(ry); + params.emplace_back(rz); + params.emplace_back(scale_difference); + bool bFound9606 = false; + bool bFound9607 = false; + for (const auto &row : d->run(sql, params)) { + if (row[0] == "9606") { + bFound9606 = true; + } else if (row[0] == "9607") { + bFound9607 = true; + } + } + if (bFound9607 && !bFound9606) { + params.clear(); + params.emplace_back(tx); + params.emplace_back(ty); + params.emplace_back(tz); + params.emplace_back(-rx); + params.emplace_back(-ry); + params.emplace_back(-rz); + params.emplace_back(scale_difference); + if (d->run(sql, params).empty()) { + if (d->pjCtxt()) { + pj_log(d->pjCtxt(), PJ_LOG_ERROR, + "Auto-correcting wrong sign of rotation terms of " + "TOWGS84 clause from %s,%s,%s,%s,%s,%s,%s to " + "%s,%s,%s,%s,%s,%s,%s", + internal::toString(tx).c_str(), + internal::toString(ty).c_str(), + internal::toString(tz).c_str(), + internal::toString(rx).c_str(), + internal::toString(ry).c_str(), + internal::toString(rz).c_str(), + internal::toString(scale_difference).c_str(), + internal::toString(tx).c_str(), + internal::toString(ty).c_str(), + internal::toString(tz).c_str(), + internal::toString(-rx).c_str(), + internal::toString(-ry).c_str(), + internal::toString(-rz).c_str(), + internal::toString(scale_difference).c_str()); + } + rx = -rx; + ry = -ry; + rz = -rz; + return true; + } + } + return false; +} + //! @endcond // --------------------------------------------------------------------------- diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 43e957371c..cfadadf099 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -2457,6 +2457,15 @@ GeodeticReferenceFrameNNPtr WKTParser::Private::buildGeodeticReferenceFrame( for (const auto &child : TOWGS84Children) { toWGS84Parameters_.push_back(asDouble(child)); } + + if (TOWGS84Size == 7 && dbContext_) { + dbContext_->toWGS84AutocorrectWrongValues( + toWGS84Parameters_[0], toWGS84Parameters_[1], + toWGS84Parameters_[2], toWGS84Parameters_[3], + toWGS84Parameters_[4], toWGS84Parameters_[5], + toWGS84Parameters_[6]); + } + for (size_t i = TOWGS84Size; i < 7; ++i) { toWGS84Parameters_.push_back(0.0); } @@ -11427,6 +11436,26 @@ PROJStringParser::Private::buildBoundOrCompoundCRSIfNeeded(int iStep, throw ParsingException("Non numerical value in towgs84 clause"); } } + + if (towgs84Values.size() == 7 && dbContext_) { + if (dbContext_->toWGS84AutocorrectWrongValues( + towgs84Values[0], towgs84Values[1], towgs84Values[2], + towgs84Values[3], towgs84Values[4], towgs84Values[5], + towgs84Values[6])) { + for (auto &pair : step.paramValues) { + if (ci_equal(pair.key, "towgs84")) { + pair.value.clear(); + for (int i = 0; i < 7; ++i) { + if (i > 0) + pair.value += ','; + pair.value += internal::toString(towgs84Values[i]); + } + break; + } + } + } + } + crs = BoundCRS::createFromTOWGS84(crs, towgs84Values); } diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp index 728a77da1d..769e6c4058 100644 --- a/test/unit/test_factory.cpp +++ b/test/unit/test_factory.cpp @@ -4935,4 +4935,106 @@ TEST(factory, getPointMotionOperationsFor) { // --------------------------------------------------------------------------- +TEST(factory, toWGS84AutocorrectWrongValues) { + auto ctxt = DatabaseContext::create(); + { + double tx = 1; + double ty = 2; + double tz = 3; + double rx = 0; + double ry = 0; + double rz = 0; + double scale_difference = 0; + EXPECT_FALSE(ctxt->toWGS84AutocorrectWrongValues(tx, ty, tz, rx, ry, rz, + scale_difference)); + EXPECT_EQ(tx, 1); + EXPECT_EQ(ty, 2); + EXPECT_EQ(tz, 3); + EXPECT_EQ(rx, 0); + EXPECT_EQ(ry, 0); + EXPECT_EQ(rz, 0); + EXPECT_EQ(scale_difference, 0); + } + { + // Incorrect parameters for EPSG:15929: WGS84 -> Belgian Lambert 72 + // Cf https://github.com/OSGeo/PROJ/issues/4170 + double tx = -106.8686; + double ty = 52.2978; + double tz = -103.7239; + double rx = -0.3366; + double ry = 0.457; + double rz = -1.8422; + double scale_difference = -1.2747; + EXPECT_TRUE(ctxt->toWGS84AutocorrectWrongValues(tx, ty, tz, rx, ry, rz, + scale_difference)); + EXPECT_EQ(tx, -106.8686); + EXPECT_EQ(ty, 52.2978); + EXPECT_EQ(tz, -103.7239); + EXPECT_EQ(rx, 0.3366); + EXPECT_EQ(ry, -0.457); + EXPECT_EQ(rz, 1.8422); + EXPECT_EQ(scale_difference, -1.2747); + } + { + // Almost incorrect parameters EPSG:15929: WGS84 -> Belgian Lambert 72 + double tx = -106; + double ty = 52.2978; + double tz = -103.7239; + double rx = -0.3366; + double ry = 0.457; + double rz = -1.8422; + double scale_difference = -1.2747; + EXPECT_FALSE(ctxt->toWGS84AutocorrectWrongValues(tx, ty, tz, rx, ry, rz, + scale_difference)); + EXPECT_EQ(tx, -106); + EXPECT_EQ(ty, 52.2978); + EXPECT_EQ(tz, -103.7239); + EXPECT_EQ(rx, -0.3366); + EXPECT_EQ(ry, 0.457); + EXPECT_EQ(rz, -1.8422); + EXPECT_EQ(scale_difference, -1.2747); + } + { + // Correct Position Vector transformation ('EPSG','15869','DHDN to WGS + // 84 (3)) + double tx = 612.4; + double ty = 77.0; + double tz = 440.2; + double rx = -0.054; + double ry = 0.057; + double rz = -2.797; + double scale_difference = 2.55; + EXPECT_FALSE(ctxt->toWGS84AutocorrectWrongValues(tx, ty, tz, rx, ry, rz, + scale_difference)); + EXPECT_EQ(tx, 612.4); + EXPECT_EQ(ty, 77.0); + EXPECT_EQ(tz, 440.2); + EXPECT_EQ(rx, -0.054); + EXPECT_EQ(ry, 0.057); + EXPECT_EQ(rz, -2.797); + EXPECT_EQ(scale_difference, 2.55); + } + { + // Correct parameters for EPSG:15929: WGS84 -> Belgian Lambert 72 + // (Coordinate Frame rotation) Cf + // https://github.com/OSGeo/PROJ/issues/4170 + double tx = -106.8686; + double ty = 52.2978; + double tz = -103.7239; + double rx = 0.3366; + double ry = -0.457; + double rz = 1.8422; + double scale_difference = -1.2747; + EXPECT_FALSE(ctxt->toWGS84AutocorrectWrongValues(tx, ty, tz, rx, ry, rz, + scale_difference)); + EXPECT_EQ(tx, -106.8686); + EXPECT_EQ(ty, 52.2978); + EXPECT_EQ(tz, -103.7239); + EXPECT_EQ(rx, 0.3366); + EXPECT_EQ(ry, -0.457); + EXPECT_EQ(rz, 1.8422); + EXPECT_EQ(scale_difference, -1.2747); + } +} + } // namespace diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 4ccaff1772..f9b0436585 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -5250,6 +5250,50 @@ TEST(wkt_parse, projcs_TOWGS84_7terms) { // --------------------------------------------------------------------------- +TEST(io, projcs_TOWGS84_7terms_autocorrect) { + // Auto-correct wrong sign for rotation terms + // Cf https://github.com/OSGeo/PROJ/issues/4170 + auto wkt = + "PROJCS[\"BD72 / Belgian Lambert 72\",\n" + " GEOGCS[\"BD72\",\n" + " DATUM[\"Reseau_National_Belge_1972\",\n" + " SPHEROID[\"International 1924\",6378388,297],\n" + " " + "TOWGS84[-106.8686,52.2978,-103.7239,-0.3366,0.457,-1.8422,-1.2747]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4313\"]],\n" + " PROJECTION[\"Lambert_Conformal_Conic_2SP\"],\n" + " PARAMETER[\"latitude_of_origin\",90],\n" + " PARAMETER[\"central_meridian\",4.36748666666667],\n" + " PARAMETER[\"standard_parallel_1\",51.1666672333333],\n" + " PARAMETER[\"standard_parallel_2\",49.8333339],\n" + " PARAMETER[\"false_easting\",150000.013],\n" + " PARAMETER[\"false_northing\",5400088.438],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Easting\",EAST],\n" + " AXIS[\"Northing\",NORTH],\n" + " AUTHORITY[\"EPSG\",\"31370\"]]"; + + auto dbContext = DatabaseContext::create(); + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + auto params = crs->transformation()->getTOWGS84Parameters(); + auto expected = std::vector{-106.8686, 52.2978, -103.7239, 0.3366, + -0.457, 1.8422, -1.2747}; + ASSERT_EQ(params.size(), expected.size()); + for (int i = 0; i < 7; i++) { + EXPECT_NEAR(params[i], expected[i], 1e-10); + } +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, WKT1_VERT_DATUM_EXTENSION) { auto wkt = "VERT_CS[\"EGM2008 geoid height\",\n" " VERT_DATUM[\"EGM2008 geoid\",2005,\n" @@ -10933,6 +10977,31 @@ TEST(io, projparse_longlat_towgs84_7_terms) { // --------------------------------------------------------------------------- +TEST(io, projparse_longlat_towgs84_7_terms_autocorrect) { + // Auto-correct wrong sign for rotation terms + // Cf https://github.com/OSGeo/PROJ/issues/4170 + auto dbContext = DatabaseContext::create(); + auto obj = createFromUserInput( + "+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 " + "+lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl " + "+towgs84=-106.8686,52.2978,-103.7239,-0.3366,0.457,-1.8422,-1.2747 " + "+units=m +no_defs +type=crs", + dbContext, true); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + "+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 " + "+lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl " + "+towgs84=-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747 " + "+units=m +no_defs +type=crs"); +} + +// --------------------------------------------------------------------------- + TEST(io, projparse_longlat_nadgrids) { auto obj = PROJStringParser().createFromPROJString( "+proj=longlat +ellps=GRS80 +nadgrids=foo.gsb +type=crs");