REF: Replace usage of _osr_from_crs with CRS.from_user_input (#2488)

rasterio/_base.pxd

@@ -32,7 +32,6 @@ cdef class DatasetBase:
 
 
 cdef const char *get_driver_name(GDALDriverH driver)
-cdef OGRSpatialReferenceH _osr_from_crs(object crs) except NULL
+
-cdef _safe_osr_release(OGRSpatialReferenceH srs)
 cdef void osr_set_traditional_axis_mapping_strategy(OGRSpatialReferenceH hSrs)
 cdef GDALDatasetH open_dataset(object filename, int mode, object allowed_drivers, object open_options, object siblings) except NULL

rasterio/_base.pyx

@@ -10,6 +10,7 @@ import os
 import warnings
 
 from libc.string cimport strncmp
+from rasterio.crs cimport CRS
 
 from rasterio._err import (
     GDALError, CPLE_BaseError, CPLE_IllegalArgError, CPLE_OpenFailedError,
@@ -1454,16 +1455,14 @@ def _transform(src_crs, dst_crs, xs, ys, zs):
     cdef double *x = NULL
     cdef double *y = NULL
     cdef double *z = NULL
-    cdef OGRSpatialReferenceH src = NULL
-    cdef OGRSpatialReferenceH dst = NULL
     cdef OGRCoordinateTransformationH transform = NULL
     cdef int i
 
     assert len(xs) == len(ys)
     assert zs is None or len(xs) == len(zs)
 
-    src = _osr_from_crs(src_crs)
+    cdef CRS src = CRS.from_user_input(src_crs)
-    dst = _osr_from_crs(dst_crs)
+    cdef CRS dst = CRS.from_user_input(dst_crs)
 
     n = len(xs)
     x = <double *>CPLMalloc(n*sizeof(double))
@@ -1478,7 +1477,7 @@ def _transform(src_crs, dst_crs, xs, ys, zs):
             z[i] = zs[i]
 
     try:
-        transform = OCTNewCoordinateTransformation(src, dst)
+        transform = OCTNewCoordinateTransformation(src._osr, dst._osr)
         transform = exc_wrap_pointer(transform)
         # OCTTransform() returns TRUE/FALSE contrary to most GDAL API functions
         exc_wrap_int(OCTTransform(transform, n, x, y, z) == 0)
@@ -1500,47 +1499,6 @@ def _transform(src_crs, dst_crs, xs, ys, zs):
         CPLFree(y)
         CPLFree(z)
         OCTDestroyCoordinateTransformation(transform)
-        _safe_osr_release(src)
-        _safe_osr_release(dst)
-
-
-cdef OGRSpatialReferenceH _osr_from_crs(object crs) except NULL:
-    """Returns a reference to memory that must be deallocated
-    by the caller."""
-    crs = CRS.from_user_input(crs)
-
-    # EPSG is a special case.
-    init = crs.get('init')
-    if init:
-        auth, val = init.strip().split(':')
-
-        if not val or auth.upper() != 'EPSG':
-            raise CRSError("Invalid CRS: {!r}".format(crs))
-        proj = 'EPSG:{}'.format(val).encode('utf-8')
-    else:
-        proj = crs.to_string().encode('utf-8')
-        log.debug("PROJ.4 to be imported: %r", proj)
-
-    cdef OGRSpatialReferenceH osr = OSRNewSpatialReference(NULL)
-    try:
-        retval = exc_wrap_int(OSRSetFromUserInput(osr, <const char *>proj))
-        if retval:
-            _safe_osr_release(osr)
-            raise CRSError("Invalid CRS: {!r}".format(crs))
-    except CPLE_BaseError as exc:
-        _safe_osr_release(osr)
-        raise CRSError(str(exc))
-    else:
-        osr_set_traditional_axis_mapping_strategy(osr)
-        return osr
-
-
-cdef _safe_osr_release(OGRSpatialReferenceH srs):
-    """Wrapper to handle OSR release when NULL."""
-
-    if srs != NULL:
-        OSRRelease(srs)
-    srs = NULL
 
 
 def _can_create_osr(crs):
@@ -1558,26 +1516,14 @@ def _can_create_osr(crs):
         True if source coordinate reference appears valid.
     """
 
-    cdef char *wkt = NULL
-    cdef OGRSpatialReferenceH osr = NULL
-
     try:
-        # Note: _osr_from_crs() has "except NULL" in its signature.
+        wkt = CRS.from_user_input(crs).to_wkt()
-        # It raises, it does not return NULL.
-        osr = _osr_from_crs(crs)
-        OSRExportToWkt(osr, &wkt)
-
         # If input was empty, WKT can be too; otherwise the conversion
         # didn't work properly and indicates an error.
-        return wkt != NULL and bool(crs) == (wkt[0] != 0)
+        return bool(wkt)
-
     except CRSError:
         return False
 
-    finally:
-        _safe_osr_release(osr)
-        CPLFree(wkt)
-
 
 cdef void osr_set_traditional_axis_mapping_strategy(OGRSpatialReferenceH hSrs):
     OSRSetAxisMappingStrategy(hSrs, OAMS_TRADITIONAL_GIS_ORDER)

rasterio/_env.pyx

@@ -17,7 +17,6 @@ import os.path
 import sys
 import threading
 
-from rasterio._base cimport _safe_osr_release
 from rasterio._err import CPLE_BaseError
 from rasterio._err cimport exc_wrap_ogrerr, exc_wrap_int
 from rasterio._filepath cimport install_filepath_plugin, uninstall_filepath_plugin
@@ -312,7 +311,8 @@ class PROJDataFinder:
         else:
             return True
         finally:
-            _safe_osr_release(osr)
+            if osr != NULL:
+                OSRRelease(osr)
 
 
     def search(self, prefix=None):

rasterio/_io.pyx

@@ -38,8 +38,7 @@ from libc.stdio cimport FILE
 from rasterio.enums import Resampling
 from rasterio.env import GDALVersion
 from rasterio.errors import ResamplingAlgorithmError, DatasetIOShapeError
-from rasterio._base cimport (
+from rasterio._base cimport get_driver_name, DatasetBase
-    _osr_from_crs, _safe_osr_release, get_driver_name, DatasetBase)
 from rasterio._err cimport exc_wrap_int, exc_wrap_pointer, exc_wrap_vsilfile
 
 cimport numpy as np
@@ -2001,15 +2000,15 @@ cdef class DatasetWriterBase(DatasetReaderBase):
                 gcplist[i].dfGCPZ = obj.z or 0.0
 
             # Try to use the primary crs if possible.
-            if not crs:
+            if crs is None:
                 crs = self.crs
-
+            else:
-            osr = _osr_from_crs(crs)
+                crs = CRS.from_user_input(crs)
-            OSRExportToWkt(osr, <char**>&srcwkt)
+            srcwkt_b = crs.to_wkt().encode('utf-8')
+            srcwkt = srcwkt_b
             GDALSetGCPs(self.handle(), len(gcps), gcplist, srcwkt)
         finally:
             CPLFree(gcplist)
-            CPLFree(srcwkt)
 
         # Invalidate cached value.
         self._gcps = None

rasterio/_warp.pyx

@@ -35,12 +35,12 @@ cimport cython
 cimport numpy as np
 from libc.math cimport HUGE_VAL
 
-from rasterio._base cimport _osr_from_crs, get_driver_name, _safe_osr_release
+from rasterio._base cimport get_driver_name
 from rasterio._err cimport exc_wrap, exc_wrap_pointer, exc_wrap_int
 from rasterio._io cimport (
     DatasetReaderBase, MemoryDataset, in_dtype_range, io_auto)
 from rasterio._features cimport GeomBuilder, OGRGeomBuilder
-
+from rasterio.crs cimport CRS
 
 log = logging.getLogger(__name__)
 
@@ -92,19 +92,13 @@ def _transform_geom(
         int precision):
     """Return a transformed geometry."""
     cdef char **options = NULL
-    cdef OGRSpatialReferenceH src = NULL
-    cdef OGRSpatialReferenceH dst = NULL
     cdef OGRCoordinateTransformationH transform = NULL
     cdef OGRGeometryFactory *factory = NULL
 
-    src = _osr_from_crs(src_crs)
+    cdef CRS src = CRS.from_user_input(src_crs)
-    dst = _osr_from_crs(dst_crs)
+    cdef CRS dst = CRS.from_user_input(dst_crs)
 
-    try:
+    transform = exc_wrap_pointer(OCTNewCoordinateTransformation(src._osr, dst._osr))
-        transform = exc_wrap_pointer(OCTNewCoordinateTransformation(src, dst))
-    finally:
-        _safe_osr_release(src)
-        _safe_osr_release(dst)
 
     # GDAL cuts on the antimeridian by default and using different
     # logic in versions >= 2.2.
@@ -632,7 +626,6 @@ def _calculate_default_transform(
     cdef int nlines = 0
     cdef double extent[4]
     cdef double geotransform[6]
-    cdef OGRSpatialReferenceH osr = NULL
     cdef char *wkt = NULL
     cdef GDALDatasetH hds = NULL
     cdef char **imgProjOptions = NULL
@@ -650,19 +643,11 @@ def _calculate_default_transform(
     else:
         transform = None
 
-    try:
+    dst_crs = CRS.from_user_input(dst_crs)
-        osr = _osr_from_crs(dst_crs)
+    wkt_b = dst_crs.to_wkt().encode('utf-8')
-        exc_wrap_int(OSRExportToWkt(osr, &wkt))
+    wkt = wkt_b
-    except CPLE_BaseError as exc:
+    if src_crs is not None:
-        raise CRSError("Could not convert to WKT. {}".format(str(exc)))
+        src_crs = CRS.from_user_input(src_crs)
-    finally:
-        _safe_osr_release(osr)
-
-    if isinstance(src_crs, str):
-        src_crs = CRS.from_string(src_crs)
-    elif isinstance(src_crs, dict):
-        src_crs = CRS(**src_crs)
-
     # The transformer at the heart of this function requires a dataset.
     # We use an in-memory VRT dataset.
     vrt_doc = _suggested_proxy_vrt_doc(
@@ -671,14 +656,12 @@ def _calculate_default_transform(
         transform=transform,
         crs=src_crs,
         gcps=gcps
-    ).decode('ascii')
+    ).decode('utf-8')
     hds = open_dataset(vrt_doc, 0x00 | 0x02 | 0x04, ['VRT'], {}, None)
-
     try:
         imgProjOptions = CSLSetNameValue(imgProjOptions, "GCPS_OK", "TRUE")
         imgProjOptions = CSLSetNameValue(imgProjOptions, "MAX_GCP_ORDER", "0")
         imgProjOptions = CSLSetNameValue(imgProjOptions, "DST_SRS", wkt)
-
         for key, val in kwargs.items():
             key = key.upper().encode('utf-8')
 
@@ -728,8 +711,6 @@ def _calculate_default_transform(
         raise CRSError(err.errmsg)
 
     finally:
-        if wkt != NULL:
-            CPLFree(wkt)
         if hTransformArg != NULL:
             GDALDestroyGenImgProjTransformer(hTransformArg)
         if hds != NULL:
@@ -963,23 +944,13 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
         self.dst_crs = CRS.from_user_input(crs) if crs is not None else self.src_crs
 
         # Convert CRSes to C WKT strings.
-        try:
+        if self.src_crs is not None:
-            if not self.src_crs:
+            src_crs_wkt_b = self.src_crs.to_wkt().encode("utf-8")
-                src_crs_wkt = NULL
+            src_crs_wkt = src_crs_wkt_b
-            else:
-                osr = _osr_from_crs(self.src_crs)
-                OSRExportToWkt(osr, &src_crs_wkt)
-        finally:
-            if osr != NULL:
-                OSRRelease(osr)
-            osr = NULL
 
         if self.dst_crs is not None:
-            try:
+            dst_crs_wkt_b = self.dst_crs.to_wkt().encode("utf-8")
-                osr = _osr_from_crs(self.dst_crs)
+            dst_crs_wkt = dst_crs_wkt_b
-                OSRExportToWkt(osr, &dst_crs_wkt)
-            finally:
-                _safe_osr_release(osr)
 
         log.debug("Exported CRS to WKT.")
 
@@ -1057,8 +1028,6 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
                     <const char **>c_warp_extras,
                 )
             finally:
-                CPLFree(dst_crs_wkt)
-                CPLFree(src_crs_wkt)
                 CSLDestroy(c_warp_extras)
                 if psWOptions != NULL:
                     GDALDestroyWarpOptions(psWOptions)
@@ -1112,8 +1081,6 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
             # If we get here it's because the tests above are buggy.
             # We raise a Python exception to indicate that.
             else:
-                CPLFree(dst_crs_wkt)
-                CPLFree(src_crs_wkt)
                 CSLDestroy(c_warp_extras)
                 if psWOptions != NULL:
                     GDALDestroyWarpOptions(psWOptions)
@@ -1152,8 +1119,6 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
                     GDALSetProjection(hds_warped, dst_crs_wkt)
 
             finally:
-                CPLFree(dst_crs_wkt)
-                CPLFree(src_crs_wkt)
                 CSLDestroy(c_warp_extras)
                 if psWOptions != NULL:
                     GDALDestroyWarpOptions(psWOptions)
@@ -1469,28 +1434,21 @@ cdef double antimeridian_max(data):
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def _transform_bounds(
-    src_crs,
+    CRS src_crs,
-    dst_crs,
+    CRS dst_crs,
     double left,
     double bottom,
     double right,
     double top,
     int densify_pts,
 ):
-    src_crs = CRS.from_user_input(src_crs)
-    dst_crs = CRS.from_user_input(dst_crs)
-
     IF (CTE_GDAL_MAJOR_VERSION, CTE_GDAL_MINOR_VERSION) >= (3, 4):
         cdef double out_left = np.inf
         cdef double out_bottom = np.inf
         cdef double out_right = np.inf
         cdef double out_top = np.inf
-        cdef OGRSpatialReferenceH src = NULL
-        cdef OGRSpatialReferenceH dst = NULL
         cdef OGRCoordinateTransformationH transform = NULL
-        src = _osr_from_crs(src_crs)
+        transform = OCTNewCoordinateTransformation(src_crs._osr, dst_crs._osr)
-        dst = _osr_from_crs(dst_crs)
-        transform = OCTNewCoordinateTransformation(src, dst)
         transform = exc_wrap_pointer(transform)
 
         # OCTTransformBounds() returns TRUE/FALSE contrary to most GDAL API functions
@@ -1506,9 +1464,6 @@ def _transform_bounds(
             exc_wrap_int(status == 0)
         finally:
             OCTDestroyCoordinateTransformation(transform)
-            _safe_osr_release(src)
-            _safe_osr_release(dst)
-
         return out_left, out_bottom, out_right, out_top
 
     ELSE:

rasterio/crs.pyx

@@ -32,8 +32,7 @@ from rasterio._err import CPLE_BaseError, CPLE_NotSupportedError
 from rasterio.errors import CRSError
 from rasterio.enums import WktVersion
 
-from rasterio._base cimport _osr_from_crs as osr_from_crs
+from rasterio._base cimport osr_set_traditional_axis_mapping_strategy
-from rasterio._base cimport _safe_osr_release, osr_set_traditional_axis_mapping_strategy
 from rasterio._err cimport exc_wrap_ogrerr, exc_wrap_int, exc_wrap_pointer
 
 
@@ -49,6 +48,14 @@ _RE_PROJ_PARAM = re.compile(r"""
 """, re.X)
 
 
+
+cdef _safe_osr_release(OGRSpatialReferenceH srs):
+    """Wrapper to handle OSR release when NULL."""
+    if srs != NULL:
+        OSRRelease(srs)
+    srs = NULL
+
+
 cdef class CRS:
     """A geographic or projected coordinate reference system.
 

rasterio/warp.py

@@ -10,6 +10,7 @@ with rasterio._loading.add_gdal_dll_directories():
     import rasterio
 
     from rasterio._base import _transform
+    from rasterio.crs import CRS
     from rasterio.enums import Resampling
     from rasterio.env import ensure_env, require_gdal_version
     from rasterio.errors import TransformError, RPCError
@@ -145,6 +146,8 @@ def transform_bounds(
     left, bottom, right, top: float
         Outermost coordinates in target coordinate reference system.
     """
+    src_crs = CRS.from_user_input(src_crs)
+    dst_crs = CRS.from_user_input(dst_crs)
     return _transform_bounds(
         src_crs,
         dst_crs,

tests/test_warp.py

@@ -285,8 +285,8 @@ def test_transform_bounds__esri_wkt():
 @pytest.mark.parametrize(
     "density,expected",
     [
-        (0, (-1684649.41338, -350356.81377, 1684649.41338, 2234551.18559)),
+        (0, (-1688721.99764, -350040.36880, 1688799.61159, 2236495.86829)),
-        (100, (-1684649.41338, -555777.79210, 1684649.41338, 2234551.18559)),
+        (100, (-1688721.99764, -555239.84875, 1688799.61159, 2236495.86829)),
     ],
 )
 def test_transform_bounds_densify(density, expected):