Update docstrings

[mullenkamp]

Hi,

Could you please update the docstring documentation for the Transformer class?
One minor typo is proj_from in the from_proj method, but the most important ones are the itransform and transform methods.

They both state that the input should be in the order of x, y when they actually are in the order of y, x. And the output is in the order of y, x. That got me quite confused for a while.

Thanks

[snowman2]

Huh, it is x,y,z. I think it is because it is swapped for some EPSG projections. @rouault @kbevers can you verify my assumption?

[snowman2]

@mullenkamp, would you mind providing specific code examples?

[snowman2]

Reference to my assumption reasoning:OSGeo/PROJ#1301

[mullenkamp]

No problem.

from pyproj import CRS, Proj, Transformer

epsg1 = 4326
epsg2 = 2193

x1 = 173.2996473
y1 = -40.60674803

x2 = 1625350
y2 = 5504853

crs1 = CRS.from_epsg(epsg1)
crs2 = CRS.from_epsg(epsg2)

trans1 = Transformer.from_crs(crs1, crs2)

set1a = trans1.transform(y1, x1)

trans2 = Transformer.from_crs(crs2, crs1)

set2a = trans2.transform(y2, x2)

These produce the desired output. Reverse them to see strange output (i.e. inf or incorrect coords).

Thanks!

[kbevers]

This entirely depend on how the source and destination CRS's are defined. By chance, both your source and destination are defined as having the first coordinate component point in a northerly direction. A general rule of thumb is that a CRS in geographic coordinates has an axis ordering of (latitude, longitude) where as most projected CRS's has axis on ordering of (easting, northing).

projinfo can help you understand the CRS's you are working with:

$ ./bin/projinfo EPSG:4326
PROJ.4 string:
+proj=longlat +datum=WGS84 +no_defs +type=crs

WKT2_2018 string:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

$ ./bin/projinfo EPSG:2193
PROJ.4 string:
+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs +type=crs

WKT2_2018 string:
PROJCRS["NZGD2000 / New Zealand Transverse Mercator 2000",
    BASEGEOGCRS["NZGD2000",
        DATUM["New Zealand Geodetic Datum 2000",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4167]],
    CONVERSION["New Zealand Transverse Mercator 2000",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",173,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",1600000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["New Zealand - onshore"],
        BBOX[-47.33,166.37,-34.1,178.63]],
    ID["EPSG",2193]]

Note the CS section where AXIS comes in the order they are defined by the CRS authority.

[mullenkamp]

Thank you very much for the explanation. I do understand how this works now.

But unfortunately that still does mean that the x and y order in pyproj is totally dependent on the crs definition rather than the class and method itself.
Is there any possibility or interest to try to standardize the ordering for pyproj?
Or should it be expected that we get that info first from the axis_info attribute then order the x and y appropriately?
Either way then docstrings will still need some updating.

I'm happy to help in anyway I can by the way.

[snowman2]

This is definitely an issue that will need some thought and careful consideration. I am not sure if this should be handled in pyproj or in PROJ.

But, as you mentioned, this should definitely be documented at a minimum.

If you have ideas for handling this, feel free to comment below. My main idea is to use the axis_info property on the CRS object to check for this special case.

[kbevers]

This is definitely an issue that will need some thought and careful consideration. I am not sure if this should be handled in pyproj or in PROJ.

First of all, I would say that pyproj should make a principal decision as a project on how it wants to treat coordinates. Respect the axis order of local authorities or go with "east is always first" convention that is used by most GIS tools. We are working on providing better tools in PROJ for supporting the latter. I expect that version 6.1 includes a simpler solution than what is currently available.

[snowman2]

@kbevers, I would think pyproj would want to follow how PROJ is doing it by default. Maybe add a kwarg for strict_axis_order for users who want to keep the axis order east/west. @jswhit, thoughts on this one?

[snowman2]

I expect that version 6.1 includes a simpler solution than what is currently available.

That will be useful.

[mullenkamp]

