Adds CRS.geodetic_crs property (#3218)

* Adds CRS.to_geodetic function

* inspired by pyproj's .geodetic_crs function to return the
  corresponding base geographic CRS for a given CRS
* helpful for projections that may not use common geographic CRSs

* switch to else rather than finally in to_geodetic

* improve docstring and remove todos

* fixed test comment

* added .geodetic_crs property following pyproj's approach

* change to throw exception, added an aux.xml to .gitignore I keep almost committing

* Update rasterio/crs.pyx

Co-authored-by: Sean Gillies <sean.gillies@gmail.com>

* Update rasterio/crs.pyx

Co-authored-by: Sean Gillies <sean.gillies@gmail.com>

---------

Co-authored-by: Sean Gillies <sean.gillies@gmail.com>

.gitignore

@@ -73,6 +73,7 @@ VERSION.txt
 rasterio/*.cpp
 rasterio/*.c
 tests/data/white-gemini-iv.zip
+tests/data/all-nodata.tif.aux.xml
 .hypothesis/
 
 # vim

rasterio/crs.pxd

@@ -7,3 +7,4 @@ cdef class CRS:
     cdef object _data
     cdef public object _epsg
     cdef public object _wkt
+    cdef object _geodetic_crs

rasterio/crs.pyx

@@ -98,6 +98,7 @@ cdef class CRS:
         self._data = {}
         self._epsg = None
         self._wkt = None
+        self._geodetic_crs = None
 
         if initialdata or kwargs:
             tmp = CRS.from_dict(initialdata=initialdata, **kwargs)
@@ -272,6 +273,32 @@ cdef class CRS:
             units_b = units_c
             return (units_b.decode('utf-8'), factor)
 
+    @property
+    def geodetic_crs(self):
+        """Get the Geographic CRS from the CRS.
+
+        Returns
+        -------
+        CRS
+
+        Raises
+        ------
+        CRSError
+
+        """
+        if self._geodetic_crs:
+            return self._geodetic_crs
+        cdef CRS obj = CRS.__new__(CRS)
+        try:
+            obj._osr = exc_wrap_pointer(OSRCloneGeogCS(self._osr))
+        except CPLE_BaseError as exc:
+            raise CRSError("Cannot determine Geodetic CRS. {}".format(exc))
+        else:
+            osr_set_traditional_axis_mapping_strategy(obj._osr)
+            self._geodetic_crs = obj
+        return self._geodetic_crs 
+        
+
     def to_dict(self, projjson=False):
         """Convert CRS to a PROJ dict.
 

rasterio/gdal.pxi

@@ -166,6 +166,7 @@ cdef extern from "ogr_srs_api.h" nogil:
                      double *y, double *z)
     void OSRCleanup()
     OGRSpatialReferenceH OSRClone(OGRSpatialReferenceH srs)
+    OGRSpatialReferenceH OSRCloneGeogCS(OGRSpatialReferenceH srs)
     int OSRExportToProj4(OGRSpatialReferenceH srs, char **params)
     int OSRExportToWkt(OGRSpatialReferenceH srs, char **params)
     const char *OSRGetAuthorityName(OGRSpatialReferenceH srs, const char *key)

tests/test_crs.py

@@ -672,7 +672,24 @@ def test_crs_compound_epsg():
     assert CRS.from_string("EPSG:4326+3855").to_wkt().startswith("COMPD")
 
 
+@pytest.mark.parametrize(
+    'crs_obj,geod_crs',
+    [
+        (CRS.from_epsg(4087), CRS.from_epsg(4326)),
+        (CRS.from_epsg(3995), CRS.from_epsg(4326)),
+        (CRS.from_epsg(3031), CRS.from_epsg(4326)),
+        (CRS.from_user_input("ESRI:102004"), CRS.from_user_input("EPSG:4269")),
+        (CRS.from_user_input("IAU_2015:49910"), CRS.from_user_input("IAU_2015:49900")),
+        (CRS.from_user_input("IAU_2015:49911"), CRS.from_user_input("IAU_2015:49901"))
+    ]
+)
+def test_construct_geodetic_crs(crs_obj, geod_crs):
+    """Test if CRS geodetic CRS matches expectations."""
+    assert crs_obj.geodetic_crs == geod_crs
+
+    
 @pytest.mark.parametrize("crs", [CRS.from_epsg(4326), CRS.from_string("EPSG:4326")])
 def test_epsg_4326_ogc_crs84(crs):
     """EPSG:4326 not equivalent to OGC:CRS84."""
     assert CRS.from_string("OGC:CRS84") != crs
+