Skip to content

Commit

Permalink
Merge pull request #2357 from rouault/fix_2356
Browse files Browse the repository at this point in the history 
Adjust createBoundCRSToWGS84IfPossible() and operation filtering (for POSGAR 2007 to WGS84 issues) (fixes #2356)
  • Loading branch information
@rouault
rouault committed Sep 18, 2020
2 parents 1e11bf6 + 75074ce commit 81aaa47
Show file tree
Hide file tree
Showing 7 changed files with 121 additions and 154 deletions.
133 changes: 49 additions & 84 deletions src/iso19111/coordinateoperation.cpp
View file
Expand Up @@ -11380,8 +11380,31 @@ struct SortFunction {
return false;
}

// Arbitrary final criterion
return a_name < b_name;
// Arbitrary final criterion. We actually return the greater element
// first, so that "Amersfoort to WGS 84 (4)" is presented before
// "Amersfoort to WGS 84 (3)", which is probably a better guess.

// Except for French NTF (Paris) to NTF, where the (1) conversion
// should be preferred because in the remarks of (2), it is mentionned
// OGP prefers value from IGN Paris (code 1467)...
if (a_name.find("NTF (Paris) to NTF (1)") != std::string::npos &&
b_name.find("NTF (Paris) to NTF (2)") != std::string::npos) {
return true;
}
if (a_name.find("NTF (Paris) to NTF (2)") != std::string::npos &&
b_name.find("NTF (Paris) to NTF (1)") != std::string::npos) {
return false;
}
if (a_name.find("NTF (Paris) to RGF93 (1)") != std::string::npos &&
b_name.find("NTF (Paris) to RGF93 (2)") != std::string::npos) {
return true;
}
if (a_name.find("NTF (Paris) to RGF93 (2)") != std::string::npos &&
b_name.find("NTF (Paris) to RGF93 (1)") != std::string::npos) {
return false;
}

return a_name > b_name;
}

