Add option in proj CLI to use a CRS (OSGeo#3825)

jjimenezshaw · web-flow · commit fde3d58e4af5 · 2023-07-26T09:20:27.000+02:00
To easily apply proj on the projection of a projected CRS, this PR adds the option `-C`.
It follows the same logic as in function proj_factors, that is to compute the projection from the CRS to the geographic underlying CRS.

Example of usage:

```
echo -110.16666 31 | proj -C EPSG:26948 -V 
#Transverse Mercator
#	Cyl, Sph&amp;Ell
#	approx
# +proj=tmerc +lat_0=31 +lon_0=-110.166666666667 +k=0.9999 +x_0=213360 +y_0=0
# +ellps=GRS80
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257222
#  squared eccentricity: 0.006694380023
Longitude: 110d9'59.976"W [ -110.16666 ]
Latitude:  31dN [ 31 ]
Easting (x):   213360.64
Northing (y):  0.00
Meridian scale (h) : 0.99990000  ( -0.01 % error )
Parallel scale (k) : 0.99990000  ( -0.01 % error )
Areal scale (s):     0.99980001  ( -0.02 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 0d0'0.012" [ 0.00000343 ]
Max-min (Tissot axis a-b) scale error: 0.99990 0.99990
```
diff --git 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** [**-beEfiIlmorsStTvVwW**] [args]] ([*+opt[=arg]* ...] | {crs}) file ...
 
-    **invproj** [**-beEfiIlmorsStTvVwW**] [args]] [*+opt[=arg]* ...] file ...
+    **invproj** [**-beEfiIlmorsStTvVwW**] [args]] ([*+opt[=arg]* ...] | {crs}) file ...
 
 
 Description
@@ -172,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 <projections>`.
 
+.. versionadded:: 9.3.0
+
+    *{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
+    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
@@ -205,6 +212,41 @@ data will appear as three lines of::
 
     460770.43     5011865.86
 
+This other example
+
+.. code-block:: console
+
+    proj EPSG:6421 -V <<EOF
+    -120 35.8
+    EOF
+
+Will perform the projection of the coordinates in "NAD83(2011) / California zone 4"
+(`EPSG:6421`) into its geographic system, "NAD83(2011)", showing the expanded annotated listing.
+The output will appear as:
+
+.. code-block:: console
+
+    #Lambert Conformal Conic
+    #	Conic, Sph&Ell
+    #	lat_1= and lat_2= or lat_0, k_0=
+    # +proj=lcc +lat_0=35.3333333333333 +lon_0=-119 +lat_1=37.25 +lat_2=36
+    # +x_0=2000000 +y_0=500000 +ellps=GRS80
+    #Final Earth figure: ellipsoid
+    #  Major axis (a): 6378137.000
+    #  1/flattening: 298.257222
+    #  squared eccentricity: 0.006694380023
+    Longitude: 120dW [ -120 ]
+    Latitude:  35d48'N [ 35.8 ]
+    Easting (x):   1909606.87
+    Northing (y):  552253.58
+    Meridian scale (h) : 1.00004382  ( 0.004382 % error )
+    Parallel scale (k) : 1.00004382  ( 0.004382 % error )
+    Areal scale (s):     1.00008765  ( 0.008765 % error )
+    Angular distortion (w): 0.000
+    Meridian/Parallel angle: 90.00000
+    Convergence : -0d35'47.714" [ -0.59658715 ]
+    Max-min (Tissot axis a-b) scale error: 1.00004 1.00004
+
 .. only:: man
 
     Other programs
diff --git 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 <ctype.h>
@@ -9,6 +10,8 @@
 #include <stdlib.h>
 #include <string.h>
 
+#include <proj/crs.hpp>
+
 #include <vector>
 
 #if defined(MSDOS) || defined(OS2) || defined(WIN32) || defined(__WIN32__)
@@ -22,7 +25,9 @@
 #define MAX_LINE 1000
 #define PJ_INVERSE(P) (P->inv ? 1 : 0)
 
-static PJ *Proj;
+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 *);
@@ -115,16 +120,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) {
@@ -281,8 +286,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;
         }
@@ -298,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",
@@ -507,8 +514,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<char *>("-");
 
     if (oform) {
         if (!validate_form_string_for_numbers(oform)) {
@@ -525,14 +530,80 @@ int main(int argc, char **argv) {
     }
     proj_context_use_proj4_init_rules(nullptr, true);
 
+    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.c_str())) {
+            const auto type = proj_get_type(P);
+            if (type == PJ_TYPE_PROJECTED_CRS) {
+                try {
+                    using namespace osgeo::proj;
+                    auto crs = dynamic_cast<const crs::ProjectedCRS *>(
+                        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);
+                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 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_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");
+            }
+            proj_destroy(P);
+        } else {
+            emess(3, "CRS is not parseable");
+        }
+    }
+    if (eargc == 0) /* if no specific files force sysin */
+        eargv[eargc++] = const_cast<char *>("-");
+
     // proj historically ignores any datum shift specifier, like nadgrids,
     // towgs84, etc
     argvVector.push_back(const_cast<char *>("break_cs2cs_recursion"));
 
-    if (!(Proj = proj_create_argv(nullptr, static_cast<int>(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<int>(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 "
@@ -616,6 +687,8 @@ int main(int argc, char **argv) {
         emess_dat.File_name = nullptr;
     }
 
+    if (ProjForFactors && ProjForFactors != Proj)
+        proj_destroy(ProjForFactors);
     if (Proj)
         proj_destroy(Proj);
 
diff --git 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} <<EOF
 2 49
 EOF
+#
+echo "Test CRS option"
+#
+$EXE EPSG:32620 -S >>${OUT} <<EOF
+-63 0
+EOF
 
 #
 # do 'diff' with distribution results
diff --git a/test/cli/testproj_out.dist b/test/cli/testproj_out.dist
@@ -1 +1,2 @@
 2 49	2.000	49.000
+500000.00	0.00	<0.9996 0.9996 0.9992 0 0.9996 0.9996>

Original file line number	Diff line number	Diff line change
`@@ -1 +1,2 @@`
`1`	`1`	`2 49 2.000 49.000`
	`2`	`+500000.00 0.00 <0.9996 0.9996 0.9992 0 0.9996 0.9996>`