Merge pull request #2750 from rouault/better_operation_selection

proj_trans/cs2cs: If two operations have the same accuracy, use the one that is contained within a larger one

4 files changed, 33 insertions, 4 deletions

diff --git a/docs/source/operations/operations_computation.rst b/docs/source/operations/operations_computation.rst
index 0bc57ff0..5c2596e6 100644
--- a/docs/source/operations/operations_computation.rst
+++ b/docs/source/operations/operations_computation.rst
@@ -59,9 +59,9 @@ From a code point of view, the entry point of the algorithm is the C++
 
 It combines several strategies:
 
-- look up in the PROJ database for available operations
-- consider the pair (source CRS, target CRS) to synthetize operations depending
-  on the nature of the source and target CRS.
+  - look up in the PROJ database for available operations
+  - consider the pair (source CRS, target CRS) to synthetize operations depending
+    on the nature of the source and target CRS.
 
 Geographic CRS to Geographic CRS, with known identifiers
 --------------------------------------------------------
@@ -156,6 +156,18 @@ performed in the order they are listed below:
     lexicographic order (e.g. "FOO to BAR (3)" will have higher precedence than
     "FOO to BAR (2)")
 
+.. note::
+
+    :c:func:`proj_trans`, on the results returned by :c:func:`proj_create_crs_to_crs`,
+    will not necessarily use the operation that
+    is listed in first position due to the above algorithm. :c:func:`proj_trans`
+    has more context, since it has the coordinate to transform, so it can compare
+    this coordinate to the area of use of operations. Typically, the above criteria
+    will favor an operation that has a larger area of use over another one with a
+    smaller area, due to it being more generally applicable. But once coordinates are known,
+    :c:func:`proj_trans` can select an operation with a smaller
+    area of use that applies to the coordinate to transform.
+
 Geodetic/geographic CRS to Geodetic/geographic CRS, without known identifiers
 -----------------------------------------------------------------------------
 
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 8b676f9b..95582e48 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -243,7 +243,15 @@ int pj_get_suggested_operation(PJ_CONTEXT*,
             // onshore. So in a general way, prefer a onshore area to a
             // offshore one.
             if( iBest < 0 ||
-                (alt.accuracy >= 0 && alt.accuracy < bestAccuracy &&
+                (alt.accuracy >= 0 &&
+                (alt.accuracy < bestAccuracy ||
+                 // If two operations have the same accuracy, use the one that
+                 // is contained within a larger one
+                 (alt.accuracy == bestAccuracy &&
+                  alt.minxSrc > opList[iBest].minxSrc &&
+                  alt.minySrc > opList[iBest].minySrc &&
+                  alt.maxxSrc < opList[iBest].maxxSrc &&
+                  alt.maxySrc < opList[iBest].maxySrc)) &&
                 !alt.isOffshore) ) {
                 iBest = i;
                 bestAccuracy = alt.accuracy;
diff --git a/test/cli/testvarious b/test/cli/testvarious
index e88cf68b..6e9cd43c 100755
--- a/test/cli/testvarious
+++ b/test/cli/testvarious
@@ -1047,6 +1047,12 @@ $EXE -E EPSG:2636 "WGS 84" >> ${OUT} <<EOF
 5540944.47  500000.001
 EOF
 
+echo  "##############################################################" >> ${OUT}
+echo  "Check that we select the operation that has the smallest area of use, when 2 have the same accuracy" >> ${OUT}
+$EXE -E EPSG:4326 "NAD83(HARN)" >> ${OUT} <<EOF
+34 -120
+EOF
+
 
 # Done!
 # do 'diff' with distribution results
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index 44926677..9bb85933 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -504,3 +504,6 @@ Check that we can use a transformation spanning the antimeridian (should use Pul
 Check that we can use a transformation spanning the antimeridian (should use Pulkovo 1942 to WGS 84 (20))
 5540944.47  499999.999	49d59'59.36"N	179d59'52.133"W 0.000
 5540944.47  500000.001	49d59'59.36"N	179d59'52.133"W 0.000
+##############################################################
+Check that we select the operation that has the smallest area of use, when 2 have the same accuracy
+34 -120	33d59'59.983"N	119d59'59.955"W 0.000