Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add option in proj CLI to use a CRS #3825

Merged
merged 7 commits into from
Jul 26, 2023
Merged

Add option in proj CLI to use a CRS #3825

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
1ac2c96
Add option in proj CLI to use a CRS
jjimenezshaw Jul 16, 2023
3dfbed7
add test for proj cli with -C option
jjimenezshaw Jul 16, 2023
05ccb72
choose more numerically stable test case
jjimenezshaw Jul 16, 2023
9adf4b5
Read CRS directly from the arguments, without option flag.
jjimenezshaw Jul 17, 2023
5301fdc
pass metre and e-n axis order to proj_factors
jjimenezshaw Jul 20, 2023
807ae5a
consider northing-easing in verbose mode
jjimenezshaw Jul 20, 2023
168a0e6
add CRS example in proj documentation
jjimenezshaw Jul 25, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 44 additions & 2 deletions docs/source/apps/proj.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
103 changes: 88 additions & 15 deletions src/apps/proj.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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>
Expand All @@ -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__)
Expand All @@ -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 *);
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
}
Expand All @@ -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",
Expand Down Expand Up @@ -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)) {
Expand All @@ -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 "
Expand Down Expand Up @@ -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);

Expand Down
6 changes: 6 additions & 0 deletions test/cli/testproj
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions test/cli/testproj_out.dist
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
2 49 2.000 49.000
500000.00 0.00 <0.9996 0.9996 0.9992 0 0.9996 0.9996>
Loading