From 1ac2c960058a5dfb47c1fc7b36b3357b6f25fbcb Mon Sep 17 00:00:00 2001 From: Javier Jimenez Shaw Date: Sun, 16 Jul 2023 12:48:16 +0200 Subject: [PATCH 1/7] Add option in proj CLI to use a CRS --- docs/source/apps/proj.rst | 12 +++++++++-- src/apps/proj.cpp | 45 ++++++++++++++++++++++++++++++++++++++- 2 files changed, 54 insertions(+), 3 deletions(-) diff --git a/docs/source/apps/proj.rst b/docs/source/apps/proj.rst index 6f4a1866d5..2552388173 100644 --- a/docs/source/apps/proj.rst +++ b/docs/source/apps/proj.rst @@ -18,9 +18,9 @@ invproj Synopsis ******** - **proj** [**-beEfiIlmorsStTvVwW**] [args]] [*+opt[=arg]* ...] file ... + **proj** [**-beEfiIlmorsStTvVwWC**] [args]] [*+opt[=arg]* ...] file ... - **invproj** [**-beEfiIlmorsStTvVwW**] [args]] [*+opt[=arg]* ...] file ... + **invproj** [**-beEfiIlmorsStTvVwWC**] [args]] [*+opt[=arg]* ...] file ... Description @@ -157,6 +157,14 @@ The following control parameters can appear in any order This option causes an expanded annotated listing of the characteristics of the projected point. :option:`-v` is implied with this option. +.. option:: -C + +.. versionadded:: 9.3.0 + + Where *crs* is a definition of a Projected CRS, as its WKT or code + for instance EPSG:32632. This option generates the projection to + the underlying geographic system. When used, it ignores *+opt* arguments. + The *+opt* run-line arguments are associated with cartographic parameters. Additional projection control parameters may be contained in two auxiliary diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp index d5afaa4909..161fea3a0e 100644 --- a/src/apps/proj.cpp +++ b/src/apps/proj.cpp @@ -1,6 +1,7 @@ /* <<<< Cartographic projection filter program >>>> */ #include "proj.h" #include "emess.h" +#include "proj_experimental.h" #include "proj_internal.h" #include "utils.h" #include @@ -48,6 +49,8 @@ static const char *oterr = "*\t*", /* output line for unprojectable input */ *usage = "%s\nusage: %s [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] " "[file ...]\n"; +static const char *ocrs = nullptr; /* CRS use case */ + static PJ_FACTORS facs; static double (*informat)(const char *, @@ -496,6 +499,11 @@ int main(int argc, char **argv) { case 's': /* reverse output */ reverseout = 1; continue; + case 'C': /* CRS use case */ + if (--argc <= 0) + goto noargument; + ocrs = *++argv; + continue; default: emess(1, "invalid option: -%c", *arg); break; @@ -525,11 +533,46 @@ int main(int argc, char **argv) { } proj_context_use_proj4_init_rules(nullptr, true); + if (ocrs) { + // logic copied from proj_factors function + if (PJ *P = proj_create(nullptr, ocrs)) { + const auto type = proj_get_type(P); + if (type == PJ_TYPE_PROJECTED_CRS) { + auto ctx = P->ctx; + auto geodetic_crs = proj_get_source_crs(ctx, P); + assert(geodetic_crs); + auto datum = proj_crs_get_datum(ctx, geodetic_crs); + auto datum_ensemble = + proj_crs_get_datum_ensemble(ctx, geodetic_crs); + auto cs = proj_create_ellipsoidal_2D_cs( + ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0); + auto temp = proj_create_geographic_crs_from_datum( + ctx, "unnamed crs", datum ? datum : datum_ensemble, cs); + proj_destroy(datum); + proj_destroy(datum_ensemble); + proj_destroy(cs); + proj_destroy(geodetic_crs); + Proj = proj_create_crs_to_crs_from_pj(ctx, temp, P, nullptr, + nullptr); + proj_destroy(temp); + } else { + emess(3, "CRS must be projected"); + } + proj_destroy(P); + } else { + emess(3, "-C argument is not parseable"); + } + if (!argvVector.empty()) { + emess(-1, "+opt arguments are ignored due to -C option"); + } + } + // proj historically ignores any datum shift specifier, like nadgrids, // towgs84, etc argvVector.push_back(const_cast("break_cs2cs_recursion")); - if (!(Proj = proj_create_argv(nullptr, static_cast(argvVector.size()), + if (!Proj && + !(Proj = proj_create_argv(nullptr, static_cast(argvVector.size()), argvVector.data()))) emess(3, "projection initialization failure\ncause: %s", proj_errno_string(proj_context_errno(nullptr))); From 3dfbed7e671203c0d6a0edb0eac086f9fcb4205b Mon Sep 17 00:00:00 2001 From: Javier Jimenez Shaw Date: Sun, 16 Jul 2023 13:19:12 +0200 Subject: [PATCH 2/7] add test for proj cli with -C option --- test/cli/testproj | 6 ++++++ test/cli/testproj_out.dist | 1 + 2 files changed, 7 insertions(+) diff --git a/test/cli/testproj b/test/cli/testproj index 9086f2d3a3..070d1e79a5 100755 --- a/test/cli/testproj +++ b/test/cli/testproj @@ -35,6 +35,12 @@ echo "doing tests into file ${OUT}, please wait" $EXE +ellps=WGS84 +proj=ob_tran +o_proj=latlon +o_lon_p=0.0 +o_lat_p=90.0 +lon_0=360.0 +to_meter=0.0174532925199433 +no_defs -E -f '%.3f' >${OUT} <>${OUT} < From 05ccb727463800b5d4956e31ec8b564fe38c536c Mon Sep 17 00:00:00 2001 From: Javier Jimenez Shaw Date: Sun, 16 Jul 2023 20:50:18 +0200 Subject: [PATCH 3/7] choose more numerically stable test case --- test/cli/testproj | 4 ++-- test/cli/testproj_out.dist | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/cli/testproj b/test/cli/testproj index 070d1e79a5..a9fe62995c 100755 --- a/test/cli/testproj +++ b/test/cli/testproj @@ -38,8 +38,8 @@ EOF # echo "Test -C option" # -$EXE -C EPSG:26948 -S >>${OUT} <>${OUT} < +500000.00 0.00 <0.9996 0.9996 0.9992 0 0.9996 0.9996> From 9adf4b5d96d83392fba1e35f791638b929f31b3d Mon Sep 17 00:00:00 2001 From: Javier Jimenez Shaw Date: Mon, 17 Jul 2023 23:48:13 +0200 Subject: [PATCH 4/7] Read CRS directly from the arguments, without option flag. --- docs/source/apps/proj.rst | 19 +++++++++---------- src/apps/proj.cpp | 24 +++++++++--------------- test/cli/testproj | 4 ++-- 3 files changed, 20 insertions(+), 27 deletions(-) diff --git a/docs/source/apps/proj.rst b/docs/source/apps/proj.rst index 2552388173..2183bd0d20 100644 --- a/docs/source/apps/proj.rst +++ b/docs/source/apps/proj.rst @@ -18,9 +18,9 @@ invproj Synopsis ******** - **proj** [**-beEfiIlmorsStTvVwWC**] [args]] [*+opt[=arg]* ...] file ... + **proj** [**-beEfiIlmorsStTvVwW**] [args]] ([*+opt[=arg]* ...] | {crs}) file ... - **invproj** [**-beEfiIlmorsStTvVwWC**] [args]] [*+opt[=arg]* ...] file ... + **invproj** [**-beEfiIlmorsStTvVwW**] [args]] ([*+opt[=arg]* ...] | {crs}) file ... Description @@ -157,14 +157,6 @@ The following control parameters can appear in any order This option causes an expanded annotated listing of the characteristics of the projected point. :option:`-v` is implied with this option. -.. option:: -C - -.. versionadded:: 9.3.0 - - Where *crs* is a definition of a Projected CRS, as its WKT or code - for instance EPSG:32632. This option generates the projection to - the underlying geographic system. When used, it ignores *+opt* arguments. - The *+opt* run-line arguments are associated with cartographic parameters. Additional projection control parameters may be contained in two auxiliary @@ -180,6 +172,13 @@ also used for supporting files like datum shift files. Usage of *+opt* varies with projection and for a complete description consult the :ref:`projection pages `. +.. versionadded:: 9.3.0 + + *{crs}* is one of the possibilities accepted by proj_create(), provided it + expresses a projected CRS, like a WKT string, an object code (like "EPSG:32632") + a PROJJSON string, etc. + The projection computed will be those of the map projection implied by + the transformation from the base geographic CRS of the projected CRS to the projected CRS. One or more files (processed in left to right order) specify the source of data to be converted. A ``-`` will specify the location of processing standard diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp index 161fea3a0e..a6ced784a9 100644 --- a/src/apps/proj.cpp +++ b/src/apps/proj.cpp @@ -49,8 +49,6 @@ static const char *oterr = "*\t*", /* output line for unprojectable input */ *usage = "%s\nusage: %s [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] " "[file ...]\n"; -static const char *ocrs = nullptr; /* CRS use case */ - static PJ_FACTORS facs; static double (*informat)(const char *, @@ -499,11 +497,6 @@ int main(int argc, char **argv) { case 's': /* reverse output */ reverseout = 1; continue; - case 'C': /* CRS use case */ - if (--argc <= 0) - goto noargument; - ocrs = *++argv; - continue; default: emess(1, "invalid option: -%c", *arg); break; @@ -515,8 +508,6 @@ int main(int argc, char **argv) { } else /* assumed to be input file name(s) */ eargv[eargc++] = *argv; } - if (eargc == 0) /* if no specific files force sysin */ - eargv[eargc++] = const_cast("-"); if (oform) { if (!validate_form_string_for_numbers(oform)) { @@ -533,9 +524,13 @@ int main(int argc, char **argv) { } proj_context_use_proj4_init_rules(nullptr, true); - if (ocrs) { + if (argvVector.empty() && eargc >= 1) { + // Consider the next arg as a CRS, not a file. + std::string ocrs = eargv[0]; + eargv++; + eargc--; // logic copied from proj_factors function - if (PJ *P = proj_create(nullptr, ocrs)) { + if (PJ *P = proj_create(nullptr, ocrs.c_str())) { const auto type = proj_get_type(P); if (type == PJ_TYPE_PROJECTED_CRS) { auto ctx = P->ctx; @@ -560,12 +555,11 @@ int main(int argc, char **argv) { } proj_destroy(P); } else { - emess(3, "-C argument is not parseable"); - } - if (!argvVector.empty()) { - emess(-1, "+opt arguments are ignored due to -C option"); + emess(3, "CRS is not parseable"); } } + if (eargc == 0) /* if no specific files force sysin */ + eargv[eargc++] = const_cast("-"); // proj historically ignores any datum shift specifier, like nadgrids, // towgs84, etc diff --git a/test/cli/testproj b/test/cli/testproj index a9fe62995c..375f53cdbf 100755 --- a/test/cli/testproj +++ b/test/cli/testproj @@ -36,9 +36,9 @@ $EXE +ellps=WGS84 +proj=ob_tran +o_proj=latlon +o_lon_p=0.0 +o_lat_p=90.0 +lon_0 2 49 EOF # -echo "Test -C option" +echo "Test CRS option" # -$EXE -C EPSG:32620 -S >>${OUT} <>${OUT} < Date: Thu, 20 Jul 2023 09:21:19 +0200 Subject: [PATCH 5/7] pass metre and e-n axis order to proj_factors --- src/apps/proj.cpp | 53 +++++++++++++++++++++++++++++++++-------------- 1 file changed, 37 insertions(+), 16 deletions(-) diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp index a6ced784a9..7b2e2a80dd 100644 --- a/src/apps/proj.cpp +++ b/src/apps/proj.cpp @@ -23,7 +23,8 @@ #define MAX_LINE 1000 #define PJ_INVERSE(P) (P->inv ? 1 : 0) -static PJ *Proj; +static PJ *Proj = nullptr; +static PJ *ProjForFactors = nullptr; static union { PJ_XY (*fwd)(PJ_LP, PJ *); PJ_LP (*inv)(PJ_XY, PJ *); @@ -116,16 +117,16 @@ static void process(FILE *fid) { data.uv.v *= fscale; } if (dofactors && !inverse) { - facs = proj_factors(Proj, coord); - facs_bad = proj_errno(Proj); + facs = proj_factors(ProjForFactors, coord); + facs_bad = proj_errno(ProjForFactors); } const auto xy = (*proj.fwd)(data.lp, Proj); data.xy = xy; if (dofactors && inverse) { - facs = proj_factors(Proj, coord); - facs_bad = proj_errno(Proj); + facs = proj_factors(ProjForFactors, coord); + facs_bad = proj_errno(ProjForFactors); } if (postscale && data.uv.u != HUGE_VAL) { @@ -282,8 +283,8 @@ static void vprocess(FILE *fid) { if (!*s && (s > line)) --s; /* assumed we gobbled \n */ coord.lp = dat_ll; - facs = proj_factors(Proj, coord); - if (proj_errno(Proj)) { + facs = proj_factors(ProjForFactors, coord); + if (proj_errno(ProjForFactors)) { emess(-1, "failed to compute factors\n\n"); continue; } @@ -541,15 +542,29 @@ int main(int argc, char **argv) { proj_crs_get_datum_ensemble(ctx, geodetic_crs); auto cs = proj_create_ellipsoidal_2D_cs( ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0); - auto temp = proj_create_geographic_crs_from_datum( + auto geogCRSNormalized = proj_create_geographic_crs_from_datum( ctx, "unnamed crs", datum ? datum : datum_ensemble, cs); proj_destroy(datum); proj_destroy(datum_ensemble); proj_destroy(cs); + Proj = proj_create_crs_to_crs_from_pj(ctx, geogCRSNormalized, P, + nullptr, nullptr); + + auto conversion = proj_crs_get_coordoperation(ctx, P); + auto projCS = proj_create_cartesian_2D_cs( + ctx, PJ_CART2D_EASTING_NORTHING, "metre", 1.0); + auto projCRSNormalized = proj_create_projected_crs( + ctx, nullptr, geodetic_crs, conversion, projCS); + assert(projCRSNormalized); proj_destroy(geodetic_crs); - Proj = proj_create_crs_to_crs_from_pj(ctx, temp, P, nullptr, - nullptr); - proj_destroy(temp); + proj_destroy(conversion); + proj_destroy(projCS); + ProjForFactors = proj_create_crs_to_crs_from_pj( + ctx, geogCRSNormalized, projCRSNormalized, nullptr, + nullptr); + + proj_destroy(geogCRSNormalized); + proj_destroy(projCRSNormalized); } else { emess(3, "CRS must be projected"); } @@ -565,11 +580,15 @@ int main(int argc, char **argv) { // towgs84, etc argvVector.push_back(const_cast("break_cs2cs_recursion")); - if (!Proj && - !(Proj = proj_create_argv(nullptr, static_cast(argvVector.size()), - argvVector.data()))) - emess(3, "projection initialization failure\ncause: %s", - proj_errno_string(proj_context_errno(nullptr))); + if (!Proj) { + if (!(Proj = + proj_create_argv(nullptr, static_cast(argvVector.size()), + argvVector.data()))) + emess(3, "projection initialization failure\ncause: %s", + proj_errno_string(proj_context_errno(nullptr))); + + ProjForFactors = Proj; + } if (!proj_angular_input(Proj, PJ_FWD)) { emess(3, "can't initialize operations that take non-angular input " @@ -653,6 +672,8 @@ int main(int argc, char **argv) { emess_dat.File_name = nullptr; } + if (ProjForFactors && ProjForFactors != Proj) + proj_destroy(ProjForFactors); if (Proj) proj_destroy(Proj); From 807ae5ac315f90d99e734422f9d8a0a65cd634cf Mon Sep 17 00:00:00 2001 From: Javier Jimenez Shaw Date: Fri, 21 Jul 2023 01:38:50 +0200 Subject: [PATCH 6/7] consider northing-easing in verbose mode --- src/apps/proj.cpp | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp index 7b2e2a80dd..c2783bd1f8 100644 --- a/src/apps/proj.cpp +++ b/src/apps/proj.cpp @@ -10,6 +10,8 @@ #include #include +#include + #include #if defined(MSDOS) || defined(OS2) || defined(WIN32) || defined(__WIN32__) @@ -25,6 +27,7 @@ static PJ *Proj = nullptr; static PJ *ProjForFactors = nullptr; +static bool swapAxisCrs = false; static union { PJ_XY (*fwd)(PJ_LP, PJ *); PJ_LP (*inv)(PJ_XY, PJ *); @@ -300,10 +303,12 @@ static void vprocess(FILE *fid) { (void)fputs(proj_rtodms2(pline, sizeof(pline), dat_ll.phi, 'N', 'S'), stdout); (void)printf(" [ %.11g ]\n", dat_ll.phi * RAD_TO_DEG); - (void)fputs("Easting (x): ", stdout); + (void)fputs(swapAxisCrs ? "Northing (y): " : "Easting (x): ", + stdout); (void)printf(oform, dat_xy.x); putchar('\n'); - (void)fputs("Northing (y): ", stdout); + (void)fputs(swapAxisCrs ? "Easting (x): " : "Northing (y): ", + stdout); (void)printf(oform, dat_xy.y); putchar('\n'); (void)printf("Meridian scale (h) : %.8f ( %.4g %% error )\n", @@ -534,6 +539,16 @@ int main(int argc, char **argv) { if (PJ *P = proj_create(nullptr, ocrs.c_str())) { const auto type = proj_get_type(P); if (type == PJ_TYPE_PROJECTED_CRS) { + try { + using namespace osgeo::proj; + auto crs = dynamic_cast( + P->iso_obj.get()); + auto dir = + crs->coordinateSystem()->axisList()[0]->direction(); + swapAxisCrs = dir == cs::AxisDirection::NORTH || + dir == cs::AxisDirection::SOUTH; + } catch (...) { + } auto ctx = P->ctx; auto geodetic_crs = proj_get_source_crs(ctx, P); assert(geodetic_crs); From 168a0e641312215b26845f45c954b60f84de6f2d Mon Sep 17 00:00:00 2001 From: Javier Jimenez Shaw Date: Tue, 25 Jul 2023 15:06:57 +0200 Subject: [PATCH 7/7] add CRS example in proj documentation --- docs/source/apps/proj.rst | 37 ++++++++++++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/docs/source/apps/proj.rst b/docs/source/apps/proj.rst index 2183bd0d20..327cc393ce 100644 --- a/docs/source/apps/proj.rst +++ b/docs/source/apps/proj.rst @@ -174,7 +174,7 @@ also used for supporting files like datum shift files. .. versionadded:: 9.3.0 - *{crs}* is one of the possibilities accepted by proj_create(), provided it + *{crs}* is one of the possibilities accepted by :c::func:`proj_create()`, provided it expresses a projected CRS, like a WKT string, an object code (like "EPSG:32632") a PROJJSON string, etc. The projection computed will be those of the map projection implied by @@ -212,6 +212,41 @@ data will appear as three lines of:: 460770.43 5011865.86 +This other example + +.. code-block:: console + + proj EPSG:6421 -V <