Bound CRS does not roundtrip through WKT

[brendan-ward]

Code Sample, a copy-pastable example if possible

Using the Bound CRS example from the geopandas docs:

>>> import pyproj
>>> crs = pyproj.CRS("+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m +no_defs")
>>> crs.to_wkt()
'BOUNDCRS[SOURCECRS[PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on International 1909 (Hayford) ellipsoid",ELLIPSOID["International 1909 (Hayford)",6378388,297,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",90,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",4.36748666666667,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",51.1666672333333,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",49.8333339,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",150000.013,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",5400088.438,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]],TARGETCRS[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["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]],ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",METHOD["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",106.869,ID["EPSG",8605]],PARAMETER["Y-axis translation",-52.2978,ID["EPSG",8606]],PARAMETER["Z-axis translation",103.724,ID["EPSG",8607]],PARAMETER["X-axis rotation",-0.33657,ID["EPSG",8608]],PARAMETER["Y-axis rotation",0.456955,ID["EPSG",8609]],PARAMETER["Z-axis rotation",-1.84218,ID["EPSG",8610]],PARAMETER["Scale difference",1.000001,ID["EPSG",8611]]]]'

Converting through WKT does not return the original CRS:

>>> crs2 = pyproj.CRS(crs.to_wkt())
>>> crs2 == crs
False
>>> crs.to_wkt() == crs2.to_wkt()
False
>>> crs2.to_wkt()
'BOUNDCRS[SOURCECRS[PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on International 1909 (Hayford) ellipsoid",ELLIPSOID["International 1909 (Hayford)",6378388,297,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",90,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",4.36748666666667,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",51.1666672333333,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",49.8333339,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",150000.013,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",5400088.438,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]],TARGETCRS[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]],ID["EPSG",4326]]],ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",METHOD["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",106.869,ID["EPSG",8605]],PARAMETER["Y-axis translation",-52.2978,ID["EPSG",8606]],PARAMETER["Z-axis translation",103.724,ID["EPSG",8607]],PARAMETER["X-axis rotation",-0.33657,ID["EPSG",8608]],PARAMETER["Y-axis rotation",0.456955,ID["EPSG",8609]],PARAMETER["Z-axis rotation",-1.84218,ID["EPSG",8610]],PARAMETER["Scale difference",1.000001,ID["EPSG",8611]]]]'

The WKT derived from each is a bit different:

>>> len(crs.to_wkt())
2020
>>> len(crs2.to_wkt())
2050

Expected Output

Bound CRS information should roundtrip through WKT, right? I didn't see anything in the docs to suggest that this should be lossy.

Environment Information

pyproj info:
    pyproj: 2.6.0
      PROJ: 6.3.1
  data dir: /Users/bcward/opt/miniconda3/envs/geopandas/lib/python3.8/site-packages/pyproj-2.6.0-py3.8-macosx-10.9-x86_64.egg/pyproj/proj_dir/share/proj

System:
    python: 3.8.2 (default, Mar 26 2020, 10:43:30)  [Clang 4.0.1 (tags/RELEASE_401/final)]
executable: /Users/bcward/opt/miniconda3/envs/geopandas/bin/python
   machine: macOS-10.15.3-x86_64-i386-64bit

Installation method

conda

[snowman2]

It shouldn't be lossy. Not sure why the two are not identified as equivalent.

I did notice this minor difference though.

Original CRS:

    TARGETCRS[
        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["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],

CRS from WKT:

    TARGETCRS[
        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]],
            ID["EPSG",4326]]],

Looks like the difference is in the name of the AXIS.

Example:
latitude -> geodetic latitude (Lat)

[snowman2]

Yep, that appears to be the only difference:

