diff --git a/rasterio/_base.pxd b/rasterio/_base.pxd index beb7a6be..60e8288c 100644 --- a/rasterio/_base.pxd +++ b/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 diff --git a/rasterio/_base.pyx b/rasterio/_base.pyx index 3b5a2a89..c409ebb0 100644 --- a/rasterio/_base.pyx +++ b/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) - dst = _osr_from_crs(dst_crs) + cdef CRS src = CRS.from_user_input(src_crs) + cdef CRS dst = CRS.from_user_input(dst_crs) n = len(xs) x = 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, 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. - # It raises, it does not return NULL. - osr = _osr_from_crs(crs) - OSRExportToWkt(osr, &wkt) - + wkt = CRS.from_user_input(crs).to_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) diff --git a/rasterio/_env.pyx b/rasterio/_env.pyx index 2e886c21..ad8839c5 100644 --- a/rasterio/_env.pyx +++ b/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): diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx index d6c2dc4e..bee85770 100644 --- a/rasterio/_io.pyx +++ b/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 ( - _osr_from_crs, _safe_osr_release, get_driver_name, DatasetBase) +from rasterio._base cimport 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 - - osr = _osr_from_crs(crs) - OSRExportToWkt(osr, &srcwkt) + else: + crs = CRS.from_user_input(crs) + 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 diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx index 9404beae..f219bda9 100644 --- a/rasterio/_warp.pyx +++ b/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) - dst = _osr_from_crs(dst_crs) + cdef CRS src = CRS.from_user_input(src_crs) + cdef CRS dst = CRS.from_user_input(dst_crs) - try: - transform = exc_wrap_pointer(OCTNewCoordinateTransformation(src, dst)) - finally: - _safe_osr_release(src) - _safe_osr_release(dst) + transform = exc_wrap_pointer(OCTNewCoordinateTransformation(src._osr, dst._osr)) # 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: - osr = _osr_from_crs(dst_crs) - exc_wrap_int(OSRExportToWkt(osr, &wkt)) - except CPLE_BaseError as exc: - raise CRSError("Could not convert to WKT. {}".format(str(exc))) - 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) - + dst_crs = CRS.from_user_input(dst_crs) + wkt_b = dst_crs.to_wkt().encode('utf-8') + wkt = wkt_b + if src_crs is not None: + src_crs = CRS.from_user_input(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 not self.src_crs: - src_crs_wkt = NULL - else: - osr = _osr_from_crs(self.src_crs) - OSRExportToWkt(osr, &src_crs_wkt) - finally: - if osr != NULL: - OSRRelease(osr) - osr = NULL + if self.src_crs is not None: + src_crs_wkt_b = self.src_crs.to_wkt().encode("utf-8") + src_crs_wkt = src_crs_wkt_b if self.dst_crs is not None: - try: - osr = _osr_from_crs(self.dst_crs) - OSRExportToWkt(osr, &dst_crs_wkt) - finally: - _safe_osr_release(osr) + dst_crs_wkt_b = self.dst_crs.to_wkt().encode("utf-8") + dst_crs_wkt = dst_crs_wkt_b log.debug("Exported CRS to WKT.") @@ -1057,8 +1028,6 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase): 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, - dst_crs, + CRS src_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) - dst = _osr_from_crs(dst_crs) - transform = OCTNewCoordinateTransformation(src, dst) + transform = OCTNewCoordinateTransformation(src_crs._osr, dst_crs._osr) 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: diff --git a/rasterio/crs.pyx b/rasterio/crs.pyx index ff10d18f..7cb5bc74 100644 --- a/rasterio/crs.pyx +++ b/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 _safe_osr_release, osr_set_traditional_axis_mapping_strategy +from rasterio._base cimport 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. diff --git a/rasterio/warp.py b/rasterio/warp.py index 59894733..534653c2 100644 --- a/rasterio/warp.py +++ b/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, diff --git a/tests/test_warp.py b/tests/test_warp.py index 0047e601..9e6b9e6b 100644 --- a/tests/test_warp.py +++ b/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)), - (100, (-1684649.41338, -555777.79210, 1684649.41338, 2234551.18559)), + (0, (-1688721.99764, -350040.36880, 1688799.61159, 2236495.86829)), + (100, (-1688721.99764, -555239.84875, 1688799.61159, 2236495.86829)), ], ) def test_transform_bounds_densify(density, expected):