From 0b4c8b0a76c65d8b1146224bf314cad53e9f6607 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 28 Nov 2020 15:51:27 +0100 Subject: cs2cs: add --area and --bbox options to restrict candidate coordinate operations (fixes #2423) --- docs/source/apps/cs2cs.rst | 29 ++++++++- src/apps/cs2cs.cpp | 154 +++++++++++++++++++++++++++++++++++++++++++-- test/cli/td_out.dist | 12 ++++ test/cli/testdatumfile | 21 +++++++ 4 files changed, 210 insertions(+), 6 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/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 #include +#include #include +#include +#include +#include + #include // 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( + 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 @@ -59,6 +59,27 @@ $EXE +proj=latlong +ellps=clrk66 '+nadgrids="./dir with "" space/myconus"' \ 111d00'00.000"W 39d00'00.000"N 0.0 EOF +echo "##############################################################" >> ${OUT} +echo "Test --area Canada NAD27 NAD83 (using ntv1_can)" >> ${OUT} +$EXE --area Canada NAD27 NAD83 >>${OUT} <> ${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} <> ${OUT} +echo "Test --area \"USA - CONUS - onshore\" NAD27 NAD83 (using conus)" >> ${OUT} +$EXE --area "USA - CONUS - onshore" NAD27 NAD83 >>${OUT} <> ${OUT} echo Test MD used where available >> ${OUT} # -- cgit v1.2.3 From ebe3425bf66287e004958eb53976d3837f88b9e1 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 28 Nov 2020 15:57:17 +0100 Subject: proj_create_crs_to_crs_from_pj(): do not use PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION if area is specified --- src/4D_api.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) 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) ? -- cgit v1.2.3