From aeb024f2df5fb1c4d83959b6b9edcb62e8385226 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 17 Sep 2020 18:34:06 +0200 Subject: Adjust createBoundCRSToWGS84IfPossible() and operation filtering (for POSGAR 2007 to WGS84 issues) (fixes #2356) - We make createBoundCRSToWGS84IfPossible() more restrictive. If there are more than one Helmert transformation from the CRS to WGS 84 covering the area of use of the CRS, we do not create a BoundCRS / +towgs84 - In createOperations() filtering, we are less aggressive in discarding operations that have the same area of use but worse accuracy. We do it only if they involve more transformation steps. We now get: ``` $ projinfo EPSG:5340 -o PROJ PROJ.4 string: +proj=longlat +ellps=GRS80 +no_defs +type=crs $ projinfo -s EPSG:5340 -t EPSG:4326 --spatial-test intersects --summary Candidate operations found: 2 EPSG:9264, POSGAR 2007 to WGS 84 (2), 0.5 m, Argentina EPSG:5351, POSGAR 2007 to WGS 84 (1), 1.0 m, Argentina ``` --- src/iso19111/coordinateoperation.cpp | 106 ++++++++--------------------------- src/iso19111/crs.cpp | 25 +++++++-- 2 files changed, 44 insertions(+), 87 deletions(-) (limited to 'src') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 10b1baf4..0e61d2ee 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -11407,6 +11407,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(op.get()); + size_t stepCount = 1; + if (concat) { + stepCount = 0; + for (const auto &subOp : concat->operations()) { + if (dynamic_cast(subOp.get()) == nullptr) { + stepCount++; + } + } + } + return stepCount; +} + +// --------------------------------------------------------------------------- + static bool isNullTransformation(const std::string &name) { if (name.find(" + ") != std::string::npos) return false; @@ -11442,8 +11459,7 @@ struct FilterResults { // results // ... removeSyntheticNullTransforms(); - if (context->getDiscardSuperseded()) - removeUninterestingOps(); + removeUninterestingOps(); removeDuplicateOps(); removeSyntheticNullTransforms(); return *this; @@ -11758,50 +11774,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 resTemp; metadata::ExtentPtr lastExtent; double lastAccuracy = -1; - bool lastHasGrids = false; - bool lastGridsAvailable = true; - std::set> 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 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())) { @@ -11812,66 +11801,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); } diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 88420c8a..e96b3cc9 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -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. */ @@ -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( @@ -466,8 +468,13 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible( } catch (const std::exception &) { continue; } - return util::nn_static_pointer_cast(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( @@ -499,15 +506,23 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible( } catch (const std::exception &) { continue; } - return util::nn_static_pointer_cast( + 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 &) { } } -- cgit v1.2.3 From 75074ce863e36299ec39b39decdd435eed15ad63 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 17 Sep 2020 19:46:07 +0200 Subject: createOperations(): tune sorting of transformations so that the ones with greater 'version numbers' are prefered over other ones, when all other comparison criteria are equal. Helps with Amersfoort RD New to EPSG:4326 --- src/iso19111/coordinateoperation.cpp | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 0e61d2ee..65ef2b77 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -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, -- cgit v1.2.3