Merge pull request #2488 from rouault/crs_to_crs_improved_filtering

Improved coordinate operation filtering

11 files changed, 178 insertions, 15 deletions

diff --git a/docs/source/apps/cs2cs.rst b/docs/source/apps/cs2cs.rst
index 7df8890f..4706459a 100644
--- a/docs/source/apps/cs2cs.rst
+++ b/docs/source/apps/cs2cs.rst
@@ -13,7 +13,7 @@ Synopsis
 
     | **cs2cs** [**-eEfIlrstvwW** [args]]
     |           [[--area <name_or_code>] | [--bbox <west_long,south_lat,east_long,north_lat>]]
-    |           [--authority <name>]
+    |           [--authority <name>] [--no-ballpark] [--accuracy <accuracy>]
     |           ([*+opt[=arg]* ...] [+to *+opt[=arg]* ...] | {source_crs} {target_crs})
     |           file ...
 
@@ -166,6 +166,19 @@ The following control parameters can appear in any order:
     `south_lat` and `north_lat` in the [-90,90]. `west_long` is generally lower than
     `east_long`, except in the case where the area of interest crosses the antimeridian.
 
+.. option:: --no-ballpark
+
+    .. versionadded:: 8.0.0
+
+    Disallow any coordinate operation that is, or contains, a
+    :term:`Ballpark transformation`
+
+.. option:: --accuracy <accuracy>
+
+    .. versionadded:: 8.0.0
+
+    Sets the minimum desired accuracy for candidate coordinate operations.
+
 .. option:: --authority <name>
 
     .. versionadded:: 8.0.0
diff --git a/docs/source/apps/projinfo.rst b/docs/source/apps/projinfo.rst
index 803c0a65..f0f66620 100644
--- a/docs/source/apps/projinfo.rst
+++ b/docs/source/apps/projinfo.rst
@@ -22,7 +22,7 @@ Synopsis
     |    [--crs-extent-use none|both|intersection|smallest]
     |    [--grid-check none|discard_missing|sort|known_available]
     |    [--pivot-crs always|if_no_direct_transformation|never|{auth:code[,auth:code]*}]
-    |    [--show-superseded] [--hide-ballpark]
+    |    [--show-superseded] [--hide-ballpark] [--accuracy {accuracy}]
     |    [--allow-ellipsoidal-height-as-vertical-crs]
     |    [--boundcrs-to-wgs84]
     |    [--main-db-path path] [--aux-db-path path]*
@@ -213,6 +213,14 @@ The following control parameters can appear in any order:
 
     .. note:: only used for coordinate operation computation
 
+.. option:: --accuracy {accuracy}
+
+    .. versionadded:: 8.0
+
+    Sets the minimum desired accuracy for returned coordinate operations.
+
+    .. note:: only used for coordinate operation computation
+
 .. option:: --allow-ellipsoidal-height-as-vertical-crs
 
     .. versionadded:: 8.0
diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst
index 08b3e669..34fe59ed 100644
--- a/docs/source/development/reference/functions.rst
+++ b/docs/source/development/reference/functions.rst
@@ -195,6 +195,11 @@ paragraph for more details.
       If authority is a non-empty string different of ``any``, then coordinate operations
       will be searched only in that authority namespace (e.g ``EPSG``).
 
+    - ACCURACY=value: to set the minimum desired accuracy (in metres) of the
+      candidate coordinate operations.
+
+    - ALLOW_BALLPARK=YES/NO: can be set to NO to disallow the use of
+      :term:`Ballpark transformation` in the candidate coordinate operations.
 
 .. doxygenfunction:: proj_normalize_for_visualization
    :project: doxygen_api
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 7423240f..909c3c32 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1302,10 +1302,24 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
     }
 
     const char* authority = nullptr;