>>> from pyproj import CRS
>>> from difflib import ndiff
>>> crs = CRS("+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m +no_defs")
>>> print("".join(ndiff(crs.to_wkt(pretty=True).splitlines(1), CRS(crs).to_wkt(pretty=True).splitlines(1))))
  BOUNDCRS[
      SOURCECRS[
          PROJCRS["unknown",
              BASEGEOGCRS["unknown",
                  DATUM["Unknown based on International 1909 (Hayford) ellipsoid",
                      ELLIPSOID["International 1909 (Hayford)",6378388,297,
                          LENGTHUNIT["metre",1,
                              ID["EPSG",9001]]]],
                  PRIMEM["Greenwich",0,
                      ANGLEUNIT["degree",0.0174532925199433],
                      ID["EPSG",8901]]],
              CONVERSION["unknown",
                  METHOD["Lambert Conic Conformal (2SP)",
                      ID["EPSG",9802]],
                  PARAMETER["Latitude of false origin",90,
                      ANGLEUNIT["degree",0.0174532925199433],
                      ID["EPSG",8821]],
                  PARAMETER["Longitude of false origin",4.36748666666667,
                      ANGLEUNIT["degree",0.0174532925199433],
                      ID["EPSG",8822]],
                  PARAMETER["Latitude of 1st standard parallel",51.1666672333333,
                      ANGLEUNIT["degree",0.0174532925199433],
                      ID["EPSG",8823]],
                  PARAMETER["Latitude of 2nd standard parallel",49.8333339,
                      ANGLEUNIT["degree",0.0174532925199433],
                      ID["EPSG",8824]],
                  PARAMETER["Easting at false origin",150000.013,
                      LENGTHUNIT["metre",1],
                      ID["EPSG",8826]],
                  PARAMETER["Northing at false origin",5400088.438,
                      LENGTHUNIT["metre",1],
                      ID["EPSG",8827]]],
              CS[Cartesian,2],
                  AXIS["(E)",east,
                      ORDER[1],
                      LENGTHUNIT["metre",1,
                          ID["EPSG",9001]]],
                  AXIS["(N)",north,
                      ORDER[2],
                      LENGTHUNIT["metre",1,
                          ID["EPSG",9001]]]]],
      TARGETCRS[
          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["latitude",north,
+                 AXIS["geodetic latitude (Lat)",north,
?                       +++++++++        ++++++
                      ORDER[1],
                      ANGLEUNIT["degree",0.0174532925199433]],
-                 AXIS["longitude",east,
+                 AXIS["geodetic longitude (Lon)",east,
?                       +++++++++         ++++++
                      ORDER[2],
                      ANGLEUNIT["degree",0.0174532925199433]],
              ID["EPSG",4326]]],
      ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
          METHOD["Position Vector transformation (geog2D domain)",
              ID["EPSG",9606]],
          PARAMETER["X-axis translation",106.869,
              ID["EPSG",8605]],
          PARAMETER["Y-axis translation",-52.2978,
              ID["EPSG",8606]],
          PARAMETER["Z-axis translation",103.724,
              ID["EPSG",8607]],
          PARAMETER["X-axis rotation",-0.33657,
              ID["EPSG",8608]],
          PARAMETER["Y-axis rotation",0.456955,
              ID["EPSG",8609]],
          PARAMETER["Z-axis rotation",-1.84218,
              ID["EPSG",8610]],
          PARAMETER["Scale difference",1.000001,
              ID["EPSG",8611]]]]

[snowman2]

This is interesting:

>>> crs == CRS(crs)
False
>>> crs.target_crs == CRS(crs).target_crs
True
>>> crs.source_crs == CRS(crs).source_crs
True
>>> crs.coordinate_operation == CRS(crs).coordinate_operation
True

$ python -m pyproj -v
pyproj info:
    pyproj: 2.6.0
      PROJ: 7.0.0
  data dir: /home/snowal/miniconda/envs/proj700/share/proj
...

[snowman2]

Looks like once you convert to WKT it compares just fine:

>>> from pyproj import CRS
>>> crs = CRS("+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m +no_defs")
>>> crswkt = CRS(crs)
>>> crswkt == CRS(crswkt)
True

[snowman2]

Forwarded to: OSGeo/PROJ#2149

[brendan-ward]

Thanks for investigating @snowman2 !