Merge pull request #3160 from OSGeo/backport-3158-to-9.0

[Backport 9.0] createBoundCRSToWGS84IfPossible(): improve selection logic to generate +towgs84= taking into account extent
OSGeo · Apr 8, 2022 · 187108b · 187108b
2 parents 39397ed + 1bdeb23
commit 187108b
Show file tree

Hide file tree

Showing 3 changed files with 63 additions and 12 deletions.
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
@@ -574,6 +574,8 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
                 operation::CoordinateOperationFactory::create()
                     ->createOperations(NN_NO_CHECK(geodCRS), hubCRS, ctxt);
             CRSPtr candidateBoundCRS;
+            int candidateCount = 0;
+            bool candidateHasExactlyMachingExtent = false;
             for (const auto &op : list) {
                 auto transf =
                     util::nn_dynamic_pointer_cast<operation::Transformation>(
@@ -584,13 +586,31 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
                     } catch (const std::exception &) {
                         continue;
                     }
+                    bool unused = false;
+                    auto opExtent =
+                        getExtent(NN_NO_CHECK(transf), false, unused);
+                    const bool exactlyMatchingExtent =
+                        opExtent && extentResolved &&
+                        opExtent->contains(NN_NO_CHECK(extentResolved)) &&
+                        extentResolved->contains(NN_NO_CHECK(opExtent));
                     if (candidateBoundCRS) {
-                        candidateBoundCRS = nullptr;
-                        break;
+                        if (exactlyMatchingExtent &&
+                            !candidateHasExactlyMachingExtent) {
+                            candidateBoundCRS = nullptr;
+                        } else if (exactlyMatchingExtent ==
+                                   candidateHasExactlyMachingExtent) {
+                            candidateCount++;
+                        }
+                    }
+                    if (candidateBoundCRS == nullptr) {
+                        candidateCount = 1;
+                        candidateHasExactlyMachingExtent =
+                            exactlyMatchingExtent;
+                        candidateBoundCRS =
+                            BoundCRS::create(thisAsCRS, hubCRS,
+                                             NN_NO_CHECK(transf))
+                                .as_nullable();
                     }
-                    candidateBoundCRS =
-                        BoundCRS::create(thisAsCRS, hubCRS, NN_NO_CHECK(transf))
-                            .as_nullable();
                 } else {
                     auto concatenated =
                         dynamic_cast<const operation::ConcatenatedOperation *>(
@@ -621,21 +641,42 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
                                     } catch (const std::exception &) {
                                         continue;
                                     }
+                                    bool unused = false;
+                                    auto opExtent = getExtent(
+                                        NN_NO_CHECK(transf), false, unused);
+                                    const bool exactlyMatchingExtent =
+                                        opExtent && extentResolved &&
+                                        opExtent->contains(
+                                            NN_NO_CHECK(extentResolved)) &&
+                                        extentResolved->contains(
+                                            NN_NO_CHECK(opExtent));
                                     if (candidateBoundCRS) {
-                                        candidateBoundCRS = nullptr;
-                                        break;
+                                        if (exactlyMatchingExtent &&
+                                            !candidateHasExactlyMachingExtent) {
+                                            candidateBoundCRS = nullptr;
+                                        } else if (
+                                            exactlyMatchingExtent ==
+                                            candidateHasExactlyMachingExtent) {
+                                            candidateCount++;
+                                        }
+                                    }
+                                    if (candidateBoundCRS == nullptr) {
+                                        candidateCount = 1;
+                                        candidateHasExactlyMachingExtent =
+                                            exactlyMatchingExtent;
+                                        candidateBoundCRS =
+                                            BoundCRS::create(
+                                                thisAsCRS, hubCRS,
+                                                NN_NO_CHECK(transf))
+                                                .as_nullable();
                                     }
-                                    candidateBoundCRS =
-                                        BoundCRS::create(thisAsCRS, hubCRS,
-                                                         NN_NO_CHECK(transf))
-                                            .as_nullable();
                                 }
                             }
                         }
                     }
                 }
             }
-            if (candidateBoundCRS) {
+            if (candidateCount == 1 && candidateBoundCRS) {
                 return NN_NO_CHECK(candidateBoundCRS);
             }
         } catch (const std::exception &) {

diff --git a/test/cli/testprojinfo b/test/cli/testprojinfo
@@ -376,6 +376,13 @@ echo 'Testing --list-crs projected --area France | grep "RGF93 v1 / Lambert-93\|
 $EXE --list-crs projected --area France | grep "RGF93 v1 / Lambert-93\|UTM zone 11N" | sort >> ${OUT}
 echo "" >>${OUT}
 
+# Test North-Macedonia CRS for which there are several MGI 1901 -> WGS 84 Helmert transformations
+# but only one matches the extent of the CRS
+echo 'Testing EPSG:9945 -o PROJ -q' >> ${OUT}
+$EXE EPSG:9945 -o PROJ -q >> ${OUT}
+echo "" >>${OUT}
+
+
 # do 'diff' with distribution results
 echo "diff ${OUT} with testprojinfo_out.dist"
 diff -u ${OUT} ${TEST_CLI_DIR}/testprojinfo_out.dist

diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
@@ -1683,3 +1683,6 @@ EPSG:6512 "NAD83(2011) / Missouri East"
 Testing --list-crs projected --area France | grep "RGF93 v1 / Lambert-93\|UTM zone 11N" | sort
 EPSG:2154 "RGF93 v1 / Lambert-93"
 
+Testing EPSG:9945 -o PROJ -q
++proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=500000 +y_0=-4000000 +ellps=bessel +towgs84=521.748,229.489,590.921,4.029,4.488,-15.521,-9.78 +units=m +no_defs +type=crs
+