+    double accuracy = -1;
+    bool allowBallparkTransformations = true;
     for (auto iter = options; iter && iter[0]; ++iter) {
         const char *value;
         if ((value = getOptionValue(*iter, "AUTHORITY="))) {
             authority = value;
+        } else if ((value = getOptionValue(*iter, "ACCURACY="))) {
+            accuracy = pj_atof(value);
+        } else if ((value = getOptionValue(*iter, "ALLOW_BALLPARK="))) {
+            if( ci_equal(value, "yes") )
+                allowBallparkTransformations = true;
+            else if( ci_equal(value, "no") )
+                allowBallparkTransformations = false;
+            else {
+                ctx->logger(ctx->logger_app_data, PJ_LOG_ERROR,
+                            "Invalid value for ALLOW_BALLPARK option.");
+                return nullptr;
+            }
         } else {
             std::string msg("Unknown option :");
             msg += *iter;
@@ -1319,6 +1333,14 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
         return nullptr;
     }
 
+    proj_operation_factory_context_set_allow_ballpark_transformations(
+        ctx, operation_ctx, allowBallparkTransformations);
+
+    if( accuracy >= 0 ) {
+        proj_operation_factory_context_set_desired_accuracy(ctx, operation_ctx,
+                                                            accuracy);
+    }
+
     if( area && area->bbox_set ) {
         proj_operation_factory_context_set_area_of_interest(
                                             ctx,
diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
index 409a5ef3..5542a282 100644
--- a/src/apps/cs2cs.cpp
+++ b/src/apps/cs2cs.cpp
@@ -38,6 +38,7 @@
 #include <cassert>
 #include <iostream>
 #include <string>
+#include <vector>
 
 #include <proj/io.hpp>
 #include <proj/metadata.hpp>
@@ -77,7 +78,7 @@ static const char *oterr = "*\t*"; /* output line for unprojectable input */
 static const char *usage =
     "%s\nusage: %s [-dDeEfIlrstvwW [args]]\n"
     "              [[--area name_or_code] | [--bbox west_long,south_lat,east_long,north_lat]]\n"
-    "              [--authority {name}]\n"
+    "              [--authority {name}] [--accuracy {accuracy}] [--no-ballpark]\n"
     "              [+opt[=arg] ...] [+to +opt[=arg] ...] [file ...]\n";
 
 static double (*informat)(const char *,
@@ -374,6 +375,8 @@ int main(int argc, char **argv) {
     ExtentPtr bboxFilter;
     std::string area;
     const char* authority = nullptr;
+    double accuracy = -1;
+    bool allowBallpark = true;
 
     /* process run line arguments */
     while (--argc > 0) { /* collect run line arguments */
@@ -412,6 +415,15 @@ int main(int argc, char **argv) {
                 std::exit(1);
             }
         }
+        else if (strcmp(*argv, "--accuracy") == 0 ) {
+            ++argv;
+            --argc;
+            if( argc == 0 ) {
+                emess(1, "missing argument for --accuracy");
+                std::exit(1);
+            }
+            accuracy = c_locale_stod(*argv);
+        }
         else if (strcmp(*argv, "--authority") == 0 ) {
             ++argv;
             --argc;
@@ -421,6 +433,9 @@ int main(int argc, char **argv) {
             }
             authority = *argv;
         }
+        else if (strcmp(*argv, "--no-ballpark") == 0 ) {
+            allowBallpark = false;
+        }
         else if (**argv == '-') {
             for (arg = *argv;;) {
                 switch (*++arg) {
@@ -773,14 +788,24 @@ int main(int argc, char **argv) {
     }
 
     std::string authorityOption; /* keep this variable in this outer scope ! */
-    const char* options[2] = { nullptr, nullptr };
+    std::string accuracyOption; /* keep this variable in this outer scope ! */
+    std::vector<const char*> options;
     if( authority ) {
         authorityOption = "AUTHORITY=";
         authorityOption += authority;
-        options[0] = authorityOption.data();
+        options.push_back(authorityOption.data());
+    }
+    if( accuracy >= 0 ) {
+        accuracyOption = "ACCURACY=";
+        accuracyOption += toString(accuracy);
+        options.push_back(accuracyOption.data());
+    }
+    if( !allowBallpark ) {
+        options.push_back("ALLOW_BALLPARK=NO");
     }
+    options.push_back(nullptr);
     transformation = proj_create_crs_to_crs_from_pj(nullptr, src, dst,
-                                                    pj_area, options);
+                                                    pj_area, options.data());
 
     proj_destroy(src);
     proj_destroy(dst);
diff --git a/src/apps/projinfo.cpp b/src/apps/projinfo.cpp
index da885fbb..8d389019 100644
--- a/src/apps/projinfo.cpp
+++ b/src/apps/projinfo.cpp
@@ -95,7 +95,9 @@ static void usage() {
         << std::endl
         << "                [--pivot-crs always|if_no_direct_transformation|"
         << "never|{auth:code[,auth:code]*}]" << std::endl
-        << "                [--show-superseded] [--hide-ballpark]" << std::endl
+        << "                [--show-superseded] [--hide-ballpark] "
+           "[--accuracy {accuracy}]"
+        << std::endl
         << "                [--allow-ellipsoidal-height-as-vertical-crs]"
         << std::endl
         << "                [--boundcrs-to-wgs84]" << std::endl
@@ -672,8 +674,8 @@ static void outputOperations(
     CoordinateOperationContext::IntermediateCRSUse allowUseIntermediateCRS,
     const std::vector<std::pair<std::string, std::string>> &pivots,
     const std::string &authority, bool usePROJGridAlternatives,
-    bool showSuperseded, bool promoteTo3D, const OutputOptions &outputOpt,
-    bool summary) {
+    bool showSuperseded, bool promoteTo3D, double minimumAccuracy,
+    const OutputOptions &outputOpt, bool summary) {
     auto sourceObj =
         buildObject(dbContext, sourceCRSStr, "crs", "source CRS", false,
                     CoordinateOperationContext::IntermediateCRSUse::NEVER,
@@ -715,6 +717,9 @@ static void outputOperations(
         ctxt->setUsePROJAlternativeGridNames(usePROJGridAlternatives);
         ctxt->setDiscardSuperseded(!showSuperseded);
         ctxt->setAllowBallparkTransformations(outputOpt.ballparkAllowed);
+        if (minimumAccuracy >= 0) {
+            ctxt->setDesiredAccuracy(minimumAccuracy);
+        }
         list = CoordinateOperationFactory::create()->createOperations(
             nnSourceCRS, nnTargetCRS, ctxt);
         if (!spatialCriterionExplicitlySpecified &&
@@ -819,6 +824,7 @@ int main(int argc, char **argv) {
     bool identify = false;
     bool showSuperseded = false;
     bool promoteTo3D = false;
+    double minimumAccuracy = -1;
 
     for (int i = 1; i < argc; i++) {
         std::string arg(argv[i]);
@@ -934,6 +940,9 @@ int main(int argc, char **argv) {
                           << ", " << e.what() << std::endl;
                 usage();
             }
+        } else if (arg == "--accuracy" && i + 1 < argc) {
+            i++;
+            minimumAccuracy = c_locale_stod(argv[i]);
         } else if (arg == "--area" && i + 1 < argc) {
             i++;
             area = argv[i];
@@ -1338,12 +1347,12 @@ int main(int argc, char **argv) {
         }
 
         try {
-            outputOperations(dbContext, sourceCRSStr, targetCRSStr, bboxFilter,
-                             spatialCriterion,
-                             spatialCriterionExplicitlySpecified, crsExtentUse,
-                             gridAvailabilityUse, allowUseIntermediateCRS,
-                             pivots, authority, usePROJGridAlternatives,
-                             showSuperseded, promoteTo3D, outputOpt, summary);
+            outputOperations(
+                dbContext, sourceCRSStr, targetCRSStr, bboxFilter,
+                spatialCriterion, spatialCriterionExplicitlySpecified,
+                crsExtentUse, gridAvailabilityUse, allowUseIntermediateCRS,
+                pivots, authority, usePROJGridAlternatives, showSuperseded,
+                promoteTo3D, minimumAccuracy, outputOpt, summary);
         } catch (const std::exception &e) {
             std::cerr << "outputOperations() failed with: " << e.what()
                       << std::endl;
diff --git a/test/cli/testprojinfo b/test/cli/testprojinfo
index c31cfef0..ee1b27f5 100755
--- a/test/cli/testprojinfo
+++ b/test/cli/testprojinfo
@@ -213,6 +213,10 @@ fi
 rm testprojinfo_out_remotedata.txt
 unset PROJ_NETWORK
 
+echo 'Testing --accuracy 0.05 -s EPSG:4326 -t EPSG:4258' >> ${OUT}
+$EXE --accuracy 0.05 -s EPSG:4326 -t EPSG:4258 >>${OUT} 2>&1
+echo "" >>${OUT}
+
 ######################
 # NZGD2000 -> ITRFxx #
 ######################
diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
index 8e8ef294..829c914c 100644
--- a/test/cli/testprojinfo_out.dist
+++ b/test/cli/testprojinfo_out.dist
@@ -1384,6 +1384,9 @@ DATUM["World Geodetic System 1984",
         LENGTHUNIT["metre",1]],
     ID["EPSG",6326]]
 
+Testing --accuracy 0.05 -s EPSG:4326 -t EPSG:4258
+Candidate operations found: 0
+
 Testing -s NZGD2000 -t ITRF96 -o PROJ -q
 +proj=pipeline
   +step +proj=unitconvert +xy_in=deg +xy_out=rad
diff --git a/test/cli/testvarious b/test/cli/testvarious
index 82be4992..a121393c 100755
--- a/test/cli/testvarious
+++ b/test/cli/testvarious
@@ -1021,6 +1021,18 @@ $EXE --authority EPSG -E +proj=latlong +datum=WGS84 +no_defs +to +init=epsg:6342
 -105 40
 EOF
 
+echo  "##############################################################" >> ${OUT}
+echo  "Test effect of --accuracy" >> ${OUT}
+$EXE -E --accuracy 0.05 EPSG:4326 EPSG:4258 >> ${OUT} <<EOF
+49 2
+EOF
+
+echo  "##############################################################" >> ${OUT}
+echo  "Test effect of --no-ballpark" >> ${OUT}
+$EXE -E --no-ballpark EPSG:4267 EPSG:4258 >> ${OUT} <<EOF
+49 2
+EOF
+
 
 # Done!
 # do 'diff' with distribution results
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index fe1aa452..59129d99 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -492,3 +492,7 @@ The first result should use the 'WGS_1984_(ITRF08)_To_NAD_1983_2011' (ESRI:10836
 and the second one a no-op
 -105 40	500000.86	4427756.50 0.00
 -105 40	500000.00	4427757.22 0.00
+##############################################################
+Test effect of --accuracy
+##############################################################
+Test effect of --no-ballpark
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index c417371d..e54c6fa7 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -4196,6 +4196,64 @@ TEST_F(CApi, proj_create_crs_to_crs_from_pj) {
 
 // ---------------------------------------------------------------------------
 
+TEST_F(CApi, proj_create_crs_to_crs_from_pj_accuracy_filter) {
+
+    auto src = proj_create(m_ctxt, "EPSG:4326"); // WGS 84
+    ObjectKeeper keeper_src(src);
+    ASSERT_NE(src, nullptr);
+
+    auto dst = proj_create(m_ctxt, "EPSG:4258"); // ETRS89
+    ObjectKeeper keeper_dst(dst);
+    ASSERT_NE(dst, nullptr);
+
+    // No options
+    {
+        auto P =
+            proj_create_crs_to_crs_from_pj(m_ctxt, src, dst, nullptr, nullptr);
+        ObjectKeeper keeper_P(P);
+        ASSERT_NE(P, nullptr);
+    }
+
+    {
+        const char *const options[] = {"ACCURACY=0.05", nullptr};
+        auto P =
+            proj_create_crs_to_crs_from_pj(m_ctxt, src, dst, nullptr, options);
+        ObjectKeeper keeper_P(P);
+        ASSERT_EQ(P, nullptr);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_create_crs_to_crs_from_pj_ballpark_filter) {
+
+    auto src = proj_create(m_ctxt, "EPSG:4267"); // NAD 27
+    ObjectKeeper keeper_src(src);
+    ASSERT_NE(src, nullptr);
+
+    auto dst = proj_create(m_ctxt, "EPSG:4258"); // ETRS89
+    ObjectKeeper keeper_dst(dst);
+    ASSERT_NE(dst, nullptr);
+
+    // No options
+    {
+        auto P =
+            proj_create_crs_to_crs_from_pj(m_ctxt, src, dst, nullptr, nullptr);
+        ObjectKeeper keeper_P(P);
+        ASSERT_NE(P, nullptr);
+    }
+
+    {
+        const char *const options[] = {"ALLOW_BALLPARK=NO", nullptr};
+        auto P =
+            proj_create_crs_to_crs_from_pj(m_ctxt, src, dst, nullptr, options);
+        ObjectKeeper keeper_P(P);
+        ASSERT_EQ(P, nullptr);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 static void
 check_axis_is_latitude(PJ_CONTEXT *ctx, PJ *cs, int axis_number,
                        const char *unit_name = "degree",