Consistency of crs.is_geographic between different EPSG representations

[sfinkens]

When creating an EPSG:4326 projection using a proj4 string, the value of the proj.crs.is_geographic property depends on whether definition_string() or crs.to_proj4(4) was used for initialization.

In [1]: import pyproj                                                                                                                                                                                                                        

In [2]: pyproj.Proj('EPSG:4326').crs.is_geographic                                                                                                                                                                                           
Out[2]: True

In [4]: pyproj.Proj('EPSG:4326').crs.to_proj4(4)                                                                                                                                                                                             
Out[4]: '+proj=longlat +datum=WGS84 +no_defs +type=crs'

In [5]: pyproj.Proj('EPSG:4326').definition_string()                                                                                                                                                                                         
Out[5]: 'proj=longlat datum=WGS84 no_defs ellps=WGS84 towgs84=0,0,0'

In [6]: pyproj.Proj(pyproj.Proj('EPSG:4326').crs.to_proj4(4)).crs.is_geographic                                                                                                                                                              
Out[6]: True

In [7]: pyproj.Proj(pyproj.Proj('EPSG:4326').definition_string()).crs.is_geographic                                                                                                                                                          
Out[7]: False

It seems like the extra towgs84=0,0,0 parameter in the definition_string causes the problem. Shouldn't this be consistent? Or is this just another example of losing information when converting CRSes to proj4 strings?

pyproj version is 2.1.3.

[jorisvandenbossche]

I have also been running into this, and it is due to the fact that, when a towgs84 entry is included in the proj4 string, this is interpreted as a BoundedCRS.

You can eg see this when displaying as WKT:

In [13]: crs = pyproj.Proj(pyproj.Proj('EPSG:4326').definition_string()).crs                                                                                  

In [14]: print(crs.to_wkt(pretty=True))                                                                                                                       
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["unknown",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]],
                ID["EPSG",6326]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8901]],
            CS[ellipsoidal,2],
                AXIS["longitude",east,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]],
                AXIS["latitude",north,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]]]],
    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["Geocentric translations (geog2D domain)",
            ID["EPSG",9603]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]]]]

[snowman2]

With changes in #263

>>> import pyproj
>>> crs = pyproj.Proj(pyproj.Proj('EPSG:4326').definition_string()).crs
>>> crs.is_geographic
False
>>> crs.source_crs.is_geographic
True

[snowman2]

The logic I have follows: https://github.com/OSGeo/gdal/blob/60ab2831abf9446c629b3a417a16cb35e7f41e63/gdal/ogr/ogrspatialreference.cpp#L7671-L7703

But, I wonder if there should also be a check if the original projection is bound and return geographic if it is bound and not a compound CRS. @rouault @kbevers

[snowman2]

Proposed solution: d69af25

[snowman2]

Actually, that appears to be the way to go. I see d->demoteFromBoundCRS(); in the code referenced above. So, the solution added should fix this. Thanks for raising this issue!

[sfinkens]

@snowman2 Thanks for looking into it!

[snowman2]

Not a problem 👍. The solution is in master now if you want to give it a try.