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.

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

Add option in proj CLI to use a CRS #3825

Merged
merged 7 commits into from
Jul 26, 2023
Merged

Conversation

jjimenezshaw
Copy link
Contributor

@jjimenezshaw jjimenezshaw commented Jul 16, 2023

Following #3824 , 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.

This PR does NOT fix #3824.

Example of usage:

echo -110.16666 31 | proj -C EPSG:26948 -V 
#Transverse Mercator
#	Cyl, Sph&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
  • Tests added
  • Added clear title that can be used to generate release notes
  • Fully documented, including updating docs/source/*.rst for new API

@jmckenna
Copy link
Contributor

thanks @jjimenezshaw very useful!

@jjimenezshaw
Copy link
Contributor Author

BTW, I chose -C just due to "CRS". If there is a more consistent option among the CLI programs, I can use a different one. Just let me know.

@busstoptaktik
Copy link
Member

Perhaps skip the flag and just parse whether it looks like a CRS, rather than a PROJ definition?

@rouault
Copy link
Member

rouault commented Jul 17, 2023

Perhaps skip the flag and just parse whether it looks like a CRS, rather than a PROJ definition?

I was just going to suggest the same thing. cs2cs has logic to be able to deal with things like "cs2cs EPSG:foo EPSG:bar" or "cs2cs +proj=... +to + proj=...." automatically.

Here I believe that the criterion would be "argvVector.empty() && eargc >= 1", and consider eargv[0] as the CRS code instead of being the first file

@rouault
Copy link
Member

rouault commented Jul 17, 2023

Hum while fixing #3824, I'm wondering what we should do if allowing to pass CRS by codes when they have northing/easting axis order:

  • should we normalize to easting/northing order ? That way in verbose mode, the "Easting (x)" and "Northing (y)" labels will be OK
  • or should we output with axis axis order of the CRS ? If so, then we should switch the values in verbose mode.
    At the very least this will deserve documentation.

Note that the proj utility didn't honour the (now deprecated) +axis= modifier which was a cs2cs / pj_transform() thing and ignored by the pj_fwd() logic .

@jjimenezshaw
Copy link
Contributor Author

what we should do if allowing to pass CRS by codes when they have northing/easting axis order

I was making myself the same question. Somehow I have a preference to keep the CRS order (to be consistent with the printed pipeline). But it is not very strong opinion.

See that the problem is not only with CRS codes (or wkt). If you pass the pipeline with the units or axis order (or anything that is not the "standar" way probably) #3824 still fails (probably you realized that already).

@jjimenezshaw
Copy link
Contributor Author

In case it is not clear: the CRS definition can be anything that works with "from user input". Can be a code, a WKT, a PROJJSON, etc.
I want to study later how it works with a derived projected (enabling that case in the if statement).

@jjimenezshaw jjimenezshaw force-pushed the crs-in-proj-cli branch 2 times, most recently from 60c2624 to 8ac1a2e Compare July 17, 2023 21:50
@jjimenezshaw
Copy link
Contributor Author

I applied the same criterion as in #3826 for the proj_factors calls (not for the projection calls).
I am still thinking about the axis order. Is there any function to easily know if the PJ* object has swapped axis?

@rouault
Copy link
Member

rouault commented Jul 20, 2023

Is there any function to easily know if the PJ* object has swapped axis?

With the C++ API: CRS::mustAxisOrderBeSwitchedForVisualization()

@jjimenezshaw
Copy link
Contributor Author

mustAxisOrderBeSwitchedForVisualization

Unfortunately the linker does not find it. Anyhow, the check was easy with the C++ object.

The output of a northing-easting CRS is something like this:

echo 9 0 | ./bin/proj EPSG:3044 -V
#Transformation pipeline manager
# +proj=pipeline +step +proj=utm +zone=32 +ellps=GRS80 +step +proj=axisswap
# +order=2,1
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257222
#  squared eccentricity: 0.006694380023
Longitude: 9dE [ 9 ]
Latitude:  0dN [ 0 ]
Northing (y):  0.00
Easting (x):   500000.00
Meridian scale (h) : 0.99960000  ( -0.04 % error )
Parallel scale (k) : 0.99960000  ( -0.04 % error )
Areal scale (s):     0.99920016  ( -0.07998 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 0d [ -0.00000000 ]
Max-min (Tissot axis a-b) scale error: 0.99960 0.99960

@jjimenezshaw
Copy link
Contributor Author

Is there anything missing in this PR? (behavior, feature or documentation)

Copy link
Member

@kbevers kbevers left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me, but a few minor things should be improved in the docs

@@ -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 proj_create(), provided it
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
*{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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Further down the page an example using proj with a CRS would be nice as well

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

@kbevers kbevers merged commit fde3d58 into OSGeo:master Jul 26, 2023
21 checks passed
@kbevers kbevers added this to the 9.3.0 milestone Jul 26, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

proj_factors produces wrong results with EPSG:2222 (ft)
5 participants