From cceff3e49e5d45ff2ed6d446b254a0735d30e2f0 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, aud autocorrect, 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");