Can anyone think of a particular use-case when someone is building upon pyproj for a universal crs transformer in their application that they would want the axis order to vary?
If you can, then we should build the option into pyproj as described by @snowman2 . But I'd say make the default in a fixed standardized order (i.e. x, y) and provide the option to order it according to the PROJ.
For me, building tools on top of pyproj (which is what I'm doing/fixing now), I'd want it consistent. Otherwise I'd have to look first at the axis_info attribute for every tool that I'd build, which would make me want to build a wrapper on top of pyproj so that I wouldn't be duplicating code all the time.

[kbevers]

Can anyone think of a particular use-case when someone is building upon pyproj for a universal crs transformer in their application that they would want the axis order to vary?

Anyone who wants to represent a given CRS as the local authority intended. Or really, if you want to claim that something is CRS xyz then you should represent it as CRS xyz is defined, and not how you think it should be defined. This can seem as useless nitpicking on the surface, but once you dig a bit deeper the real problems begin to show up. For most coordinates systems it makes much sense to put the east/west-component first because that looks like the x-axis in a normal coordinate system. The same for the north/south component that in most cases fit very well with the y-axis. But then we have coordinate systems that are rotated and have no obvious mapping to the usual x and y direction - what do you do then? Or geocentric cartesian coordinate systems which has an x, y and z axis where the z-axis points towards north?

And then we have sailors (and many others) that will always give you a position as latitude then longitude. PROJ and pyproj respects that now, the previous behaviour did not.

Personally I would really love to live in a world where this stuff was simple, but it isn't - sometimes for good reasons - and the best we can do is try to adapt and deal with it as best as possible. This is mostly a problem now because GIS was mostly invented by computer programmers and not geodesists and geographers. The geodesist and geographers are slowly catching up and this is what is now stirring the pond. Sorry for the inconvenience - everything will work out eventually :-)

[mullenkamp]

Can anyone think of a particular use-case when someone is building upon pyproj for a universal crs transformer in their application that they would want the axis order to vary?

Anyone who wants to represent a given CRS as the local authority intended. Or really, if you want to claim that something is CRS xyz then you should represent it as CRS xyz is defined, and not how you think it should be defined. This can seem as useless nitpicking on the surface, but once you dig a bit deeper the real problems begin to show up. For most coordinates systems it makes much sense to put the east/west-component first because that looks like the x-axis in a normal coordinate system. The same for the north/south component that in most cases fit very well with the y-axis. But then we have coordinate systems that are rotated and have no obvious mapping to the usual x and y direction - what do you do then? Or geocentric cartesian coordinate systems which has an x, y and z axis where the z-axis points towards north?

And then we have sailors (and many others) that will always give you a position as latitude then longitude. PROJ and pyproj respects that now, the previous behaviour did not.

Personally I would really love to live in a world where this stuff was simple, but it isn't - sometimes for good reasons - and the best we can do is try to adapt and deal with it as best as possible. This is mostly a problem now because GIS was mostly invented by computer programmers and not geodesists and geographers. The geodesist and geographers are slowly catching up and this is what is now stirring the pond. Sorry for the inconvenience - everything will work out eventually :-)

Thanks for the explanation on the background of CRS and proj4.
Sounds like we should provide the option for either if there is still appetite for it. It shouldn't be too hard and I'm happy to provide the functionality since I need to do it anyway.
How does that sound?

[snowman2]

If you would like to take a stab at it and make a PR to discuss it, that's fine with me 👍

[mullenkamp]

Thanks @snowman2 .
Oh and I'll update the documentation while I'm at it ;)

[snowman2]

Docs already in progress: #229

[snowman2]

@mullenkamp it seems PROJ already found a solution. It will be included in version 6.1.0. See #238.

[mullenkamp]

Does this mean I shouldn't make any changes, or try to implement version 6.1.0 into pyproj?

[snowman2]

If you would like to implement it from the 6.1.0 version, that would be great.

[snowman2]

@mullenkamp, I plan to start an implementation on this in the next couple of days so it can be in the next release.

[snowman2]

I have an implementation in #286.

[mullenkamp]

Thanks! Sorry that I haven't had time recently to make the implementation myself.

[snowman2]

No worries. I wasn't sure when you were planning on getting to it and was thinking it would be an important feature to include in the next release (that I thought was going to happen 3 days ago when I started on it :)).

[snowman2]

Closing as this has been merged into master.