bool operator()(const CoordinateOperationNNPtr &a,
Expand All @@ -11407,6 +11430,23 @@ static size_t getStepCount(const CoordinateOperationNNPtr &op) {

// ---------------------------------------------------------------------------

// Return number of steps that are transformations (and not conversions)
static size_t getTransformationStepCount(const CoordinateOperationNNPtr &op) {
auto concat = dynamic_cast<const ConcatenatedOperation *>(op.get());
size_t stepCount = 1;
if (concat) {
stepCount = 0;
for (const auto &subOp : concat->operations()) {
if (dynamic_cast<const Conversion *>(subOp.get()) == nullptr) {
stepCount++;
}
}
}
return stepCount;
}

// ---------------------------------------------------------------------------

static bool isNullTransformation(const std::string &name) {
if (name.find(" + ") != std::string::npos)
return false;
Expand Down Expand Up @@ -11442,8 +11482,7 @@ struct FilterResults {
// results
// ...
removeSyntheticNullTransforms();
if (context->getDiscardSuperseded())
removeUninterestingOps();
removeUninterestingOps();
removeDuplicateOps();
removeSyntheticNullTransforms();
return *this;
Expand Down Expand Up @@ -11758,50 +11797,23 @@ struct FilterResults {
void removeUninterestingOps() {

// Eliminate operations that bring nothing, ie for a given area of use,
// do not keep operations that have greater accuracy. Actually we must
// be a bit more subtle than that, and take into account grid
// availability
// do not keep operations that have similar or worse accuracy, but
// involve more (non conversion) steps
std::vector<CoordinateOperationNNPtr> resTemp;
metadata::ExtentPtr lastExtent;
double lastAccuracy = -1;
bool lastHasGrids = false;
bool lastGridsAvailable = true;
std::set<std::set<std::string>> setOfSetOfGrids;
size_t lastStepCount = 0;
CoordinateOperationPtr lastOp;

bool first = true;
const auto gridAvailabilityUse = context->getGridAvailabilityUse();
for (const auto &op : res) {
const auto curAccuracy = getAccuracy(op);
bool dummy = false;
const auto curExtent = getExtent(op, true, dummy);
bool curHasGrids = false;
bool curGridsAvailable = true;
std::set<std::string> curSetOfGrids;

const auto curStepCount = getStepCount(op);

if (context->getAuthorityFactory()) {
const auto gridsNeeded = op->gridsNeeded(
context->getAuthorityFactory()->databaseContext(),
gridAvailabilityUse ==
CoordinateOperationContext::GridAvailabilityUse::
KNOWN_AVAILABLE);
for (const auto &gridDesc : gridsNeeded) {
curHasGrids = true;
curSetOfGrids.insert(gridDesc.shortName);
if (!gridDesc.available) {
curGridsAvailable = false;
}
}
}
const auto curStepCount = getTransformationStepCount(op);

if (first) {
resTemp.emplace_back(op);

lastHasGrids = curHasGrids;
lastGridsAvailable = curGridsAvailable;
first = false;
} else {
if (lastOp->_isEquivalentTo(op.get())) {
Expand All @@ -11812,66 +11824,19 @@ struct FilterResults {
(curExtent && lastExtent &&
curExtent->contains(NN_NO_CHECK(lastExtent)) &&
lastExtent->contains(NN_NO_CHECK(curExtent))));
if (((curAccuracy > lastAccuracy && lastAccuracy >= 0) ||
if (((curAccuracy >= lastAccuracy && lastAccuracy >= 0) ||
(curAccuracy < 0 && lastAccuracy >= 0)) &&
sameExtent) {
// If that set of grids has always been used for that
// extent,
// no need to add them again
if (setOfSetOfGrids.find(curSetOfGrids) !=
setOfSetOfGrids.end()) {
continue;
}

const bool sameNameOrEmptyName =
((!curExtent && !lastExtent) ||
(curExtent && lastExtent &&
!curExtent->description()->empty() &&
*(curExtent->description()) ==
*(lastExtent->description())));

// If we have already found a operation without grids for
// that extent, no need to add any lower accuracy operation
if (!lastHasGrids && sameNameOrEmptyName) {
continue;
}
// If we had only operations involving grids, but one
// past operation had available grids, no need to add
// the new one.
if (curHasGrids && curGridsAvailable &&
lastGridsAvailable) {
continue;
}
} else if (curAccuracy == lastAccuracy && sameExtent) {
if (curStepCount > lastStepCount) {
continue;
}
sameExtent && curStepCount > lastStepCount) {
continue;
}

resTemp.emplace_back(op);

if (sameExtent) {
if (!curHasGrids) {
lastHasGrids = false;
}
if (curGridsAvailable) {
lastGridsAvailable = true;
}
} else {
setOfSetOfGrids.clear();

lastHasGrids = curHasGrids;
lastGridsAvailable = curGridsAvailable;
}
}

lastOp = op.as_nullable();
lastStepCount = curStepCount;
lastExtent = curExtent;
lastAccuracy = curAccuracy;
if (!curSetOfGrids.empty()) {
setOfSetOfGrids.insert(curSetOfGrids);
}
}
res = std::move(resTemp);
}
Expand Down
25 changes: 20 additions & 5 deletions src/iso19111/crs.cpp
View file
Expand Up @@ -378,7 +378,8 @@ VerticalCRSPtr CRS::extractVerticalCRS() const {
* a +towgs84 parameter or a WKT1:GDAL string with a TOWGS node.
*
* This method will fetch the GeographicCRS of this CRS and find a
* transformation to EPSG:4326 using the domain of the validity of the main CRS.
* transformation to EPSG:4326 using the domain of the validity of the main CRS,
* and there's only one Helmert transformation.
*
* @return a CRS.
*/
Expand Down Expand Up @@ -456,6 +457,7 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
auto list =
operation::CoordinateOperationFactory::create()
->createOperations(NN_NO_CHECK(geodCRS), hubCRS, ctxt);
CRSPtr candidateBoundCRS;
for (const auto &op : list) {
auto transf =
util::nn_dynamic_pointer_cast<operation::Transformation>(
Expand All @@ -466,8 +468,13 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
} catch (const std::exception &) {
continue;
}
return util::nn_static_pointer_cast<CRS>(BoundCRS::create(
thisAsCRS, hubCRS, NN_NO_CHECK(transf)));
if (candidateBoundCRS) {
candidateBoundCRS = nullptr;
break;
}
candidateBoundCRS =
BoundCRS::create(thisAsCRS, hubCRS, NN_NO_CHECK(transf))
.as_nullable();
} else {
auto concatenated =
dynamic_cast<const operation::ConcatenatedOperation *>(
Expand Down Expand Up @@ -499,15 +506,23 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
} catch (const std::exception &) {
continue;
}
return util::nn_static_pointer_cast<CRS>(
if (candidateBoundCRS) {
candidateBoundCRS = nullptr;
break;
}
candidateBoundCRS =
BoundCRS::create(thisAsCRS, hubCRS,
NN_NO_CHECK(transf)));
NN_NO_CHECK(transf))
.as_nullable();
}
}
}
}
}
}
if (candidateBoundCRS) {
return NN_NO_CHECK(candidateBoundCRS);
}
} catch (const std::exception &) {
}
}
Expand Down
4 changes: 2 additions & 2 deletions test/cli/proj_outIGNF.dist
View file
Expand Up @@ -8,7 +8,7 @@
311552.5340 1906457.4840 0.0000 358799.172 6342652.486 0.000
960488.4138 1910172.8812 0.0000 1007068.686 6340907.237 0.000
600000.0000 1699510.8340 0.0000 645204.279 6133556.746 0.000
1203792.5981 626873.17210 0.0000 * * inf
1203792.5981 626873.17210 0.0000 1238837.253 5057451.037 0.000
+init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX
600000.0000 2600545.4523 0.0000 179040.148 5610495.275 0.000
135638.3592 2418760.4094 0.0000 -303729.363 5410118.356 0.000
Expand All @@ -17,7 +17,7 @@
311552.5340 1906457.4840 0.0000 -96825.465 4909184.136 0.000
960488.4138 1910172.8812 0.0000 523880.019 4909191.141 0.000
600000.0000 1699510.8340 0.0000 179047.633 4708817.007 0.000
1203792.5981 626873.17210 0.0000 * * inf
1203792.5981 626873.17210 0.0000 658259.467 3623786.764 0.000
+init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX
2d20'11.4239243" 50d23'59.7718445" 0.0 179040.151 5610495.281 0.000
-3d57'49.4051448" 48d35'59.7121716" 0.0 -303729.365 5410118.352 0.000
Expand Down
10 changes: 7 additions & 3 deletions test/cli/testprojinfo_out.dist
View file
Expand Up @@ -1044,25 +1044,29 @@ Testing -s +proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs -t EPSG:432
+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=@foo.gtx +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

Testing -s "GDA94" -t "WGS 84 (G1762)" --spatial-test intersects --summary. Should include transformations through ITRF2008 and GDA2020
Candidate operations found: 6
Candidate operations found: 7
unknown id, GDA94 to GDA2020 (1) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.21 m, Australia - GDA
unknown id, GDA94 to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 5 m, Australia - GDA
unknown id, Conversion from GDA94 (geog2D) to GDA94 (geocentric) + Inverse of ITRF2008 to GDA94 (1) + Inverse of WGS 84 (G1762) to ITRF2008 (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.04 m, Australia - onshore and EEZ
unknown id, GDA94 to GDA2020 (2) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore
unknown id, GDA94 to GDA2020 (3) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore
unknown id, GDA94 to GDA2020 (2) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore
unknown id, GDA94 to GDA2020 (5) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Cocos (Keeling) Islands - onshore
unknown id, GDA94 to GDA2020 (4) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Christmas Island - onshore

Testing -s "AGD66" -t "WGS 84 (G1762)" --spatial-test intersects --summary. Should include a transformation through GDA2020
Candidate operations found: 11
Candidate operations found: 14
unknown id, AGD66 to WGS 84 (18) + WGS 84 to WGS 84 (G1762), 5 m, Australia - offshore
unknown id, AGD66 to WGS 84 (16) + WGS 84 to WGS 84 (G1762), 7 m, Australia - onshore
unknown id, AGD66 to WGS 84 (20) + WGS 84 to WGS 84 (G1762), 11 m, Australia - onshore
unknown id, AGD66 to WGS 84 (15) + WGS 84 to WGS 84 (G1762), 3 m, Australia - Northern Territory
unknown id, AGD66 to WGS 84 (13) + WGS 84 to WGS 84 (G1762), 3 m, Australia - New South Wales and Victoria
unknown id, AGD66 to WGS 84 (21) + WGS 84 to WGS 84 (G1762), 7 m, Papua New Guinea - mainland onshore
unknown id, AGD66 to WGS 84 (14) + WGS 84 to WGS 84 (G1762), 3 m, Australia - Tasmania
unknown id, AGD66 to WGS 84 (19) + WGS 84 to WGS 84 (G1762), 4 m, Papua New Guinea - PFTB
unknown id, AGD66 to WGS 84 (22) + WGS 84 to WGS 84 (G1762), 6 m, Papua New Guinea - PFTB
unknown id, AGD66 to WGS 84 (23) + WGS 84 to WGS 84 (G1762), 6 m, Papua New Guinea - North Fly
unknown id, AGD66 to GDA2020 (1) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - Australian Capital Territory
unknown id, AGD66 to WGS 84 (12) + WGS 84 to WGS 84 (G1762), 3 m, Australia - Australian Capital Territory
unknown id, AGD66 to WGS 84 (17) + WGS 84 to WGS 84 (G1762), 3 m, Australia - onshore
unknown id, Ballpark geographic offset from AGD66 to WGS 84 (G1762), unknown accuracy, World, has ballpark transformation

Expand Down
15 changes: 6 additions & 9 deletions test/unit/test_c_api.cpp
View file
Expand Up @@ -540,7 +540,7 @@ TEST_F(CApi, proj_as_proj_string_approx_tmerc_option_yes) {
// ---------------------------------------------------------------------------

TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) {
auto crs = proj_create_from_database(m_ctxt, "EPSG", "3844",
auto crs = proj_create_from_database(m_ctxt, "EPSG", "4807",
PJ_CATEGORY_CRS, false, nullptr);
ObjectKeeper keeper(crs);
ASSERT_NE(crs, nullptr);
Expand All @@ -552,10 +552,8 @@ TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) {
auto proj_4 = proj_as_proj_string(m_ctxt, res, PJ_PROJ_4, nullptr);
ASSERT_NE(proj_4, nullptr);
EXPECT_EQ(std::string(proj_4),
"+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 "
"+y_0=500000 +ellps=krass "
"+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69 "
"+units=m +no_defs +type=crs");
"+proj=longlat +ellps=clrk80ign +pm=paris "
"+towgs84=-168,-60,320,0,0,0,0 +no_defs +type=crs");

auto base_crs = proj_get_source_crs(m_ctxt, res);
ObjectKeeper keeper_base_crs(base_crs);
Expand All @@ -572,8 +570,7 @@ TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) {
std::vector<double> values(7, 0);
EXPECT_TRUE(proj_coordoperation_get_towgs84_values(m_ctxt, transf,
values.data(), 7, true));
auto expected = std::vector<double>{2.329, -147.042, -92.08, -0.309,
0.325, 0.497, 5.69};
auto expected = std::vector<double>{-168, -60, 320, 0, 0, 0, 0};
EXPECT_EQ(values, expected);

auto res2 = proj_crs_create_bound_crs(m_ctxt, base_crs, hub_crs, transf);
Expand Down Expand Up @@ -1485,7 +1482,7 @@ TEST_F(CApi, proj_create_operations_discard_superseded) {
ASSERT_NE(res, nullptr);
ObjListKeeper keeper_res(res);

EXPECT_EQ(proj_list_get_count(res), 2);
EXPECT_EQ(proj_list_get_count(res), 4);
}

// ---------------------------------------------------------------------------
Expand Down Expand Up @@ -1594,7 +1591,7 @@ TEST_F(CApi, proj_create_operations_with_pivot) {
auto res = proj_create_operations(m_ctxt, source_crs, target_crs, ctxt);
ASSERT_NE(res, nullptr);
ObjListKeeper keeper_res(res);
EXPECT_EQ(proj_list_get_count(res), 7);
EXPECT_EQ(proj_list_get_count(res), 8);
auto op = proj_list_get(m_ctxt, res, 0);
ASSERT_NE(op, nullptr);
ObjectKeeper keeper_op(op);
Expand Down

0 comments on commit 81aaa47

Please sign in to comment.