Merge pull request #2466 from rouault/fix_2423

cs2cs: add --area and --bbox  options to restrict candidate coordinate operations (fixes #2423)

5 files changed, 214 insertions, 8 deletions

diff --git a/docs/source/apps/cs2cs.rst b/docs/source/apps/cs2cs.rst
index e214a5c0..7fe570d9 100644
--- a/docs/source/apps/cs2cs.rst
+++ b/docs/source/apps/cs2cs.rst
@@ -11,11 +11,15 @@ cs2cs
 Synopsis
 ********
 
-    **cs2cs** [**-eEfIlrstvwW** [args]] [*+opt[=arg]* ...] [+to *+opt[=arg]* ...] file ...
+    **cs2cs** [**-eEfIlrstvwW** [args]]
+              [[--area name_or_code] | [--bbox west_long,south_lat,east_long,north_lat]]
+              [*+opt[=arg]* ...] [+to *+opt[=arg]* ...] file ...
 
 or
 
-    **cs2cs** [**-eEfIlrstvwW** [args]] {source_crs} {target_crs} file ...
+    **cs2cs** [**-eEfIlrstvwW** [args]] [--area name_or_code]
+              [[--area name_or_code] | [--bbox west_long,south_lat,east_long,north_lat]]
+              {source_crs} {target_crs} file ...
 
     where {source_crs} or {target_crs} is one of the possibilities accepted
     by :c:func:`proj_create`, provided it expresses a CRS
@@ -144,6 +148,27 @@ The following control parameters can appear in any order:
     Causes a listing of cartographic control parameters tested for and used by
     the program to be printed prior to input data.
 
+.. option:: --area name_or_code
+
+    .. versionadded:: 8.0.0
+
+    Specify an area of interest to restrict the results when researching
+    coordinate operations between 2 CRS. The area of interest can be specified either
+    as a name (e.g "Denmark - onshore") or a AUTHORITY:CODE (EPSG:3237)
+
+    This option is mutually exclusive with :option:`--bbox`.
+
+.. option:: --bbox west_long,south_lat,east_long,north_lat
+
+    .. versionadded:: 8.0.0
+
+    Specify an area of interest to restrict the results when researching
+    coordinate operations between 2 CRS. The area of interest is specified as a
+    bounding box with geographic coordinates, expressed in degrees in a
+    unspecified geographic CRS.
+    `west_long` and `east_long` should be in the [-180,180] range, and
+    `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.
 
 .. only:: man
 
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 15bc73c8..9231a7ed 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1287,9 +1287,11 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
                                             area->east_lon_degree,
                                             area->north_lat_degree);
     }
+    else {
+        proj_operation_factory_context_set_spatial_criterion(
+            ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
+    }
 
-    proj_operation_factory_context_set_spatial_criterion(
-        ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
     proj_operation_factory_context_set_grid_availability_use(
         ctx, operation_ctx,
         proj_context_is_network_enabled(ctx) ?
diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
index 12e045bb..dd65baf4 100644
--- a/src/apps/cs2cs.cpp
+++ b/src/apps/cs2cs.cpp
@@ -36,8 +36,13 @@
 #include <string.h>
 
 #include <cassert>
+#include <iostream>
 #include <string>
 
+#include <proj/io.hpp>
+#include <proj/metadata.hpp>
+#include <proj/util.hpp>
+
 #include <proj/internal/internal.hpp>
 
 // PROJ include order is sensitive
@@ -70,12 +75,18 @@ static const char *oform =
 static char oform_buffer[16]; /* buffer for oform when using -d */
 static const char *oterr = "*\t*"; /* output line for unprojectable input */
 static const char *usage =
-    "%s\nusage: %s [-dDeEfIlrstvwW [args]] [+opt[=arg] ...]\n"
-    "                   [+to +opt[=arg] ...] [file ...]\n";
+    "%s\nusage: %s [-dDeEfIlrstvwW [args]]\n"
+    "              [[--area name_or_code] | [--bbox west_long,south_lat,east_long,north_lat]]\n"
+    "              [+opt[=arg] ...] [+to +opt[=arg] ...] [file ...]\n";
 
 static double (*informat)(const char *,
                           char **); /* input data deformatter function */
 
+using namespace NS_PROJ::io;
+using namespace NS_PROJ::metadata;
+using namespace NS_PROJ::util;
+using namespace NS_PROJ::internal;
+
 /************************************************************************/
 /*                              process()                               */
 /*                                                                      */
@@ -359,9 +370,47 @@ int main(int argc, char **argv) {
         }
     }
 
+    ExtentPtr bboxFilter;
+    std::string area;
+
     /* process run line arguments */
     while (--argc > 0) { /* collect run line arguments */
-        if (**++argv == '-') {
+        ++argv;
+        if (strcmp(*argv, "--area") == 0 ) {
+            ++argv;
+            --argc;
+            if( argc == 0 ) {
+                emess(1, "missing argument for --area");
+                std::exit(1);
+            }
+            area = *argv;
+        }
+        else if (strcmp(*argv, "--bbox") == 0) {
+            ++argv;
+            --argc;
+            if( argc == 0 ) {
+                emess(1, "missing argument for --bbox");
+                std::exit(1);
+            }
+            auto bboxStr(*argv);
+            auto bbox(split(bboxStr, ','));
+            if (bbox.size() != 4) {
+                std::cerr << "Incorrect number of values for option --bbox: "
+                          << bboxStr << std::endl;
+                std::exit(1);
+            }
+            try {
+                bboxFilter = Extent::createFromBBOX(
+                                 c_locale_stod(bbox[0]), c_locale_stod(bbox[1]),
+                                 c_locale_stod(bbox[2]), c_locale_stod(bbox[3]))
+                                 .as_nullable();
+            } catch (const std::exception &e) {
+                std::cerr << "Invalid value for option --bbox: " << bboxStr
+                          << ", " << e.what() << std::endl;
+                std::exit(1);
+            }
+        }
+        else if (**argv == '-') {
             for (arg = *argv;;) {
                 switch (*++arg) {
                 case '\0': /* position of "stdin" */
@@ -536,6 +585,102 @@ int main(int argc, char **argv) {
         }
     }
 
+    if (bboxFilter && !area.empty()) {
+        std::cerr << "ERROR: --bbox and --area are exclusive" << std::endl;
+        std::exit(1);
+    }
+
+    PJ_AREA* pj_area = nullptr;
+    if (!area.empty()) {
+
+        DatabaseContextPtr dbContext;
+        try {
+            dbContext =
+                DatabaseContext::create().as_nullable();
+        } catch (const std::exception &e) {
+            std::cerr << "ERROR: Cannot create database connection: "
+                    << e.what() << std::endl;
+            std::exit(1);
+        }
+
+        // Process area of use
+        try {
+            if (area.find(' ') == std::string::npos &&
+                area.find(':') != std::string::npos) {
+                auto tokens = split(area, ':');
+                if (tokens.size() == 2) {
+                    const std::string &areaAuth = tokens[0];
+                    const std::string &areaCode = tokens[1];
+                    bboxFilter = AuthorityFactory::create(
+                                        NN_NO_CHECK(dbContext), areaAuth)
+                                        ->createExtent(areaCode)
+                                        .as_nullable();
+                }
+            }
+            if (!bboxFilter) {
+                auto authFactory = AuthorityFactory::create(
+                    NN_NO_CHECK(dbContext), std::string());
+                auto res = authFactory->listAreaOfUseFromName(area, false);
+                if (res.size() == 1) {
+                    bboxFilter =
+                        AuthorityFactory::create(NN_NO_CHECK(dbContext),
+                                                    res.front().first)
+                            ->createExtent(res.front().second)
+                            .as_nullable();
+                } else {
+                    res = authFactory->listAreaOfUseFromName(area, true);
+                    if (res.size() == 1) {
+                        bboxFilter =
+                            AuthorityFactory::create(NN_NO_CHECK(dbContext),
+                                                        res.front().first)
+                                ->createExtent(res.front().second)
+                                .as_nullable();
+                    } else if (res.empty()) {
+                        std::cerr << "No area of use matching provided name"
+                                    << std::endl;
+                        std::exit(1);
+                    } else {
+                        std::cerr << "Several candidates area of use "
+                                        "matching provided name :"
+                                    << std::endl;
+                        for (const auto &candidate : res) {
+                            auto obj =
+                                AuthorityFactory::create(
+                                    NN_NO_CHECK(dbContext), candidate.first)
+                                    ->createExtent(candidate.second);
+                            std::cerr << "  " << candidate.first << ":"
+                                        << candidate.second << " : "
+                                        << *obj->description() << std::endl;
+                        }
+                        std::exit(1);
+                    }
+                }
+            }
+        } catch (const std::exception &e) {
+            std::cerr << "Area of use retrieval failed: " << e.what()
+                        << std::endl;
+            std::exit(1);
+        }
+    }
+
+    if (bboxFilter) {
+        auto geogElts = bboxFilter->geographicElements();
+        if (geogElts.size() == 1)
+        {
+            auto bbox = std::dynamic_pointer_cast<GeographicBoundingBox>(
+                geogElts[0].as_nullable());
+            if (bbox)
+            {
+                pj_area = proj_area_create();
+                proj_area_set_bbox(pj_area,
+                                   bbox->westBoundLongitude(),
+                                   bbox->southBoundLatitude(),
+                                   bbox->eastBoundLongitude(),
+                                   bbox->northBoundLatitude());
+            }
+        }
+    }
+
     /*
      * If the user has requested inverse, then just reverse the
      * coordinate systems.
@@ -617,10 +762,11 @@ int main(int argc, char **argv) {
     }
 
     transformation = proj_create_crs_to_crs_from_pj(nullptr, src, dst,
-                                                    nullptr, nullptr);
+                                                    pj_area, nullptr);
 
     proj_destroy(src);
     proj_destroy(dst);
+    proj_area_destroy(pj_area);
 
     if (!transformation) {
         emess(3, "cannot initialize transformation\ncause: %s",
diff --git a/test/cli/td_out.dist b/test/cli/td_out.dist
index ab0c0911..6bad8e57 100644
--- a/test/cli/td_out.dist
+++ b/test/cli/td_out.dist
@@ -7,6 +7,18 @@ As above, but without ntv1 everything goes through conus file.
 111d00'00.000"W 44d00'00.000"N 0.0	111d0'2.788"W	43d59'59.725"N 0.000
 111d00'00.000"W 39d00'00.000"N 0.0	111d0'2.604"W	38d59'59.912"N 0.000
 ##############################################################
+Test --area Canada NAD27 NAD83 (using ntv1_can)
+43d59'59.732"N	111d0'3.208"W 0.000
+*	* inf
+##############################################################
+Test --bbox -141.01,40.04,-47.74,86.46 NAD27 NAD83 (using ntv1_can)
+43d59'59.732"N	111d0'3.208"W 0.000
+*	* inf
+##############################################################
+Test --area "USA - CONUS - onshore" NAD27 NAD83 (using conus)
+43d59'59.725"N	111d0'2.788"W 0.000
+38d59'59.912"N	111d0'2.604"W 0.000
+##############################################################
 Test MD used where available
 79d58'00.000"W 37d02'00.000"N 0.0	79d58'0.005"W	37d1'59.998"N 0.000
 79d58'00.000"W 36d58'00.000"N 0.0	79d57'59.128"W	36d58'0.501"N 0.000
diff --git a/test/cli/testdatumfile b/test/cli/testdatumfile
index 16e4bbc3..5a013f12 100755
--- a/test/cli/testdatumfile
+++ b/test/cli/testdatumfile
@@ -60,6 +60,27 @@ $EXE +proj=latlong +ellps=clrk66 '+nadgrids="./dir with "" space/myconus"' \
 EOF
 
 echo "##############################################################" >> ${OUT}
+echo "Test --area Canada NAD27 NAD83 (using ntv1_can)" >> ${OUT}
+$EXE --area Canada NAD27 NAD83 >>${OUT} <<EOF
+44 -111
+39 -111
+EOF
+
+echo "##############################################################" >> ${OUT}
+echo "Test --bbox -141.01,40.04,-47.74,86.46 NAD27 NAD83 (using ntv1_can)" >> ${OUT}
+$EXE --bbox -141.01,40.04,-47.74,86.46 NAD27 NAD83 >>${OUT} <<EOF
+44 -111
+39 -111
+EOF
+
+echo "##############################################################" >> ${OUT}
+echo "Test --area \"USA - CONUS - onshore\" NAD27 NAD83 (using conus)" >> ${OUT}
+$EXE --area "USA - CONUS - onshore" NAD27 NAD83 >>${OUT} <<EOF
+44 -111
+39 -111
+EOF
+
+echo "##############################################################" >> ${OUT}
 echo Test MD used where available >> ${OUT}
 #
 $EXE +proj=latlong +ellps=clrk66 +nadgrids=MD,conus \