WKT1 format for CRS84 from GDAL results in EPSG:4326 with lon/lat axis order

[jorisvandenbossche]

From Toblerity/Fiona#816 (comment), I noticed that the WKT that was produced by Fiona / GDAL for the "urn:ogc:def:crs:OGC:1.3:CRS84" has a strange behaviour when read in pyproj:

In [11]: crs = pyproj.CRS('GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]')                                     

In [12]: crs                                                                                                                                                                                                       
Out[12]: 
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- lon[east]: Longitude (Degree)
- lat[north]: Latitude (Degree)
Area of Use:
- undefined
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

In [13]: crs.to_authority()                                                                                                                                                                                        
Out[13]: ('EPSG', '4326')

So you can see that the axis order is correct I think (lon, lat), as this is coming from CRS84 (which should be EPSG:4326 but with swapped axis). But it still directly recognizes this as EPSG:4326 (it's even in the first line of the repr).
Which seemingly gives a CRS that contradicts itself? (this probably is a behaviour defined by PROJ though)

---

In [16]: pyproj.show_versions()   
pyproj info:
    pyproj: 2.4.0
      PROJ: 6.2.0
  data dir: /home/joris/miniconda3/envs/test-fiona-189/share/proj

System:
    python: 3.7.3 | packaged by conda-forge | (default, Jul  1 2019, 21:52:21)  [GCC 7.3.0]
executable: /home/joris/miniconda3/envs/test-fiona-189/bin/python
   machine: Linux-4.15.0-1059-oem-x86_64-with-debian-buster-sid

Python deps:
       pip: 19.3.1
setuptools: 41.6.0.post20191101
    Cython: None

[snowman2]

Interesting.

this probably is a behaviour defined by PROJ though

Yes, that is the case. Probably need to raise this issue upstream.

[jorisvandenbossche]

cc @rouault

[rouault]

Yes, this is a fuzzy area... "GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]" is the canonical representation of EPSG:4326 in ESRI WKT. So identification to it is OK. For some reason I don't really remember, I preferred that on import of such a WKT without an explicit axis order, we use traditional axis order "long", "lat". Not really sure what would be the impact of changing this

[jorisvandenbossche]

I think what I find most unexpected (without knowing the background you just sketched) is that it gives EPSG:4326 with full confidence, while the axis order does not match:

In [73]: crs.to_authority(min_confidence=100)  
Out[73]: ('EPSG', '4326')

While for a similar object (reconstructed by roundtripping to WKT2), needs a min_confidence of 20 to ignore the axis order to return EPSG:4326:

In [74]: crs2 = pyproj.CRS(crs.to_wkt()) 

In [75]: crs2.to_authority(min_confidence=100) 

In [76]: crs2.to_authority(min_confidence=70)
Out[76]: ('OGC', 'CRS84')

In [77]: crs2.to_epsg(min_confidence=20)   
Out[77]: 4326

[jorisvandenbossche]

But thanks for the explanation of the context!

In any case, it's thus not a pyproj issue, so closing here.