.travis.yml

@@ -14,9 +14,9 @@ env:
     - GDALVERSION="1.11.5"
     - GDALVERSION="2.0.3"
     - GDALVERSION="2.1.4"
-    - GDALVERSION="2.2.1"
-    - GDALVERSION="2.2.2"
-    - GDALVERSION="2.2.3"
+    - GDALVERSION="2.2.4"
+    - GDALVERSION="2.3.3"
+    - GDALVERSION="2.4.0"
     - GDALVERSION="master"
 
 matrix:
@@ -43,7 +43,7 @@ before_script: # configure a headless display to test matplotlib
 before_install:
   - python -m pip install -U pip
   - python -m pip install wheel
-  - travis_wait . ./scripts/travis_gdal_install.sh
+  - travis_wait 20 bash ./scripts/travis_gdal_install.sh
   - export PATH=$GDALINST/gdal-$GDALVERSION/bin:$PATH
   - export LD_LIBRARY_PATH=$GDALINST/gdal-$GDALVERSION/lib:$LD_LIBRARY_PATH
 install:

AUTHORS.txt

@@ -68,3 +68,4 @@ Authors
 * Vincent Schut
 * Alan D.
 * grovduck
+* Dan Little

CHANGES.txt

 Changes
 =======
 
+1.0.14 (2019-01-22)
+-------------------
+
+- The _CRS class has been refactored so that a WKT representation, rather than
+  PROJ4 representation, is the canonical form. This resolves issues #1397 and
+  #1587 specifically, and an entire category of issues discussed in
+  https://rasterio.groups.io/g/dev/message/68.
+- Native support for Google Cloud Storage using "gs://" URLs has been added
+  (#1577).
+- On entering a dataset context (DatasetBase.__enter__) a new anonymous GDAL
+  environment is created if needed and is entered. This makes `with
+  rasterio.open(...) as dataset:` roughly equivalent to `with
+  rasterio.open(...) as dataset, Env():`. This helps prevent bugs when datasets
+  are created and then used later or are used in different scopes.
+
 1.0.13 (2018-12-14)
 -------------------
 

debian/changelog

-rasterio (1.0.13-2) UNRELEASED; urgency=medium
+rasterio (1.0.14-1) unstable; urgency=medium
 
+  * Team upload.
+  * New upstream release.
   * Bump Standards-Version to 4.3.0, no changes.
 
- -- Bas Couwenberg <sebastic@debian.org>  Tue, 25 Dec 2018 23:11:57 +0100
+ -- Bas Couwenberg <sebastic@debian.org>  Wed, 23 Jan 2019 06:59:47 +0100
 
 rasterio (1.0.13-1) unstable; urgency=medium
 

docs/quickstart.rst

@@ -202,7 +202,7 @@ The coordinates of the center of the image can be computed like this.
 
 .. code-block:: pycon
 
-    >>> dataset.xy(dataset.width // 2, dataset.height // 2)
+    >>> dataset.xy(dataset.height // 2, dataset.width // 2)
     (476550.0, 4149150.0)
 
 Creating data

environment.yml

 name: _rasterio
 channels:
 - conda-forge
-- defaults
 dependencies:
-- python>=3.5
-- cython
-- libgdal
-- numpy
+- python=3.6.6
+- libgdal=2.3.2

rasterio/__init__.py

@@ -42,7 +42,7 @@ import rasterio.path
 
 
 __all__ = ['band', 'open', 'pad', 'Env']
-__version__ = "1.0.13"
+__version__ = "1.0.14"
 __gdal_version__ = gdal_version()
 
 # Rasterio attaches NullHandler to the 'rasterio' logger and its

rasterio/_base.pxd

@@ -27,6 +27,7 @@ cdef class DatasetBase:
     cdef public object _offsets
     cdef public object _read
     cdef public object _gcps
+    cdef public object _env
     cdef GDALDatasetH handle(self) except NULL
     cdef GDALRasterBandH band(self, int bidx) except NULL
 

rasterio/_base.pyx

@@ -23,7 +23,7 @@ from rasterio.coords import BoundingBox
 from rasterio.crs import CRS
 from rasterio.enums import (
     ColorInterp, Compression, Interleaving, MaskFlags, PhotometricInterp)
-from rasterio.env import Env
+from rasterio.env import Env, env_ctx_if_needed
 from rasterio.errors import (
     RasterioIOError, CRSError, DriverRegistrationError, NotGeoreferencedWarning,
     RasterBlockError, BandOverviewError)
@@ -262,56 +262,19 @@ cdef class DatasetBase(object):
 
     def _handle_crswkt(self, wkt):
         """Return the GDAL dataset's stored CRS"""
-        cdef OGRSpatialReferenceH osr = NULL
-        cdef const char *auth_key = NULL
-        cdef const char *auth_val = NULL
-
-        if not wkt:
-            log.debug("No projection detected.")
-            return None
-
-        wkt_b = wkt.encode('utf-8')
-        cdef const char *wkt_c = wkt_b
-
-        try:
-
-            osr = exc_wrap_pointer(OSRNewSpatialReference(wkt_c))
-            log.debug("Got coordinate system")
-
-            retval = OSRAutoIdentifyEPSG(osr)
-
-            if retval > 0:
-                log.debug("Failed to auto identify EPSG: %d", retval)
-
-            else:
-                log.debug("Auto identified EPSG: %d", retval)
-
-                try:
-                    auth_key = OSRGetAuthorityName(osr, NULL)
-                    auth_val = OSRGetAuthorityCode(osr, NULL)
-
-                except CPLE_NotSupportedError as exc:
-                    log.debug("{}".format(exc))
-
-                if auth_key != NULL and auth_val != NULL:
-                    return CRS({'init': u'{}:{}'.format(auth_key.lower(), auth_val)})
-
+        # No dialect morphing, if the dataset was created using software
+        # "speaking" the Esri dialect, we will read Esri WKT.
+        if wkt:
             return CRS.from_wkt(wkt)
-
-        except CPLE_BaseError as exc:
-            raise CRSError("{}".format(exc))
-
-        finally:
-             _safe_osr_release(osr)
+        else:
+            return CRS()
 
     def read_crs(self):
         """Return the GDAL dataset's stored CRS"""
-        cdef const char *wkt_b = NULL
-
-        wkt_b = GDALGetProjectionRef(self._hds)
-        if wkt_b == NULL:
-            raise ValueError("Unexpected NULL spatial reference")
+        cdef const char *wkt_b = GDALGetProjectionRef(self.handle())
         wkt = wkt_b
+        if wkt == NULL:
+            raise ValueError("Unexpected NULL spatial reference")
         return self._handle_crswkt(wkt)
 
     def read_transform(self):
@@ -333,20 +296,19 @@ cdef class DatasetBase(object):
         if self._hds != NULL:
             GDALClose(self._hds)
         self._hds = NULL
-        log.debug("Dataset %r has been stopped.", self)
 
     def close(self):
         self.stop()
         self._closed = True
-        log.debug("Dataset %r has been closed.", self)
 
     def __enter__(self):
-        log.debug("Entering Dataset %r context.", self)
+        self._env = env_ctx_if_needed()
+        self._env.__enter__()
         return self
 
     def __exit__(self, type, value, traceback):
+        self._env.__exit__()
         self.close()
-        log.debug("Exited Dataset %r context.", self)
 
     def __dealloc__(self):
         if self._hds != NULL:
@@ -1329,6 +1291,7 @@ cdef OGRSpatialReferenceH _osr_from_crs(object crs) except NULL:
         if retval:
             _safe_osr_release(osr)
             raise CRSError("Invalid CRS: {!r}".format(crs))
+        exc_wrap_int(OSRMorphFromESRI(osr))
     except CPLE_BaseError as exc:
         _safe_osr_release(osr)
         raise CRSError(str(exc))

rasterio/_crs.pxd

+# _CRS class definition
+
+include "gdal.pxi"
+
+
+cdef class _CRS:
+
+    cdef OGRSpatialReferenceH _osr

rasterio/_crs.pyx

 """Coordinate reference systems, class and functions.
 """
 
-include "gdal.pxi"
-
-import json
 import logging
 
 from rasterio._err import CPLE_BaseError, CPLE_NotSupportedError
-from rasterio.compat import UserDict, string_types
+from rasterio.compat import string_types
 from rasterio.errors import CRSError
 from rasterio.env import env_ctx_if_needed
 
 from rasterio._base cimport _osr_from_crs as osr_from_crs
 from rasterio._base cimport _safe_osr_release
-from rasterio._err cimport exc_wrap_int
+from rasterio._err cimport exc_wrap_ogrerr, exc_wrap_int, exc_wrap_pointer
 
 
 log = logging.getLogger(__name__)
 
 
-class _CRS(UserDict):
-    """CRS base class."""
+cdef class _CRS(object):
+    """Cython extension class"""
+
+    def __cinit__(self):
+        self._osr = OSRNewSpatialReference(NULL)
+
+    def __dealloc__(self):
+        _safe_osr_release(self._osr)
 
     @property
     def is_geographic(self):
@@ -29,16 +32,12 @@ class _CRS(UserDict):
         Returns
         -------
         bool
-        """
-        cdef OGRSpatialReferenceH osr_crs = NULL
-        cdef int retval
 
+        """
         try:
-            osr_crs = osr_from_crs(self)
-            retval = OSRIsGeographic(osr_crs)
-            return bool(retval == 1)
-        finally:
-            _safe_osr_release(osr_crs)
+            return bool(OSRIsGeographic(self._osr) == 1)
+        except CPLE_BaseError as exc:
+            raise CRSError("{}".format(exc))
 
     @property
     def is_projected(self):
@@ -47,59 +46,58 @@ class _CRS(UserDict):
         Returns
         -------
         bool
-        """
-        cdef OGRSpatialReferenceH osr_crs = NULL
-        cdef int retval
 
+        """
         try:
-            osr_crs = osr_from_crs(self)
-            retval = OSRIsProjected(osr_crs)
-            return bool(retval == 1)
-        finally:
-            _safe_osr_release(osr_crs)
+            return bool(OSRIsProjected(self._osr) == 1)
+        except CPLE_BaseError as exc:
+            raise CRSError("{}".format(exc))
 
     def __eq__(self, other):
-        cdef OGRSpatialReferenceH osr_crs1 = NULL
-        cdef OGRSpatialReferenceH osr_crs2 = NULL
-        cdef int retval
+        cdef OGRSpatialReferenceH osr_s = NULL
+        cdef OGRSpatialReferenceH osr_o = NULL
+        cdef _CRS crs_o = other
 
         try:
-            if (
-                isinstance(other, self.__class__) and
-                self.data == other.data
-            ):
-                return True
-
-            if not self or not other:
-                return not self and not other
-
-            osr_crs1 = osr_from_crs(self)
-            osr_crs2 = osr_from_crs(other)
-            retval = OSRIsSame(osr_crs1, osr_crs2)
-            return bool(retval == 1)
+            osr_s = exc_wrap_pointer(OSRClone(self._osr))
+            exc_wrap_ogrerr(OSRMorphFromESRI(osr_s))
+            osr_o = exc_wrap_pointer(OSRClone(crs_o._osr))
+            exc_wrap_ogrerr(OSRMorphFromESRI(osr_o))
+            return bool(OSRIsSame(osr_s, osr_o) == 1)
 
         finally:
-            _safe_osr_release(osr_crs1)
-            _safe_osr_release(osr_crs2)
+            _safe_osr_release(osr_s)
+            _safe_osr_release(osr_o)
 
-    @property
-    def wkt(self):
+    def to_wkt(self, morph_to_esri_dialect=False):
         """An OGC WKT representation of the CRS
 
+        Parameters
+        ----------
+        morph_to_esri_dialect : bool, optional
+            Whether or not to morph to the Esri dialect of WKT
+
         Returns
         -------
         str
+
         """
-        cdef char *srcwkt = NULL
-        cdef OGRSpatialReferenceH osr = NULL
+        cdef char *conv_wkt = NULL
 
         try:
-            osr = osr_from_crs(self)
-            OSRExportToWkt(osr, &srcwkt)
-            return srcwkt.decode('utf-8')
+            if morph_to_esri_dialect:
+                exc_wrap_ogrerr(OSRMorphToESRI(self._osr))
+
+            exc_wrap_ogrerr(OSRExportToWkt(self._osr, &conv_wkt))
+
+        except CPLE_BaseError as exc:
+            raise CRSError("Cannot convert to WKT. {}".format(exc))
+
+        else:
+            return conv_wkt.decode('utf-8')
+
         finally:
-            CPLFree(srcwkt)
-            _safe_osr_release(osr)
+            CPLFree(conv_wkt)
 
     def to_epsg(self):
         """The epsg code of the CRS
@@ -107,11 +105,13 @@ class _CRS(UserDict):
         Returns
         -------
         int
+
         """
         cdef OGRSpatialReferenceH osr = NULL
 
         try:
-            osr = osr_from_crs(self)
+            osr = exc_wrap_pointer(OSRClone(self._osr))
+            exc_wrap_ogrerr(OSRMorphFromESRI(osr))
             if OSRAutoIdentifyEPSG(osr) == 0:
                 epsg_code = OSRGetAuthorityCode(osr, NULL)
                 return int(epsg_code.decode('utf-8'))
@@ -136,149 +136,173 @@ class _CRS(UserDict):
         Returns
         -------
         CRS
+
         """
         if int(code) <= 0:
             raise CRSError("EPSG codes are positive integers")
-        return cls(init="epsg:%s" % code, no_defs=True)
+        return cls.from_proj4('+init=epsg:{}'.format(code))
 
-    @classmethod
-    def from_string(cls, s):
-        """Make a CRS from an EPSG, PROJ, or WKT string
+    @staticmethod
+    def from_proj4(proj):
+        """Make a CRS from a PROJ4 string
 
         Parameters
         ----------
-        s : str
-            An EPSG, PROJ, or WKT string.
+        proj : str
+            A PROJ4 string like "+proj=longlat ..."
 
         Returns
         -------
         CRS
+
         """
-        if not s:
-            raise CRSError("CRS is empty or invalid: {!r}".format(s))
+        cdef _CRS obj = _CRS.__new__(_CRS)
 
-        elif s.strip().upper().startswith('EPSG:'):
-            auth, val = s.strip().split(':')
-            if not val:
-                raise CRSError("Invalid CRS: {!r}".format(s))
-            return cls.from_epsg(val)
+        # Filter out nonsensical items.
+        items_filtered = []
+        items = proj.split()
+        for item in items:
+            parts = item.split('=')
+            if len(parts) == 2 and parts[1] in ('false', 'False'):
+                continue
+            items_filtered.append(item)
+
+        proj = ' '.join(items_filtered)
+        proj_b = proj.encode('utf-8')
 
-        elif '{' in s:
-            # may be json, try to decode it
         try:
-                val = json.loads(s, strict=False)
-            except ValueError:
-                raise CRSError('CRS appears to be JSON but is not valid')
+            exc_wrap_ogrerr(exc_wrap_int(OSRImportFromProj4(obj._osr, <const char *>proj_b)))
+
+        except CPLE_BaseError as exc:
+            raise CRSError("The PROJ4 dict could not be understood. {}".format(exc))
 
-            if not val:
-                raise CRSError("CRS is empty JSON")
         else:
-                return cls(**val)
+            return obj
 
-        elif '+' in s and '=' in s:
+    @staticmethod
+    def from_dict(initialdata=None, **kwargs):
+        """Make a CRS from a PROJ dict
 
-            parts = [o.lstrip('+') for o in s.strip().split()]
+        Parameters
+        ----------
+        initialdata : mapping, optional
+            A dictionary or other mapping
+        kwargs : mapping, optional
+            Another mapping. Will be overlaid on the initialdata.
 
-            def parse(v):
-                if v in ('True', 'true'):
-                    return True
-                elif v in ('False', 'false'):
-                    return False
-                else:
-                    try:
-                        return int(v)
-                    except ValueError:
-                        pass
-                    try:
-                        return float(v)
-                    except ValueError:
-                        return v
+        Returns
+        -------
+        CRS
 
-            items = map(
-                lambda kv: len(kv) == 2 and (kv[0], parse(kv[1])) or (kv[0], True),
-                (p.split('=') for p in parts))
+        """
+        data = dict(initialdata or {})
+        data.update(**kwargs)
+        data = {k: v for k, v in data.items() if k in all_proj_keys}
 
-            out = cls((k, v) for k, v in items if k in all_proj_keys)
+        # always use lowercase 'epsg'.
+        if 'init' in data:
+            data['init'] = data['init'].replace('EPSG:', 'epsg:')
 
-            if not out:
-                raise CRSError("CRS is empty or invalid: {}".format(s))
+        proj = ' '.join(['+{}={}'.format(key, val) for key, val in data.items()])
+        b_proj = proj.encode('utf-8')
 
-            return out
+        cdef _CRS obj = _CRS.__new__(_CRS)
 
+        try:
+            exc_wrap_ogrerr(OSRImportFromProj4(obj._osr, <const char *>b_proj))
+        except CPLE_BaseError as exc:
+            raise CRSError("The PROJ4 dict could not be understood. {}".format(exc))
         else:
-            return cls.from_wkt(s)
+            return obj
 
-    @classmethod
-    def from_wkt(cls, s):
+    @staticmethod
+    def from_wkt(wkt, morph_from_esri_dialect=False):
         """Make a CRS from a WKT string
 
         Parameters
         ----------
-        s : str
+        wkt : str
             A WKT string.
+        morph_from_esri_dialect : bool, optional
+            If True, items in the input using Esri's dialect of WKT
+            will be replaced by OGC standard equivalents.
 
         Returns
         -------
         CRS
+
         """
-        cdef char *prj = NULL
-        cdef OGRSpatialReferenceH osr = OSRNewSpatialReference(NULL)
+        cdef char *wkt_c = NULL
 
-        if isinstance(s, string_types):
-            b_s = s.encode('utf-8')
-            log.debug("Encoded WKT: %r", b_s)
+        if not isinstance(wkt, string_types):
+            raise ValueError("A string is expected")
 
-        try:
-            exc_wrap_int(OSRSetFromUserInput(osr, <const char *>b_s))
-        except CPLE_NotSupportedError as exc:
-            log.debug("{}".format(exc))
-        except CPLE_BaseError as exc:
-            _safe_osr_release(osr)
-            raise CRSError(str(exc))
+        wkt_b= wkt.encode('utf-8')
+        wkt_c = wkt_b
 
-        try:
-            exc_wrap_int(OSRExportToProj4(osr, &prj))
-        except CPLE_NotSupportedError as exc:
-            log.debug("{}".format(exc))
+        cdef _CRS obj = _CRS.__new__(_CRS)
 
         try:
-            return cls.from_string(prj.decode('utf-8'))
-        except CRSError:
-            if OSRMorphFromESRI(osr) == 0:
-                OSRExportToProj4(osr, &prj)
-                return cls.from_string(prj.decode('utf-8'))
-            else:
-                raise
-        finally:
-            CPLFree(prj)
-            _safe_osr_release(osr)
+            errcode = exc_wrap_ogrerr(OSRImportFromWkt(obj._osr, &wkt_c))
 
-    @classmethod
-    def from_user_input(cls, value):
-        """Make a CRS from various input
+            if morph_from_esri_dialect:
+                exc_wrap_ogrerr(OSRMorphFromESRI(obj._osr))
 
-        Dispatches to from_epsg, from_proj, or from_string
+        except CPLE_BaseError as exc:
+            raise CRSError("The WKT could not be parsed. {}".format(exc))
 
-        Parameters
-        ----------
-        value : obj
-            A Python int, dict, or str.
+        else:
+            return obj
+
+    def to_dict(self):
+        """Convert CRS to a PROJ4 dict
 
         Returns
         -------
-        CRS
+        dict
+
         """
-        if isinstance(value, _CRS):
-            return value
-        elif isinstance(value, int):
-            return cls.from_epsg(value)
-        elif isinstance(value, dict):
-            return cls(**value)
-        elif isinstance(value, string_types):
-            return cls.from_string(value)
+        cdef OGRSpatialReferenceH osr = NULL
+        cdef char *proj_c = NULL
+
+        try:
+            osr = exc_wrap_pointer(OSRClone(self._osr))
+            exc_wrap_ogrerr(OSRMorphFromESRI(osr))
+            exc_wrap_ogrerr(OSRExportToProj4(osr, &proj_c))
+
+        except CPLE_BaseError as exc:
+            raise CRSError("The WKT could not be parsed. {}".format(exc))
+
+        else:
+            proj_b = proj_c
+            proj = proj_b.decode('utf-8')
+
+        finally:
+            CPLFree(proj_c)
+            _safe_osr_release(osr)
+
+        parts = [o.lstrip('+') for o in proj.strip().split()]
+
+        def parse(v):
+            if v in ('True', 'true'):
+                return True
+            elif v in ('False', 'false'):
+                return False
             else:
-            raise CRSError("CRS is invalid: {!r}".format(value))
+                try:
+                    return int(v)
+                except ValueError:
+                    pass
+                try:
+                    return float(v)
+                except ValueError:
+                    return v
+
+        items = map(
+            lambda kv: len(kv) == 2 and (kv[0], parse(kv[1])) or (kv[0], True),
+            (p.split('=') for p in parts))
 
+        return {k: v for k, v in items if k in all_proj_keys and v is not False}
 
 # Below is the big list of PROJ4 parameters from
 # http://trac.osgeo.org/proj/wiki/GenParms.

rasterio/_err.pxd

@@ -4,7 +4,12 @@ cdef extern from "cpl_vsi.h":
 
     ctypedef FILE VSILFILE
 
+cdef extern from "ogr_core.h":
+
+    ctypedef int OGRErr
+
 
 cdef int exc_wrap_int(int retval) except -1
+cdef OGRErr exc_wrap_ogrerr(OGRErr retval) except -1
 cdef void *exc_wrap_pointer(void *ptr) except NULL
 cdef VSILFILE *exc_wrap_vsilfile(VSILFILE *f) except NULL

rasterio/_err.pyx

@@ -183,6 +183,17 @@ cdef int exc_wrap_int(int err) except -1:
     return err
 
 
+cdef OGRErr exc_wrap_ogrerr(OGRErr err) except -1:
+    """Wrap a function that returns OGRErr but does not use the
+    CPL error stack.
+
+    """
+    if err == 0:
+        return err
+    else:
+        raise CPLE_BaseError(3, err, "OGR Error code {}".format(err))
+
+
 cdef void *exc_wrap_pointer(void *ptr) except NULL:
     """Wrap a GDAL/OGR function that returns GDALDatasetH etc (void *)
 

rasterio/_io.pyx

@@ -1209,59 +1209,11 @@ cdef class DatasetWriterBase(DatasetReaderBase):
 
     def _set_crs(self, crs):
         """Writes a coordinate reference system to the dataset."""
-        cdef char *proj_c = NULL
-        cdef char *wkt = NULL
-        cdef OGRSpatialReferenceH osr = NULL
-
-        osr = OSRNewSpatialReference(NULL)
-        if osr == NULL:
-            raise ValueError("Null spatial reference")
-        params = []
-
-        log.debug("Input CRS: %r", crs)
-
-        if isinstance(crs, dict):
-            crs = CRS(crs)
-        if isinstance(crs, CRS):
-            # EPSG is a special case.
-            init = crs.get('init')
-            if init:
-                auth, val = init.split(':')
-                if auth.upper() == 'EPSG':
-                    OSRImportFromEPSG(osr, int(val))
-            else:
-                for k, v in crs.items():
-                    if v is True or (k in ('no_defs', 'wktext') and v):
-                        params.append("+%s" % k)
-                    else:
-                        params.append("+%s=%s" % (k, v))
-                proj = " ".join(params)
-                log.debug("PROJ.4 to be imported: %r", proj)
-                proj_b = proj.encode('utf-8')
-                proj_c = proj_b
-                OSRImportFromProj4(osr, proj_c)
-        # Fall back for CRS strings like "EPSG:3857."
-        elif isinstance(crs, str) or crs is None:
-            if not crs:
-                crs = ''
-            proj_b = crs.encode('utf-8')
-            proj_c = proj_b
-            OSRSetFromUserInput(osr, proj_c)
-        else:
-            raise CRSError(
-                "{!r} does not define a valid CRS".format(crs))
-
-        # Fixup, export to WKT, and set the GDAL dataset's projection.
-        OSRFixup(osr)
-        OSRExportToWkt(osr, <char**>&wkt)
-        wkt_b = wkt
-        log.debug("Exported WKT: %s", wkt_b.decode('utf-8'))
-        GDALSetProjection(self._hds, wkt)
-
-        CPLFree(wkt)
-        _safe_osr_release(osr)
+        crs = CRS.from_user_input(crs)
+        wkt_b = crs.to_wkt().encode('utf-8')
+        cdef const char *wkt_c = wkt_b
+        exc_wrap_int(GDALSetProjection(self.handle(), wkt_c))
         self._crs = crs
-        log.debug("Self CRS: %r", self._crs)
 
     def _set_all_descriptions(self, value):
         """Supports the descriptions property setter"""

rasterio/_warp.pyx

@@ -14,7 +14,7 @@ import numpy as np
 import rasterio
 from rasterio._base import gdal_version
 from rasterio._err import (
-    CPLE_IllegalArgError, CPLE_NotSupportedError,
+    CPLE_BaseError, CPLE_IllegalArgError, CPLE_NotSupportedError,
     CPLE_AppDefinedError, CPLE_OpenFailedError)
 from rasterio import dtypes
 from rasterio.control import GroundControlPoint
@@ -539,8 +539,12 @@ def _calculate_default_transform(src_crs, dst_crs, width, height,
     else:
         transform = None
 
+    try:
         osr = _osr_from_crs(dst_crs)
-    OSRExportToWkt(osr, &wkt)
+        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):

rasterio/crs.py

@@ -3,44 +3,208 @@
 Notes
 -----
 
-In Rasterio 1.0, coordinate reference system support is limited to the
-CRS that can be described by PROJ parameters.
+In Rasterio versions <= 1.0.13, coordinate reference system support was limited
+to the CRS that can be described by PROJ parameters. This limitation is gone in
+versions >= 1.0.14. Any CRS that can be defined using WKT (version 1) may be
+used.
 
 """
 
+import collections
+import json
+
 from rasterio._crs import _CRS, all_proj_keys
 from rasterio.compat import string_types
+from rasterio.errors import CRSError
 
 
-class CRS(_CRS):
-    """A container class for coordinate reference system info
+class CRS(collections.Mapping):
+    """A geographic or projected coordinate reference system
 
     CRS objects may be created by passing PROJ parameters as keyword
-    arguments to the standard constructor or by passing EPSG codes,
-    PROJ strings, or WKT strings to the from_epsg and from_string
-    class methods.
+    arguments to the standard constructor or by passing EPSG codes, PROJ
+    mappings, PROJ strings, or WKT strings to the from_epsg, from_dict,
+    from_string, or from_wkt class methods or static methods.
 
     Examples
     --------
 
-    The constructor takes PROJ parameters as keyword arguments.
+    The from_dict method takes PROJ parameters as keyword arguments.
 
-    >>> crs = CRS(init='epsg:3005')
+    >>> crs = CRS.from_dict(init='epsg:3005')
 
-    EPSG codes may be used with the from_epsg class method.
+    EPSG codes may be used with the from_epsg method.
 
     >>> crs = CRS.from_epsg(3005)
 
+    The from_string method takes a variety of input.
+
+    >>> crs = CRS.from_string('EPSG:3005')
+
+    """
+    def __init__(self, initialdata=None, **kwargs):
+        """Make a CRS from a PROJ dict or mapping
+
+        Parameters
+        ----------
+        initialdata : mapping, optional
+            A dictionary or other mapping
+        kwargs : mapping, optional
+            Another mapping. Will be overlaid on the initialdata.
+
+        Returns
+        -------
+        CRS
+
+        """
+        self._wkt = None
+        self._data = None
+        self._crs = None
+
+        if initialdata or kwargs:
+
+            data = dict(initialdata or {})
+            data.update(**kwargs)
+            data = {k: v for k, v in data.items() if k in all_proj_keys}
+
+            # always use lowercase 'epsg'.
+            if 'init' in data:
+                data['init'] = data['init'].replace('EPSG:', 'epsg:')
+
+            proj = ' '.join(['+{}={}'.format(key, val) for key, val in data.items()])
+            self._crs = _CRS.from_proj4(proj)
+
+        else:
+            self._crs = _CRS()
+
+    def __getitem__(self, item):
+        return self.data[item]
+
+    def __iter__(self):
+        return iter(self.data)
+
+    def __len__(self):
+        return len(self.data)
+
+    def __bool__(self):
+        return bool(self.wkt)
+
+    __nonzero__ = __bool__
+
+    def __eq__(self, other):
+        other = CRS.from_user_input(other)
+        return (self._crs == other._crs)
+
+    def to_proj4(self):
+        """Convert CRS to a PROJ4 string
+
+        Returns
+        -------
+        str
+
         """
+        return ' '.join(['+{}={}'.format(key, val) for key, val in self.data.items()])
+
+    def to_wkt(self, morph_to_esri_dialect=False):
+        """Convert CRS to its OGC WKT representation
+
+        Parameters
+        ----------
+        morph_to_esri_dialect : bool, optional
+            Whether or not to morph to the Esri dialect of WKT
+
+        Returns
+        -------
+        str
+
+        """
+        return self._crs.to_wkt(morph_to_esri_dialect=morph_to_esri_dialect)
+
+    @property
+    def wkt(self):
+        """An OGC WKT representation of the CRS
+
+        Returns
+        -------
+        str
+
+        """
+        if not self._wkt:
+            self._wkt = self.to_wkt()
+        return self._wkt
+
+    def to_epsg(self):
+        """The epsg code of the CRS
+
+        Returns None if there is no corresponding EPSG code.
+
+        Returns
+        -------
+        int
+
+        """
+        return self._crs.to_epsg()
+
+    def to_dict(self):
+        """Convert CRS to a PROJ4 dict
+
+        Notes
+        -----
+        If there is a corresponding EPSG code, it will be used.
+
+        Returns
+        -------
+        dict
+
+        """
+        epsg_code = self.to_epsg()
+        if epsg_code:
+            return {'init': 'epsg:{}'.format(epsg_code)}
+        else:
+            return self._crs.to_dict()
+
+    @property
+    def data(self):
+        """A PROJ4 dict representation of the CRS"""
+        if not self._data:
+            self._data = self.to_dict()
+        return self._data
+
+    @property
+    def is_geographic(self):
+        """Test that the CRS is a geographic CRS
+
+        Returns
+        -------
+        bool
+
+        """
+        return self._crs.is_geographic
+
+    @property
+    def is_projected(self):
+        """Test that the CRS is a projected CRS
+
+        Returns
+        -------
+        bool
+
+        """
+        return self._crs.is_projected
 
     @property
     def is_valid(self):
-        """Test if the CRS is a valid geographic or projected
-        coordinate reference system
+        """Test that the CRS is a geographic or projected CRS
+
+        Notes
+        -----
+        There are other types of CRS, such as compound or local or
+        engineering CRS, but these are not supported in Rasterio 1.0.
 
         Returns
         -------
         bool
+
         """
         return self.is_geographic or self.is_projected
 
@@ -51,48 +215,195 @@ class CRS(_CRS):
         Returns
         -------
         bool
+
         """
-        for val in self.values():
-            if isinstance(val, string_types) and val.lower().startswith('epsg'):
-                return True
+        try:
+            return bool(self.to_epsg())
+        except CRSError:
             return False
 
     def to_string(self):
-        """Turn CRS into a PROJ string
+        """Convert CRS to a PROJ4 or WKT string
 
         Notes
         -----
 
-        Mapping keys are tested against the ``all_proj_keys`` list. Values of
-        ``True`` are omitted, leaving the key bare: {'no_defs': True} -> "+no_defs"
-        and items where the value is otherwise not a str, int, or float are
-        omitted.
+        Mapping keys are tested against the ``all_proj_keys`` list.
+        Values of ``True`` are omitted, leaving the key bare:
+        {'no_defs': True} -> "+no_defs" and items where the value is
+        otherwise not a str, int, or float are omitted.
 
         Returns
         -------
         str
+
         """
-        items = []
-        for k, v in sorted(filter(
-                lambda x: x[0] in all_proj_keys and x[1] is not False and (
-                    isinstance(x[1], (bool, int, float)) or
-                    isinstance(x[1], string_types)),
-                self.items())):
-            items.append("+" + "=".join(map(str, filter(
-                lambda y: (y or y == 0) and y is not True, (k, v)))))
-        return " ".join(items)
+        epsg_code = self.to_epsg()
+        if epsg_code:
+            return 'EPSG:{}'.format(epsg_code)
+        else:
+            return self.to_wkt() or self.to_proj4()
+
+    __str__ = to_string
 
     def __repr__(self):
-        return "CRS({})".format(dict.__repr__(self.data))
+        epsg_code = self.to_epsg()
+        if epsg_code:
+            return "CRS.from_dict(init='epsg:{}')".format(epsg_code)
+        else:
+            return "CRS.from_wkt('{}')".format(self.wkt)
 
-    def __str__(self):
-        return self.to_string()
+    @classmethod
+    def from_epsg(cls, code):
+        """Make a CRS from an EPSG code
 
-    def to_dict(self):
-        """Turn CRS object into a dict
+        Parameters
+        ----------
+        code : int or str
+            An EPSG code. Strings will be converted to integers.
+
+        Notes
+        -----
+        The input code is not validated against an EPSG database.
 
         Returns
         -------
-        dict
+        CRS
+
+        """
+        obj = cls()
+        obj._crs = _CRS.from_epsg(code)
+        return obj
+
+    @classmethod
+    def from_string(cls, string, morph_from_esri_dialect=False):
+        """Make a CRS from an EPSG, PROJ, or WKT string
+
+        Parameters
+        ----------
+        string : str
+            An EPSG, PROJ, or WKT string.
+        morph_from_esri_dialect : bool, optional
+            If True, items in the input using Esri's dialect of WKT
+            will be replaced by OGC standard equivalents.
+
+        Returns
+        -------
+        CRS
+
+        """
+        if not string:
+            raise CRSError("CRS is empty or invalid: {!r}".format(string))
+
+        elif string.strip().upper().startswith('EPSG:'):
+            auth, val = string.strip().split(':')
+            if not val:
+                raise CRSError("Invalid CRS: {!r}".format(string))
+            return cls.from_epsg(val)
+
+        elif string.startswith('{') or string.startswith('['):
+            # may be json, try to decode it
+            try:
+                val = json.loads(string, strict=False)
+            except ValueError:
+                raise CRSError('CRS appears to be JSON but is not valid')
+
+            if not val:
+                raise CRSError("CRS is empty JSON")
+            else:
+                return cls.from_dict(**val)
+
+        elif '+' in string and '=' in string:
+            return cls.from_proj4(string)
+
+        else:
+            return cls.from_wkt(string, morph_from_esri_dialect=morph_from_esri_dialect)
+
+    @classmethod
+    def from_proj4(cls, proj):
+        """Make a CRS from a PROJ4 string
+
+        Parameters
+        ----------
+        proj : str
+            A PROJ4 string like "+proj=longlat ..."
+
+        Returns
+        -------
+        CRS
+
+        """
+        obj = cls()
+        obj._crs = _CRS.from_proj4(proj)
+        return obj
+
+    @classmethod
+    def from_dict(cls, initialdata=None, **kwargs):
+        """Make a CRS from a PROJ dict
+
+        Parameters
+        ----------
+        initialdata : mapping, optional
+            A dictionary or other mapping
+        kwargs : mapping, optional
+            Another mapping. Will be overlaid on the initialdata.
+
+        Returns
+        -------
+        CRS
+
+        """
+        obj = cls()
+        obj._crs = _CRS.from_dict(initialdata, **kwargs)
+        return obj
+
+    @classmethod
+    def from_wkt(cls, wkt, morph_from_esri_dialect=False):
+        """Make a CRS from a WKT string
+
+        Parameters
+        ----------
+        wkt : str
+            A WKT string.
+        morph_from_esri_dialect : bool, optional
+            If True, items in the input using Esri's dialect of WKT
+            will be replaced by OGC standard equivalents.
+
+        Returns
+        -------
+        CRS
+
+        """
+        obj = cls()
+        obj._crs = _CRS.from_wkt(wkt, morph_from_esri_dialect=morph_from_esri_dialect)
+        return obj
+
+    @classmethod
+    def from_user_input(cls, value, morph_from_esri_dialect=False):
+        """Make a CRS from various input
+
+        Dispatches to from_epsg, from_proj, or from_string
+
+        Parameters
+        ----------
+        value : obj
+            A Python int, dict, or str.
+        morph_from_esri_dialect : bool, optional
+            If True, items in the input using Esri's dialect of WKT
+            will be replaced by OGC standard equivalents.
+
+        Returns
+        -------
+        CRS
+
         """
-        return self.data
+        if isinstance(value, cls):
+            return value
+        elif isinstance(value, int):
+            return cls.from_epsg(value)
+        elif isinstance(value, dict):
+            return cls(**value)
+        elif isinstance(value, string_types):
+            return cls.from_string(value, morph_from_esri_dialect=morph_from_esri_dialect)
+        else:
+            raise CRSError("CRS is invalid: {!r}".format(value))

rasterio/gdal.pxi

@@ -72,6 +72,7 @@ cdef extern from "cpl_vsi.h" nogil:
 
 cdef extern from "ogr_srs_api.h" nogil:
 
+    ctypedef int OGRErr
     ctypedef void * OGRCoordinateTransformationH
     ctypedef void * OGRSpatialReferenceH
 
@@ -84,6 +85,7 @@ cdef extern from "ogr_srs_api.h" nogil:
                      double *y, double *z)
     int OSRAutoIdentifyEPSG(OGRSpatialReferenceH srs)
     int OSRMorphFromESRI(OGRSpatialReferenceH srs)
+    int OSRMorphToESRI(OGRSpatialReferenceH srs)
     void OSRCleanup()
     OGRSpatialReferenceH OSRClone(OGRSpatialReferenceH srs)
     int OSRExportToProj4(OGRSpatialReferenceH srs, char **params)
@@ -93,13 +95,14 @@ cdef extern from "ogr_srs_api.h" nogil:
     const char *OSRGetAuthorityCode(OGRSpatialReferenceH srs, const char *key)
     int OSRImportFromEPSG(OGRSpatialReferenceH srs, int code)
     int OSRImportFromProj4(OGRSpatialReferenceH srs, const char *proj)
+    int OSRImportFromWkt(OGRSpatialReferenceH srs, char **wkt)
     int OSRIsGeographic(OGRSpatialReferenceH srs)
     int OSRIsProjected(OGRSpatialReferenceH srs)
     int OSRIsSame(OGRSpatialReferenceH srs1, OGRSpatialReferenceH srs2)
     OGRSpatialReferenceH OSRNewSpatialReference(const char *wkt)
     void OSRRelease(OGRSpatialReferenceH srs)
     int OSRSetFromUserInput(OGRSpatialReferenceH srs, const char *input)
-
+    OGRErr OSRValidate(OGRSpatialReferenceH srs)
 
 cdef extern from "gdal.h" nogil:
 

rasterio/io.py

@@ -11,7 +11,7 @@ from rasterio._io import (
     DatasetReaderBase, DatasetWriterBase, BufferedDatasetWriterBase,
     MemoryFileBase)
 from rasterio.windows import WindowMethodsMixin
-from rasterio.env import ensure_env
+from rasterio.env import ensure_env, env_ctx_if_needed
 from rasterio.transform import TransformMethodsMixin
 from rasterio.path import UnparsedPath
 
@@ -136,9 +136,12 @@ class MemoryFile(MemoryFileBase):
                           nodata=nodata, **kwargs)
 
     def __enter__(self):
+        self._env = env_ctx_if_needed()
+        self._env.__enter__()
         return self
 
     def __exit__(self, *args, **kwargs):
+        self._env.__exit__()
         self.close()
 
 

rasterio/rio/info.py

@@ -6,7 +6,6 @@ import json
 import click
 
 import rasterio
-import rasterio.crs
 from rasterio.rio import options
 
 
@@ -70,13 +69,13 @@ def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx,
         info['bounds'] = src.bounds
 
         if src.crs:
-            proj4 = src.crs.to_string()
-            if proj4.startswith('+init=epsg'):
-                proj4 = proj4.split('=')[1].upper()
-            info['crs'] = proj4
+            epsg = src.crs.to_epsg()
+            if epsg:
+                info['crs'] = 'EPSG:{}'.format(epsg)
+            else:
+                info['crs'] = src.crs.to_string()
         else:
             info['crs'] = None
-            proj4 = ''
 
         info['res'] = src.res
         info['colorinterp'] = [ci.name for ci in src.colorinterp]
@@ -86,26 +85,30 @@ def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx,
         info['mask_flags'] = [[
             flag.name for flag in flags] for flags in src.mask_flag_enums]
 
-        if proj4 != '':
+        if src.crs:
             info['lnglat'] = src.lnglat()
 
+        gcps, gcps_crs = src.gcps
+
+        if gcps:
+            info['gcps'] = {'points': [p.asdict() for p in gcps]}
+            if crs:
+                epsg = crs.to_epsg()
+                if epsg:
+                    info['gcps']['crs'] = 'EPSG:{}'.format(epsg)
+                else:
+                    info['gcps']['crs'] = src.crs.to_string()
+            else:
+                info['gcps']['crs'] = None
+
         if verbose:
             stats = [{'min': float(b.min()),
                       'max': float(b.max()),
                       'mean': float(b.mean())
                       } for b in src.read(masked=masked)]
             info['stats'] = stats
-
             info['checksum'] = [src.checksum(i) for i in src.indexes]
 
-            gcps, crs = src.gcps
-            proj4 = crs.to_string()
-            if proj4.startswith('+init=epsg'):
-                proj4 = proj4.split('=')[1].upper()
-            if gcps:
-                info['gcps'] = {
-                    'crs': proj4, 'points': [p.asdict() for p in gcps]}
-
         if aspect == 'meta':
             if meta_member == 'subdatasets':
                 for name in src.subdatasets:

rasterio/vrt.py

@@ -6,6 +6,7 @@ import rasterio
 from rasterio._warp import WarpedVRTReaderBase
 from rasterio.dtypes import _gdal_typename
 from rasterio.enums import MaskFlags
+from rasterio.env import env_ctx_if_needed
 from rasterio.path import parse_path, vsi_path
 from rasterio.transform import TransformMethodsMixin
 from rasterio.windows import WindowMethodsMixin
@@ -56,10 +57,13 @@ class WarpedVRT(WarpedVRTReaderBase, WindowMethodsMixin,
             self.closed and 'closed' or 'open', self.name, self.mode)
 
     def __enter__(self):
+        self._env = env_ctx_if_needed()
+        self._env.__enter__()
         self.start()
         return self
 
     def __exit__(self, *args, **kwargs):
+        self._env.__exit__()
         self.close()
 
     def __del__(self):