Skip to content

Commit

Permalink
Add toWGS84AutocorrectWrongValues() method and use it in PROJ.4 and W…
Browse files Browse the repository at this point in the history 
…KT1 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
  • Loading branch information
@rouault
rouault committed Jun 15, 2024
1 parent e0e1482 commit cceff3e
Show file tree
Hide file tree
Showing 6 changed files with 298 additions and 0 deletions.
5 changes: 5 additions & 0 deletions include/proj/io.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -960,6 +960,11 @@ class PROJ_GCC_DLL DatabaseContext {
PROJ_DLL std::vector<std::string>
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:
Expand Down
1 change: 1 addition & 0 deletions scripts/reference_exported_symbols.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::shared_ptr<osgeo::proj::common::IdentifiedObject> > 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&)
Expand Down
92 changes: 92 additions & 0 deletions src/iso19111/factory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

// ---------------------------------------------------------------------------
Expand Down
29 changes: 29 additions & 0 deletions src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand Down Expand Up @@ -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);
}

Expand Down
102 changes: 102 additions & 0 deletions test/unit/test_factory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
69 changes: 69 additions & 0 deletions test/unit/test_io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<BoundCRS>(obj);
ASSERT_TRUE(crs != nullptr);

auto params = crs->transformation()->getTOWGS84Parameters();
auto expected = std::vector<double>{-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"
Expand Down Expand Up @@ -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<BoundCRS>(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");
Expand Down

0 comments on commit cceff3e

Please sign in to comment.