Skip to content

Commit

Permalink
Merge pull request #4172 from rouault/toWGS84AutocorrectWrongValues
Browse files Browse the repository at this point in the history 
Add toWGS84AutocorrectWrongValues() method and use it in PROJ.4 and WKT1 CRS import
  • Loading branch information
@rouault
rouault authored Jun 18, 2024
2 parents d978fdc + cceff3e commit 68c678b
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 @@ -2479,6 +2479,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 @@ -11466,6 +11475,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 @@ -5280,6 +5280,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 @@ -10963,6 +11007,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 68c678b

Please sign in to comment.