[rasterio] 01/03: Imported Upstream version 0.30.0

Johan Van de Wauw johanvdw-guest at moszumanska.debian.org
Wed Nov 18 12:56:52 UTC 2015


This is an automated email from the git hooks/post-receive script.

johanvdw-guest pushed a commit to branch master
in repository rasterio.

commit aa9e1e45994f1339acf0021f63f587a22fc86438
Author: Johan Van de Wauw <johan at vandewauw.be>
Date:   Wed Nov 18 12:42:40 2015 +0100

    Imported Upstream version 0.30.0
---
 .coveragerc                      |  13 +
 .gitignore                       |   1 +
 .travis.yml                      |   5 +-
 CHANGES.txt                      |  15 +-
 docs/windowed-rw.rst             |  49 ++
 rasterio/__init__.py             |  87 +++-
 rasterio/_base.pyx               |  25 +-
 rasterio/_features.pxd           |   1 -
 rasterio/_features.pyx           | 158 ++++---
 rasterio/_io.pyx                 | 179 +++++++-
 rasterio/dtypes.py               | 106 ++++-
 rasterio/enums.py                |   7 +
 rasterio/features.py             | 164 +++----
 rasterio/rio/convert.py          |   4 -
 rasterio/rio/features.py         |  25 +-
 rasterio/rio/info.py             |   3 -
 rasterio/rio/main.py             |   1 +
 rasterio/rio/merge.py            |   1 -
 rasterio/rio/options.py          |  29 +-
 rasterio/rio/sample.py           |   4 -
 rasterio/tools/merge.py          |   2 +
 rasterio/transform.py            |   2 +
 rasterio/vfs.py                  |  43 ++
 rasterio/warnings.py             |   9 +
 requirements-dev.txt             |   7 +-
 setup.cfg                        |  10 +-
 setup.py                         |  19 +-
 tests/conftest.py                | 164 +++++++
 tests/data/files.zip             | Bin 0 -> 854723 bytes
 tests/test_band_masks.py         | 108 ++++-
 tests/test_blocks.py             |   4 +-
 tests/test_colormap.py           |   2 -
 tests/test_dtypes.py             |  53 ++-
 tests/test_features.py           | 728 ++++++++++++++++++++++++++++++
 tests/test_features_bounds.py    |  65 ---
 tests/test_features_rasterize.py | 181 --------
 tests/test_features_shapes.py    | 144 ------
 tests/test_features_sieve.py     | 185 --------
 tests/test_indexing.py           | 125 ++++++
 tests/test_mask_creation.py      |   9 +
 tests/test_profile.py            |  31 +-
 tests/test_rio_features.py       | 947 +++++++++++++++++++++++----------------
 tests/test_rio_options.py        |  63 +++
 tests/test_vfs.py                | 101 +++++
 tests/test_warnings.py           |   8 +
 45 files changed, 2644 insertions(+), 1243 deletions(-)

diff --git a/.coveragerc b/.coveragerc
new file mode 100644
index 0000000..d4a1f13
--- /dev/null
+++ b/.coveragerc
@@ -0,0 +1,13 @@
+[run]
+plugins = Cython.Coverage
+source = rasterio
+omit = *.pxd
+branch = True
+
+[report]
+exclude_lines =
+    pragma: no cover
+    def __repr__
+    raise AssertionError
+    raise NotImplementedError
+    if __name__ == .__main__.:
diff --git a/.gitignore b/.gitignore
index 0b87ebd..bd98080 100644
--- a/.gitignore
+++ b/.gitignore
@@ -42,6 +42,7 @@ htmlcov/
 nosetests.xml
 coverage.xml
 *,cover
+rasterio/_*.html
 
 # Translations
 *.mo
diff --git a/.travis.yml b/.travis.yml
index dbe6503..7301728 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -10,10 +10,11 @@ env:
     - PIP_FIND_LINKS=file://$HOME/.cache/pip/wheels
     - GDALINST=$HOME/gdalinstall
     - GDALBUILD=$HOME/gdalbuild
+    - CYTHON_COVERAGE=1
   matrix:
     - GDALVERSION = "1.9.2"
     - GDALVERSION = "1.11.2"
-      #- GDALVERSION = "2.0.1"
+    - GDALVERSION = "2.0.1"
 addons:
   apt:
     packages:
@@ -37,7 +38,7 @@ install:
   - "pip wheel -r requirements-dev.txt"
   - "pip install -r requirements-dev.txt"
   - "pip install --upgrade --force-reinstall --global-option=build_ext --global-option='-I$GDALINST/gdal-$GDALVERSION/include' --global-option='-L$GDALINST/gdal-$GDALVERSION/lib' --global-option='-R$GDALINST/gdal-$GDALVERSION/lib' -e ."
-  - "pip install coveralls"
+  - "pip install coveralls>=1.1"
   - "pip install -e ."
 script:
   - py.test --cov rasterio --cov-report term-missing
diff --git a/CHANGES.txt b/CHANGES.txt
index 4f28341..1b691b5 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,11 +1,24 @@
 Changes
 =======
 
+0.30.0 (2015-11-16)
+-------------------
+- Added window utilities: get_data_window(), window_union(),
+  window_intersection(), windows_intersect() (#496, #506).
+- Warn when an alpha band that might provide a dataset mask is shadowed by a
+  nodata attribute (#508, #523).
+- IPython is not the default interpeter for rio-insp and the documentation 
+  saying it is has been corrected (#518).
+- Guard against creating datasets with block sizes larger than the dataset
+  width and height. Such datasets are semi-broken and are likely to be 
+  mangled when read (#521).
+- Refactor of the `rasterio.features` tests (#522).
+
 0.29.0 (2015-10-22)
 -------------------
 - Fill masked arrays in rio-calc when using Numpy 1.10.x as well as with 1.8.x
   (#500).
-- When a raster dataset is not tiled, blockxszie and blockysize items are no
+- When a raster dataset is not tiled, blockxsize and blockysize items are no
   longer included in its `profile` property. This prevents meaningless block
   size parameters from stripped, not tiled, datasets from being used when
   creating new datasets (#503).
diff --git a/docs/windowed-rw.rst b/docs/windowed-rw.rst
index d91f5a9..fe5d497 100644
--- a/docs/windowed-rw.rst
+++ b/docs/windowed-rw.rst
@@ -161,6 +161,55 @@ This example also demonstrates decimation.
    :width: 500
    :height: 300
 
+
+Data windows
+------------
+
+Sometimes it is desirable to crop off an outer boundary of NODATA values around
+a dataset:
+
+.. code-block:: python
+
+    from rasterio import get_data_window
+
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        window = get_data_window(src.read(1, masked=True))
+        # window = ((3, 714), (13, 770))
+
+        kwargs = src.meta.copy()
+        del kwargs['transform']
+        kwargs.update({
+            'height': window[0][1] - window[0][0],
+            'width': window[1][1] - window[1][0],
+            'affine': src.window_transform(window)
+        })
+
+        with rasterio.open('/tmp/cropped.tif', 'w', **kwargs) as dst:
+            dst.write(src.read(window=window))
+
+
+Window utilities
+----------------
+
+Basic union and intersection operations are available for windows, to streamline
+operations across dynamically created windows for a series of bands or datasets
+with the same full extent.
+
+.. code-block:: python
+
+    from rasterio import window_union, window_intersection
+
+    # Full window is ((0, 1000), (0, 500))
+    window1 = ((100, 500), (10, 500))
+    window2 = ((10, 150), (50, 250))
+
+    outer = window_union([window1, window2])
+    # outer = ((10, 500), (10, 500))
+
+    inner = window_intersection([window1, window2])
+    # inner = ((100, 150), (50, 250))
+
+
 Blocks
 ------
 
diff --git a/rasterio/__init__.py b/rasterio/__init__.py
index 43a7b26..76ff733 100644
--- a/rasterio/__init__.py
+++ b/rasterio/__init__.py
@@ -23,7 +23,7 @@ from rasterio import _err, coords, enums
 
 __all__ = [
     'band', 'open', 'drivers', 'copy', 'pad']
-__version__ = "0.29.0"
+__version__ = "0.30.0"
 
 log = logging.getLogger('rasterio')
 class NullHandler(logging.Handler):
@@ -80,6 +80,11 @@ def open(
             a negative number i.e. -1 * pixel height
     f: Y coordinate of the top left corner of the top left pixel 
 
+    A virtual filesystem can be specified. The ``vfs`` parameter may be
+    an Apache Commons VFS style string beginning with "zip://" or
+    "tar://"". In this case, the ``path`` must be an absolute path
+    within that container.
+
     Finally, additional kwargs are passed to GDAL as driver-specific
     dataset creation parameters.
     """
@@ -167,3 +172,83 @@ def pad(array, transform, pad_width, mode=None, **kwargs):
     padded_trans[2] -= pad_width*padded_trans[0]
     padded_trans[5] -= pad_width*padded_trans[4]
     return padded_array, Affine(*padded_trans[:6])
+
+
+def get_data_window(arr, nodata=None):
+    """
+    Returns a window for the non-nodata pixels within the input array.
+
+    Parameters
+    ----------
+    arr: numpy ndarray, <= 3 dimensions
+    nodata: number
+        If None, will either return a full window if arr is not a masked
+        array, or will use the mask to determine non-nodata pixels.
+        If provided, it must be a number within the valid range of the dtype
+        of the input array.
+
+    Returns
+    -------
+    ((row_start, row_stop), (col_start, col_stop))
+
+    """
+
+    from rasterio._io import get_data_window
+    return get_data_window(arr, nodata)
+
+
+def window_union(windows):
+    """
+    Union windows and return the outermost extent they cover.
+
+    Parameters
+    ----------
+    windows: list-like of window objects
+        ((row_start, row_stop), (col_start, col_stop))
+
+    Returns
+    -------
+    ((row_start, row_stop), (col_start, col_stop))
+    """
+
+    from rasterio._io import window_union
+    return window_union(windows)
+
+
+def window_intersection(windows):
+    """
+    Intersect windows and return the innermost extent they cover.
+
+    Will raise ValueError if windows do not intersect.
+
+    Parameters
+    ----------
+    windows: list-like of window objects
+        ((row_start, row_stop), (col_start, col_stop))
+
+    Returns
+    -------
+    ((row_start, row_stop), (col_start, col_stop))
+    """
+
+    from rasterio._io import window_intersection
+    return window_intersection(windows)
+
+
+def windows_intersect(windows):
+    """
+    Test if windows intersect.
+
+    Parameters
+    ----------
+    windows: list-like of window objects
+        ((row_start, row_stop), (col_start, col_stop))
+
+    Returns
+    -------
+    boolean:
+        True if all windows intersect.
+    """
+
+    from rasterio._io import windows_intersect
+    return windows_intersect(windows)
diff --git a/rasterio/_base.pyx b/rasterio/_base.pyx
index 04f454d..b4e2fb9 100644
--- a/rasterio/_base.pyx
+++ b/rasterio/_base.pyx
@@ -2,8 +2,11 @@
 
 # cython: boundscheck=False
 
+from __future__ import absolute_import
+
 import logging
 import math
+import os
 import sys
 import warnings
 
@@ -16,6 +19,7 @@ from rasterio import dtypes
 from rasterio.coords import BoundingBox
 from rasterio.transform import Affine
 from rasterio.enums import ColorInterp, Compression, Interleaving
+from rasterio.vfs import parse_path, vsi_path
 
 
 log = logging.getLogger('rasterio')
@@ -62,7 +66,10 @@ cdef class DatasetReader(object):
             self.env = GDALEnv(False)
         self.env.start()
 
-        name_b = self.name.encode('utf-8')
+        path, archive, scheme = parse_path(self.name)
+        path = vsi_path(path, archive=archive, scheme=scheme)
+
+        name_b = path.encode('utf-8')
         cdef const char *fname = name_b
         with cpl_errs:
             self._hds = _gdal.GDALOpen(fname, 0)
@@ -196,7 +203,11 @@ cdef class DatasetReader(object):
         if self._hds == NULL:
             raise ValueError("Null dataset")
         cdef double gt[6]
-        _gdal.GDALGetGeoTransform(self._hds, gt)
+        err = _gdal.GDALGetGeoTransform(self._hds, gt)
+        if err:
+            warnings.warn("GDALGetGeoTransform failed, default invalid "
+                          "transform will be returned.")
+
         transform = [0]*6
         for i in range(6):
             transform[i] = gt[i]
@@ -316,6 +327,16 @@ cdef class DatasetReader(object):
         def __get__(self):
             return self.nodatavals[0]
 
+    property mask_flags:
+        """Mask flags for each band."""
+
+        def __get__(self):
+            flags = [0]*self.count
+            for i, j in zip(range(self.count), self.indexes):
+                hband = _gdal.GDALGetRasterBand(self._hds, j)
+                flags[i] = _gdal.GDALGetMaskFlags(hband)
+            return flags
+
     def block_windows(self, bidx=0):
         """Returns an iterator over a band's block windows and their
         indexes.
diff --git a/rasterio/_features.pxd b/rasterio/_features.pxd
index e4898ae..1ca21ba 100644
--- a/rasterio/_features.pxd
+++ b/rasterio/_features.pxd
@@ -12,7 +12,6 @@ cdef class GeomBuilder:
     cpdef _buildPolygon(self)
     cpdef _buildMultiPolygon(self)
     cdef build(self, void *geom)
-    cpdef build_wkb(self, object wkb)
 
 
 cdef class OGRGeomBuilder:
diff --git a/rasterio/_features.pyx b/rasterio/_features.pyx
index 0289303..89f9a3d 100644
--- a/rasterio/_features.pyx
+++ b/rasterio/_features.pyx
@@ -1,6 +1,5 @@
 
 import logging
-import json
 import numpy as np
 cimport numpy as np
 from rasterio._io cimport InMemoryRaster
@@ -62,10 +61,16 @@ def _shapes(image, mask, connectivity, transform):
     cdef InMemoryRaster mem_ds = None
     cdef InMemoryRaster mask_ds = None
     cdef bint is_float = np.dtype(image.dtype).kind == 'f'
-    cdef int fieldtp = 0
+    cdef int fieldtp = 2 if is_float else 0
 
-    if is_float:
-        fieldtp = 2
+    valid_dtypes = ('int16', 'int32', 'uint8', 'uint16', 'float32')
+
+    if np.dtype(image.dtype).name not in valid_dtypes:
+        raise ValueError('image dtype must be one of: %s'
+                         % (', '.join(valid_dtypes)))
+
+    if connectivity not in (4, 8):
+        raise ValueError("Connectivity Option must be 4 or 8")
 
     if dtypes.is_ndarray(image):
         mem_ds = InMemoryRaster(image, transform)
@@ -76,17 +81,22 @@ def _shapes(image, mask, connectivity, transform):
     else:
         raise ValueError("Invalid source image")
 
-    if dtypes.is_ndarray(mask):
-        # A boolean mask must be converted to uint8 for GDAL
-        mask_ds = InMemoryRaster(mask.astype('uint8'), transform)
-        hmaskband = mask_ds.band
-    elif isinstance(mask, tuple):
+    if mask is not None:
         if mask.shape != image.shape:
             raise ValueError("Mask must have same shape as image")
-        mrdr = mask.ds
-        hmaskband = mrdr.band(mask.bidx)
-    else:
-        hmaskband = NULL
+
+        if np.dtype(mask.dtype).name not in ('bool', 'uint8'):
+            raise ValueError("Mask must be dtype rasterio.bool_ or "
+                             "rasterio.uint8")
+
+        if dtypes.is_ndarray(mask):
+            # A boolean mask must be converted to uint8 for GDAL
+            mask_ds = InMemoryRaster(mask.astype('uint8'), transform)
+            hmaskband = mask_ds.band
+
+        elif isinstance(mask, tuple):
+            mrdr = mask.ds
+            hmaskband = mrdr.band(mask.bidx)
 
     # Create an in-memory feature store.
     hfdriver = _ogr.OGRGetDriverByName("Memory")
@@ -133,7 +143,7 @@ def _shapes(image, mask, connectivity, transform):
         _gdal.CSLDestroy(options)
 
 
-def _sieve(image, size, output, mask, connectivity):
+def _sieve(image, size, out, mask, connectivity):
     """
     Replaces small polygons in `image` with the value of their largest
     neighbor.  Polygons are found for each set of neighboring pixels of the
@@ -147,7 +157,7 @@ def _sieve(image, size, output, mask, connectivity):
         rasterio.uint16, or rasterio.float32.
     size : int
         minimum polygon size (number of pixels) to retain.
-    output : numpy ndarray
+    out : numpy ndarray
         Array of same shape and data type as `image` in which to store results.
     mask : numpy ndarray or rasterio Band object
         Values of False or 0 will be excluded from feature generation.
@@ -168,6 +178,29 @@ def _sieve(image, size, output, mask, connectivity):
     cdef _io.RasterUpdater udr
     cdef _io.RasterReader mask_reader
 
+    valid_dtypes = ('int16', 'int32', 'uint8', 'uint16')
+
+    if np.dtype(image.dtype).name not in valid_dtypes:
+        valid_types_str = ', '.join(('rasterio.{0}'.format(t) for t
+                                     in valid_dtypes))
+        raise ValueError('image dtype must be one of: %s' % valid_types_str)
+
+    if size <= 0:
+        raise ValueError('size must be greater than 0')
+    elif type(size) == float:
+        raise ValueError('size must be an integer number of pixels')
+    elif size > (image.shape[0] * image.shape[1]):
+        raise ValueError('size must be smaller than size of image')
+
+    if connectivity not in (4, 8):
+        raise ValueError('connectivity must be 4 or 8')
+
+    if out.shape != image.shape:
+        raise ValueError('out raster shape must be same as image shape')
+
+    if np.dtype(image.dtype).name != np.dtype(out.dtype).name:
+        raise ValueError('out raster must match dtype of image')
+
     if dtypes.is_ndarray(image):
         in_mem_ds = InMemoryRaster(image)
         in_band = in_mem_ds.band
@@ -177,27 +210,32 @@ def _sieve(image, size, output, mask, connectivity):
     else:
         raise ValueError("Invalid source image")
 
-    if dtypes.is_ndarray(output):
-        log.debug("Output array: %r", output)
-        out_mem_ds = InMemoryRaster(output)
+    if dtypes.is_ndarray(out):
+        log.debug("out array: %r", out)
+        out_mem_ds = InMemoryRaster(out)
         out_band = out_mem_ds.band
-    elif isinstance(output, tuple):
-        udr = output.ds
-        out_band = udr.band(output.bidx)
+    elif isinstance(out, tuple):
+        udr = out.ds
+        out_band = udr.band(out.bidx)
     else:
-        raise ValueError("Invalid output image")
+        raise ValueError("Invalid out image")
 
-    if dtypes.is_ndarray(mask):
-        # A boolean mask must be converted to uint8 for GDAL
-        mask_mem_ds = InMemoryRaster(mask.astype('uint8'))
-        mask_band = mask_mem_ds.band
-    elif isinstance(mask, tuple):
+    if mask is not None:
         if mask.shape != image.shape:
             raise ValueError("Mask must have same shape as image")
-        mask_reader = mask.ds
-        mask_band = mask_reader.band(mask.bidx)
-    else:
-        mask_band = NULL
+
+        if np.dtype(mask.dtype) not in ('bool', 'uint8'):
+            raise ValueError("Mask must be dtype rasterio.bool_ or "
+                             "rasterio.uint8")
+
+        if dtypes.is_ndarray(mask):
+            # A boolean mask must be converted to uint8 for GDAL
+            mask_mem_ds = InMemoryRaster(mask.astype('uint8'))
+            mask_band = mask_mem_ds.band
+
+        elif isinstance(mask, tuple):
+            mask_reader = mask.ds
+            mask_band = mask_reader.band(mask.bidx)
 
     _gdal.GDALSieveFilter(
         in_band,
@@ -210,8 +248,8 @@ def _sieve(image, size, output, mask, connectivity):
         NULL
     )
 
-    # Read from out_band into output
-    _io.io_auto(output, out_band, False)
+    # Read from out_band into out
+    _io.io_auto(out, out_band, False)
 
     if in_mem_ds is not None:
         in_mem_ds.close()
@@ -238,7 +276,7 @@ def _rasterize(shapes, image, transform, all_touched):
     all_touched : boolean, optional
         If True, all pixels touched by geometries will be burned in.
         If false, only pixels whose center is within the polygon or
-        that are selected by brezenhams line algorithm will be burned
+        that are selected by Bresenham's line algorithm will be burned
         in.
     """
 
@@ -310,24 +348,21 @@ def _bounds(geometry):
     TODO: add to Fiona.
     """
 
-    try:  
-        if 'features' in geometry:
-            xmins = []
-            ymins = []
-            xmaxs = []
-            ymaxs = []
-            for feature in geometry['features']:
-                xmin, ymin, xmax, ymax = _bounds(feature['geometry'])
-                xmins.append(xmin)
-                ymins.append(ymin)
-                xmaxs.append(xmax)
-                ymaxs.append(ymax)
-            return min(xmins), min(ymins), max(xmaxs), max(ymaxs)
-        else:
-            xyz = tuple(zip(*list(_explode(geometry['coordinates']))))
-            return min(xyz[0]), min(xyz[1]), max(xyz[0]), max(xyz[1])
-    except (KeyError, TypeError):
-        return None
+    if 'features' in geometry:
+        xmins = []
+        ymins = []
+        xmaxs = []
+        ymaxs = []
+        for feature in geometry['features']:
+            xmin, ymin, xmax, ymax = _bounds(feature['geometry'])
+            xmins.append(xmin)
+            ymins.append(ymin)
+            xmaxs.append(xmax)
+            ymaxs.append(ymax)
+        return min(xmins), min(ymins), max(xmaxs), max(ymaxs)
+    else:
+        xyz = tuple(zip(*list(_explode(geometry['coordinates']))))
+        return min(xyz[0]), min(xyz[1]), max(xyz[0]), max(xyz[1])
 
 
 # Mapping of OGR integer geometry types to GeoJSON type names.
@@ -359,16 +394,6 @@ GEOJSON2OGR_GEOMETRY_TYPES = dict(
 
 # Geometry related functions and classes follow.
 
-cdef void * _createOgrGeomFromWKB(object wkb) except NULL:
-    """Make an OGR geometry from a WKB string"""
-
-    geom_type = bytearray(wkb)[1]
-    cdef unsigned char *buffer = wkb
-    cdef void *cogr_geometry = _ogr.OGR_G_CreateGeometry(geom_type)
-    if cogr_geometry != NULL:
-        _ogr.OGR_G_ImportFromWkb(cogr_geometry, buffer, len(wkb))
-    return cogr_geometry
-
 
 cdef _deleteOgrGeom(void *cogr_geometry):
     """Delete an OGR geometry"""
@@ -445,15 +470,6 @@ cdef class GeomBuilder:
         self.geom = geom
         return getattr(self, '_build' + self.geomtypename)()
 
-    cpdef build_wkb(self, object wkb):
-        """Builds a GeoJSON object from a Well-Known Binary format (WKB)."""
-        # The only other method anyone needs to call
-        cdef object data = wkb
-        cdef void *cogr_geometry = _createOgrGeomFromWKB(data)
-        result = self.build(cogr_geometry)
-        _deleteOgrGeom(cogr_geometry)
-        return result
-
 
 cdef class OGRGeomBuilder:
     """
diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx
index 1d42712..070d0fd 100644
--- a/rasterio/_io.pyx
+++ b/rasterio/_io.pyx
@@ -1,5 +1,7 @@
 # cython: boundscheck=False
 
+from __future__ import absolute_import
+
 import logging
 import math
 import os
@@ -20,8 +22,10 @@ from rasterio import dtypes
 from rasterio.coords import BoundingBox
 from rasterio.five import text_type, string_types
 from rasterio.transform import Affine
-from rasterio.enums import ColorInterp, Resampling
+from rasterio.enums import ColorInterp, MaskFlags, Resampling
 from rasterio.sample import sample_gen
+from rasterio.vfs import parse_path
+from rasterio.warnings import NodataShadowWarning
 
 
 log = logging.getLogger('rasterio')
@@ -36,6 +40,7 @@ else:
 
 log.addHandler(NullHandler())
 
+
 cdef bint in_dtype_range(value, dtype):
     """Returns True if value is in the range of dtype, else False."""
     infos = {
@@ -1059,6 +1064,12 @@ cdef class RasterReader(_base.DatasetReader):
         gdt = dtypes.dtype_rev[dtype]
 
         if masks:
+            # Warn if nodata attribute is shadowing an alpha band.
+            if self.count == 4 and self.colorinterp(4) == ColorInterp.alpha:
+                for flags in self.mask_flags:
+                    if flags & MaskFlags.nodata:
+                        warnings.warn(NodataShadowWarning())
+
             retval = io_multi_mask(
                             self._hds, 0, xoff, yoff, width, height,
                             out, indexes_arr, indexes_count)
@@ -1249,8 +1260,7 @@ cdef class RasterUpdater(RasterReader):
         cdef void *drv = NULL
         cdef void *hband = NULL
         cdef int success
-        name_b = self.name.encode('utf-8')
-        cdef const char *fname = name_b
+
 
         # Is there not a driver manager already?
         if driver_count() == 0 and not self.env:
@@ -1259,7 +1269,16 @@ cdef class RasterUpdater(RasterReader):
         else:
             self.env = GDALEnv(False)
         self.env.start()
-        
+
+        path, archive, scheme = parse_path(self.name)
+        if scheme and scheme != 'file':
+            raise TypeError(
+                "VFS '{0}' datasets can not be created or updated.".format(
+                    scheme))
+
+        name_b = path.encode('utf-8')
+        cdef const char *fname = name_b
+
         kwds = []
 
         if self.mode == 'w':
@@ -1269,8 +1288,8 @@ cdef class RasterUpdater(RasterReader):
                 raise ValueError("only GTiffs can be opened in 'w' mode")
 
             # Delete existing file, create.
-            if os.path.exists(self.name):
-                os.unlink(self.name)
+            if os.path.exists(path):
+                os.unlink(path)
             
             driver_b = self.driver.encode('utf-8')
             drv_name = driver_b
@@ -1298,6 +1317,13 @@ cdef class RasterUpdater(RasterReader):
                     continue
                 kwds.append((k.lower(), v))
                 k, v = k.upper(), str(v).upper()
+
+                # Guard against block size that exceed image size.
+                if k == 'BLOCKXSIZE' and int(v) > self.width:
+                    raise ValueError("blockxsize exceeds raster width.")
+                if k == 'BLOCKYSIZE' and int(v) > self.height:
+                    raise ValueError("blockysize exceeds raster height.")
+
                 key_b = k.encode('utf-8')
                 val_b = v.encode('utf-8')
                 key_c = key_b
@@ -1310,6 +1336,7 @@ cdef class RasterUpdater(RasterReader):
             self._hds = _gdal.GDALCreate(
                 drv, fname, self.width, self.height, self._count,
                 gdal_dtype, options)
+
             if self._hds == NULL:
                 raise ValueError("NULL dataset")
 
@@ -1426,7 +1453,9 @@ cdef class RasterUpdater(RasterReader):
         cdef double gt[6]
         for i in range(6):
             gt[i] = transform[i]
-        retval = _gdal.GDALSetGeoTransform(self._hds, gt)
+        err = _gdal.GDALSetGeoTransform(self._hds, gt)
+        if err:
+            raise ValueError("transform not set: %s" % transform)
         self._transform = transform
 
     property transform:
@@ -2002,6 +2031,12 @@ def writer(path, mode, **kwargs):
     cdef const char *drv_name = NULL
     cdef const char *fname = NULL
 
+    path, archive, scheme = parse_path(path)
+    if scheme and scheme != 'file':
+        raise TypeError(
+            "VFS '{0}' datasets can not be created or updated.".format(
+                scheme))
+
     if mode == 'w' and 'driver' in kwargs:
         if kwargs['driver'] == 'GTiff':
             return RasterUpdater(path, mode, **kwargs)
@@ -2041,3 +2076,133 @@ def virtual_file_to_buffer(filename):
     log.debug("Buffer length: %d bytes", n)
     cdef np.uint8_t[:] buff_view = <np.uint8_t[:n]>buff
     return buff_view
+
+
+def get_data_window(arr, nodata=None):
+    """
+    Returns a window for the non-nodata pixels within the input array.
+
+    Parameters
+    ----------
+    arr: numpy ndarray, <= 3 dimensions
+    nodata: number
+        If None, will either return a full window if arr is not a masked
+        array, or will use the mask to determine non-nodata pixels.
+        If provided, it must be a number within the valid range of the dtype
+        of the input array.
+
+    Returns
+    -------
+    ((row_start, row_stop), (col_start, col_stop))
+
+    """
+
+    num_dims = len(arr.shape)
+    if num_dims > 3:
+        raise ValueError('get_data_window input array must have no more than '
+                         '3 dimensions')
+
+    if nodata is None:
+        if not hasattr(arr, 'mask'):
+            return ((0, arr.shape[-2]), (0, arr.shape[-1]))
+    else:
+        arr = np.ma.masked_array(arr, arr == nodata)
+
+    if num_dims == 2:
+        data_rows, data_cols = np.where(arr.mask == False)
+    else:
+        data_rows, data_cols = np.where(
+            np.any(np.rollaxis(arr.mask, 0, 3) == False, axis=2)
+        )
+
+    if data_rows.size:
+        row_range = (data_rows.min(), data_rows.max() + 1)
+    else:
+        row_range = (0, 0)
+
+    if data_cols.size:
+        col_range = (data_cols.min(), data_cols.max() + 1)
+    else:
+        col_range = (0, 0)
+
+    return (row_range, col_range)
+
+
+def window_union(windows):
+    """
+    Union windows and return the outermost extent they cover.
+
+    Parameters
+    ----------
+    windows: list-like of window objects
+        ((row_start, row_stop), (col_start, col_stop))
+
+    Returns
+    -------
+    ((row_start, row_stop), (col_start, col_stop))
+    """
+
+
+    stacked = np.dstack(windows)
+    return (
+        (stacked[0, 0].min(), stacked[0, 1].max()),
+        (stacked[1, 0].min(), stacked[1, 1]. max())
+    )
+
+
+def window_intersection(windows):
+    """
+    Intersect windows and return the innermost extent they cover.
+
+    Will raise ValueError if windows do not intersect.
+
+    Parameters
+    ----------
+    windows: list-like of window objects
+        ((row_start, row_stop), (col_start, col_stop))
+
+    Returns
+    -------
+    ((row_start, row_stop), (col_start, col_stop))
+    """
+
+    if not windows_intersect(windows):
+        raise ValueError('windows do not intersect')
+
+    stacked = np.dstack(windows)
+    return (
+        (stacked[0, 0].max(), stacked[0, 1].min()),
+        (stacked[1, 0].max(), stacked[1, 1]. min())
+    )
+
+
+def windows_intersect(windows):
+    """
+    Test if windows intersect.
+
+    Parameters
+    ----------
+    windows: list-like of window objects
+        ((row_start, row_stop), (col_start, col_stop))
+
+    Returns
+    -------
+    boolean:
+        True if all windows intersect.
+    """
+
+    from itertools import combinations
+
+    def intersects(range1, range2):
+        return not (
+            range1[0] > range2[1] or range1[1] < range2[0]
+        )
+
+    windows = np.array(windows)
+
+    for i in (0, 1):
+        for c in combinations(windows[:, i], 2):
+            if not intersects(*c):
+                return False
+
+    return True
diff --git a/rasterio/dtypes.py b/rasterio/dtypes.py
index 90097d8..8b2de80 100644
--- a/rasterio/dtypes.py
+++ b/rasterio/dtypes.py
@@ -74,6 +74,7 @@ def _gdal_typename(dt):
     except KeyError:
         return typename_fwd[dtype_rev[dt().dtype.name]]
 
+
 def check_dtype(dt):
     if dt not in dtype_rev:
         try:
@@ -83,32 +84,105 @@ def check_dtype(dt):
     return True
 
 
-def get_minimum_int_dtype(values):
+def get_minimum_dtype(values):
     """
-    Uses range checking to determine the minimum integer data type required
+    Uses range checking to determine the minimum integer or floating point
+    data type required
     to represent values.
 
-    :param values: numpy array
-    :return: named data type that can be later used to create a numpy dtype
+    Parameters
+    ----------
+    values: list-like
+
+
+    Returns
+    -------
+    rasterio dtype string
     """
 
+    import numpy
+
+    if not is_ndarray(values):
+        values = numpy.array(values)
+
     min_value = values.min()
     max_value = values.max()
-    
-    if min_value >= 0:
-        if max_value <= 255:
-            return uint8
-        elif max_value <= 65535:
-            return uint16
-        elif max_value <= 4294967295:
-            return uint32
-    elif min_value >= -32768 and max_value <= 32767:
-        return int16
-    elif min_value >= -2147483648 and max_value <= 2147483647:
-        return int32
+
+    if values.dtype.kind == 'i':
+        if min_value >= 0:
+            if max_value <= 255:
+                return uint8
+            elif max_value <= 65535:
+                return uint16
+            elif max_value <= 4294967295:
+                return uint32
+        elif min_value >= -32768 and max_value <= 32767:
+            return int16
+        elif min_value >= -2147483648 and max_value <= 2147483647:
+            return int32
+
+    else:
+        if min_value >= -3.4028235e+38 and max_value <= 3.4028235e+38:
+            return float32
+        return float64
 
 
 def is_ndarray(array):
     import numpy
 
     return isinstance(array, numpy.ndarray) or hasattr(array, '__array__')
+
+
+def can_cast_dtype(values, dtype):
+    """
+    Tests if values can be cast to dtype without loss of information.
+
+    Parameters
+    ----------
+    values: list-like
+    dtype: numpy dtype or string
+
+    Returns
+    -------
+    boolean
+        True if values can be cast to data type.
+    """
+
+    import numpy
+
+    if not is_ndarray(values):
+        values = numpy.array(values)
+
+    if values.dtype.name == numpy.dtype(dtype).name:
+        return True
+
+    elif values.dtype.kind == 'f':
+        return numpy.allclose(values, values.astype(dtype))
+
+    else:
+        return numpy.array_equal(values, values.astype(dtype))
+
+
+def validate_dtype(values, valid_dtypes):
+    """
+    Tests if dtype of values is one of valid_dtypes.
+
+    Parameters
+    ----------
+    values: list-like
+    valid_dtypes: list-like
+        list of valid dtype strings, e.g., ('int16', 'int32')
+
+    Returns
+    -------
+    boolean:
+        True if dtype of values is one of valid_dtypes
+    """
+
+    import numpy
+
+    if not is_ndarray(values):
+        values = numpy.array(values)
+
+    return (values.dtype.name in valid_dtypes or
+            get_minimum_dtype(values) in valid_dtypes)
\ No newline at end of file
diff --git a/rasterio/enums.py b/rasterio/enums.py
index 3e347a5..c611052 100644
--- a/rasterio/enums.py
+++ b/rasterio/enums.py
@@ -45,3 +45,10 @@ class Compression(Enum):
 class Interleaving(Enum):
     pixel='PIXEL'
     band='BAND'
+
+
+class MaskFlags(IntEnum):
+    all_valid=1
+    per_dataset=2
+    alpha=4
+    nodata=8
diff --git a/rasterio/features.py b/rasterio/features.py
index d52da1a..e49aa76 100644
--- a/rasterio/features.py
+++ b/rasterio/features.py
@@ -1,5 +1,7 @@
 """Functions for working with features in a raster dataset."""
 
+from __future__ import absolute_import
+
 import json
 import logging
 import time
@@ -10,7 +12,7 @@ import numpy as np
 import rasterio
 from rasterio._features import _shapes, _sieve, _rasterize, _bounds
 from rasterio.transform import IDENTITY, guard_transform
-from rasterio.dtypes import get_minimum_int_dtype
+from rasterio.dtypes import validate_dtype, can_cast_dtype, get_minimum_dtype
 
 
 log = logging.getLogger('rasterio')
@@ -104,18 +106,6 @@ def shapes(image, mask=None, connectivity=4, transform=IDENTITY):
 
     """
 
-    valid_dtypes = ('int16', 'int32', 'uint8', 'uint16', 'float32')
-
-    if np.dtype(image.dtype).name not in valid_dtypes:
-        raise ValueError('image dtype must be one of: %s'
-                         % (', '.join(valid_dtypes)))
-
-    if mask is not None and np.dtype(mask.dtype).name not in ('bool', 'uint8'):
-        raise ValueError("Mask must be dtype rasterio.bool_ or rasterio.uint8")
-
-    if connectivity not in (4, 8):
-        raise ValueError("Connectivity Option must be 4 or 8")
-
     transform = guard_transform(transform)
 
     with rasterio.drivers():
@@ -165,49 +155,18 @@ def sieve(image, size, out=None, output=None, mask=None, connectivity=4):
 
     """
 
-    valid_dtypes = ('int16', 'int32', 'uint8', 'uint16')
-
-    if np.dtype(image.dtype).name not in valid_dtypes:
-        valid_types_str = ', '.join(('rasterio.{0}'.format(t) for t
-                                     in valid_dtypes))
-        raise ValueError('image dtype must be one of: %s' % valid_types_str)
-
-    if size <= 0:
-        raise ValueError('size must be greater than 0')
-    elif type(size) == float:
-        raise ValueError('size must be an integer number of pixels')
-    elif size > (image.shape[0] * image.shape[1]):
-        raise ValueError('size must be smaller than size of image')
-
-    if connectivity not in (4, 8):
-        raise ValueError('connectivity must be 4 or 8')
-
-    if mask is not None:
-        if np.dtype(mask.dtype) not in ('bool', 'uint8'):
-            raise ValueError('Mask must be dtype rasterio.bool_ or '
-                             'rasterio.uint8')
-        elif mask.shape != image.shape:
-            raise ValueError('mask shape must be same as image shape')
-
     # Start moving users over to 'out'.
     if output is not None:
         warnings.warn(
             "The 'output' keyword arg has been superceded by 'out' "
             "and will be removed before Rasterio 1.0.",
             FutureWarning,
-            stacklevel=2)
+            stacklevel=2)  # pragma: no cover
     
     out = out if out is not None else output
+
     if out is None:
-        if isinstance(image, tuple):
-            out = np.zeros(image.shape, image.dtype)
-        else:
-            out = np.zeros_like(image)
-    else:
-        if np.dtype(image.dtype).name != np.dtype(out.dtype).name:
-            raise ValueError('out raster must match dtype of image')
-        elif out.shape != image.shape:
-            raise ValueError('out raster shape must be same as image shape')
+        out = np.zeros(image.shape, image.dtype)
 
     with rasterio.drivers():
         _sieve(image, size, out, mask, connectivity)
@@ -268,97 +227,92 @@ def rasterize(
 
     """
 
-    valid_dtypes = ('int16', 'int32', 'uint8', 'uint16', 'uint32', 'float32',
-                    'float64')
-
-    def get_valid_dtype(values):
-        values_dtype = values.dtype
-        if values_dtype.kind == 'i':
-            values_dtype = np.dtype(get_minimum_int_dtype(values))
-        if values_dtype.name in valid_dtypes:
-            return values_dtype
-        return None
-
-    def can_cast_dtype(values, dtype):
-        if values.dtype.name == np.dtype(dtype).name:
-            return True
-        elif values.dtype.kind == 'f':
-            return np.allclose(values, values.astype(dtype))
-        else:
-            return np.array_equal(values, values.astype(dtype))
+    valid_dtypes = (
+        'int16', 'int32', 'uint8', 'uint16', 'uint32', 'float32', 'float64'
+    )
+
+    def format_invalid_dtype(param):
+        return '{0} dtype must be one of: {1}'.format(
+            param, ', '.join(valid_dtypes)
+        )
+
+    def format_cast_error(param, dtype):
+        return '{0} cannot be cast to specified dtype: {1}'.format(param, dtype)
+
 
     if fill != 0:
         fill_array = np.array([fill])
-        if get_valid_dtype(fill_array) is None:
-            raise ValueError('fill must be one of these types: %s'
-                             % (', '.join(valid_dtypes)))
-        elif dtype is not None and not can_cast_dtype(fill_array, dtype):
-            raise ValueError('fill value cannot be cast to specified dtype')
+        if not validate_dtype(fill_array, valid_dtypes):
+            raise ValueError(format_invalid_dtype('fill'))
+
+        if dtype is not None and not can_cast_dtype(fill_array, dtype):
+            raise ValueError(format_cast_error('fill', dtype))
 
     if default_value != 1:
         default_value_array = np.array([default_value])
-        if get_valid_dtype(default_value_array) is None:
-            raise ValueError('default_value must be one of these types: %s'
-                             % (', '.join(valid_dtypes)))
-        elif dtype is not None and not can_cast_dtype(default_value_array,
-                                                      dtype):
-            raise ValueError('default_value cannot be cast to specified dtype')
+        if not validate_dtype(default_value_array, valid_dtypes):
+            raise ValueError(format_invalid_dtype('default_value'))
+
+        if dtype is not None and not can_cast_dtype(default_value_array, dtype):
+            raise ValueError(format_cast_error('default_vaue', dtype))
+
+    if dtype is not None and np.dtype(dtype).name not in valid_dtypes:
+        raise ValueError(format_invalid_dtype('dtype'))
+
 
     valid_shapes = []
     shape_values = []
     for index, item in enumerate(shapes):
-        try:
-            if isinstance(item, (tuple, list)):
-                geom, value = item
-            else:
-                geom = item
-                value = default_value
-            geom = getattr(geom, '__geo_interface__', None) or geom
-            if (not isinstance(geom, dict) or
-                'type' not in geom or 'coordinates' not in geom):
-                raise ValueError(
-                    'Object %r at index %d is not a geometry object' %
-                    (geom, index))
+        if isinstance(item, (tuple, list)):
+            geom, value = item
+        else:
+            geom = item
+            value = default_value
+        geom = getattr(geom, '__geo_interface__', None) or geom
+
+        #not isinstance(geom, dict) or
+        if 'type' in geom or 'coordinates' in geom:
             valid_shapes.append((geom, value))
             shape_values.append(value)
-        except Exception:
-            log.exception('Exception caught, skipping shape %d', index)
+
+        else:
+            raise ValueError(
+                'Invalid geometry object at index {0}'.format(index)
+            )
 
     if not valid_shapes:
-        raise ValueError('No valid shapes found for rasterize.  Shapes must be '
-                         'valid geometry objects')
+        raise ValueError('No valid geometry objects found for rasterize')
 
     shape_values = np.array(shape_values)
-    values_dtype = get_valid_dtype(shape_values)
-    if values_dtype is None:
-        raise ValueError('shape values must be one of these dtypes: %s' %
-                         (', '.join(valid_dtypes)))
+
+    if not validate_dtype(shape_values, valid_dtypes):
+        raise ValueError(format_invalid_dtype('shape values'))
 
     if dtype is None:
-        dtype = values_dtype
-    elif np.dtype(dtype).name not in valid_dtypes:
-        raise ValueError('dtype must be one of: %s' % (', '.join(valid_dtypes)))
+        dtype = get_minimum_dtype(np.append(shape_values, fill))
+
     elif not can_cast_dtype(shape_values, dtype):
-        raise ValueError('shape values could not be cast to specified dtype')
+        raise ValueError(format_cast_error('shape values', dtype))
 
     if output is not None:
         warnings.warn(
             "The 'output' keyword arg has been superceded by 'out' "
             "and will be removed before Rasterio 1.0.",
             FutureWarning,
-            stacklevel=2)
+            stacklevel=2) # pragma: no cover
+
     out = out if out is not None else output
     if out is not None:
         if np.dtype(out.dtype).name not in valid_dtypes:
-            raise ValueError('Output image dtype must be one of: %s'
-                             % (', '.join(valid_dtypes)))
+            raise ValueError(format_invalid_dtype('out'))
+
         if not can_cast_dtype(shape_values, out.dtype):
-            raise ValueError('shape values cannot be cast to dtype of output '
-                             'image')
+            raise ValueError(format_cast_error('shape values', out.dtype.name))
 
     elif out_shape is not None:
         out = np.empty(out_shape, dtype=dtype)
         out.fill(fill)
+
     else:
         raise ValueError('Either an output shape or image must be provided')
 
diff --git a/rasterio/rio/convert.py b/rasterio/rio/convert.py
index 11781f7..017b2c9 100644
--- a/rasterio/rio/convert.py
+++ b/rasterio/rio/convert.py
@@ -1,7 +1,6 @@
 """File translation command"""
 
 import logging
-import warnings
 
 import click
 from cligj import format_opt
@@ -13,9 +12,6 @@ import rasterio
 from rasterio.coords import disjoint_bounds
 
 
-warnings.simplefilter('default')
-
-
 # Clip command
 @click.command(short_help='Clip a raster to given bounds.')
 @click.argument(
diff --git a/rasterio/rio/features.py b/rasterio/rio/features.py
index 0345bbb..ebfcf19 100644
--- a/rasterio/rio/features.py
+++ b/rasterio/rio/features.py
@@ -122,14 +122,26 @@ def mask(
         bounds = geojson.get('bbox', calculate_bounds(geojson))
 
         with rasterio.open(input) as src:
-            has_disjoint_bounds = disjoint_bounds(bounds, src.bounds)
+            # If y pixel value is positive, then invert y dimension in bounds
+            invert_y = src.affine.e > 0
+
+            src_bounds = src.bounds
+            if invert_y:
+                src_bounds = [src.bounds[0], src.bounds[3],
+                              src.bounds[2], src.bounds[1]]
+
+            has_disjoint_bounds = disjoint_bounds(bounds, src_bounds)
 
             if crop:
                 if has_disjoint_bounds:
+
                     raise click.BadParameter('not allowed for GeoJSON outside '
                                              'the extent of the input raster',
                                              param=crop, param_hint='--crop')
 
+                if invert_y:
+                    bounds = (bounds[0], bounds[3], bounds[2], bounds[1])
+
                 window = src.window(*bounds)
                 transform = src.window_transform(window)
                 (r1, r2), (c1, c2) = window
@@ -170,7 +182,7 @@ def mask(
 
 # Shapes command.
 @click.command(short_help="Write shapes extracted from bands or masks.")
- at click.argument('input', type=click.Path(exists=True))
+ at options.file_in_arg
 @options.output_opt
 @precision_opt
 @indent_opt
@@ -256,6 +268,9 @@ def shapes(
 
         def __call__(self):
             with rasterio.open(input) as src:
+                if bidx is not None and bidx > src.count:
+                    raise ValueError('bidx is out of range for raster')
+
                 img = None
                 msk = None
 
@@ -533,6 +548,12 @@ def rasterize(
 
                 kwargs = template_ds.meta.copy()
                 kwargs['count'] = 1
+
+                # DEPRECATED
+                # upgrade transform to affine object or we may get an invalid
+                # transform set on output
+                kwargs['transform'] = template_ds.affine
+
                 template_ds.close()
 
             else:
diff --git a/rasterio/rio/info.py b/rasterio/rio/info.py
index 71f5a37..fa390dd 100644
--- a/rasterio/rio/info.py
+++ b/rasterio/rio/info.py
@@ -274,7 +274,6 @@ def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx,
     verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
     logger = logging.getLogger('rio')
     mode = 'r' if (verbose or meta_member == 'stats') else 'r-'
-
     try:
         with rasterio.drivers(CPL_DEBUG=(verbosity > 2)):
             with rasterio.open(input, mode) as src:
@@ -334,8 +333,6 @@ def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx,
 @click.pass_context
 def insp(ctx, input, mode, interpreter):
     """ Open the input file in a Python interpreter.
-
-    IPython will be used as the default interpreter, if available.
     """
     import rasterio.tool
     verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
diff --git a/rasterio/rio/main.py b/rasterio/rio/main.py
index 2282ae5..28bfb60 100644
--- a/rasterio/rio/main.py
+++ b/rasterio/rio/main.py
@@ -11,6 +11,7 @@ import click
 from click_plugins import with_plugins
 import cligj
 
+from . import options
 import rasterio
 
 
diff --git a/rasterio/rio/merge.py b/rasterio/rio/merge.py
index 5e637c8..6e37193 100644
--- a/rasterio/rio/merge.py
+++ b/rasterio/rio/merge.py
@@ -3,7 +3,6 @@
 import logging
 import math
 import os.path
-import warnings
 
 import click
 from cligj import files_inout_arg, format_opt
diff --git a/rasterio/rio/options.py b/rasterio/rio/options.py
index dfccd54..b5da2ae 100644
--- a/rasterio/rio/options.py
+++ b/rasterio/rio/options.py
@@ -40,14 +40,18 @@ Registry of common rio CLI options.  See cligj for more options.
 --width: width of raster.  In rio-info
 --with-nodata/--without-nodata: include nodata regions or not.  In rio-shapes.
 -v, --tell-me-more, --verbose
+--vfs: virtual file system.
 """
 
 
 # TODO: move file_in_arg and file_out_arg to cligj
 
+import os.path
 
 import click
 
+from rasterio.vfs import parse_path
+
 
 def _cb_key_val(ctx, param, value):
 
@@ -76,14 +80,29 @@ def _cb_key_val(ctx, param, value):
                 k = k.lower()
                 v = v.lower()
                 out[k] = v
-
         return out
 
 
+def file_in_handler(ctx, param, value):
+    """Normalize ordinary filesystem and VFS paths"""
+    try:
+        path, archive, scheme = parse_path(value)
+    except ValueError as exc:
+        raise click.BadParameter(str(exc))
+    path_to_check = archive or path
+    if not os.path.exists(path_to_check):
+        raise click.BadParameter(
+            "Input file {0} does not exist".format(path_to_check))
+    if archive and scheme:
+        archive = os.path.abspath(archive)
+        path = "{0}://{1}!{2}".format(scheme, archive, path)
+    else:
+        path = os.path.abspath(path)
+    return path
+
+
 # Singular input file
-file_in_arg = click.argument(
-    'INPUT',
-    type=click.Path(exists=True, resolve_path=True))
+file_in_arg = click.argument('INPUT', callback=file_in_handler)
 
 # Singular output file
 file_out_arg = click.argument(
@@ -161,4 +180,4 @@ rgb_opt = click.option(
     '--rgb', 'photometric', 
     flag_value='rgb',
     default=False,
-    help="Set RGB photometric interpretation")
+    help="Set RGB photometric interpretation.")
diff --git a/rasterio/rio/sample.py b/rasterio/rio/sample.py
index 8054bec..16cc0c6 100644
--- a/rasterio/rio/sample.py
+++ b/rasterio/rio/sample.py
@@ -1,15 +1,11 @@
 import json
 import logging
-import warnings
 
 import click
 
 import rasterio
 
 
-warnings.simplefilter('default')
-
-
 @click.command(short_help="Sample a dataset.")
 @click.argument('files', nargs=-1, required=True, metavar='FILE "[x, y]"')
 @click.option('-b', '--bidx', default=None, help="Indexes of input file bands.")
diff --git a/rasterio/tools/merge.py b/rasterio/tools/merge.py
index 0fea083..b073183 100644
--- a/rasterio/tools/merge.py
+++ b/rasterio/tools/merge.py
@@ -1,3 +1,5 @@
+from __future__ import absolute_import
+
 import logging
 import math
 import warnings
diff --git a/rasterio/transform.py b/rasterio/transform.py
index b832711..3bb073f 100644
--- a/rasterio/transform.py
+++ b/rasterio/transform.py
@@ -1,3 +1,5 @@
+from __future__ import absolute_import
+
 import warnings
 
 from affine import Affine
diff --git a/rasterio/vfs.py b/rasterio/vfs.py
new file mode 100644
index 0000000..4af1faf
--- /dev/null
+++ b/rasterio/vfs.py
@@ -0,0 +1,43 @@
+"""Implementation of Apache VFS schemes and URLs"""
+
+import os
+
+
+# NB: As not to propagate fallacies of distributed computing, Rasterio
+# does not support HTTP or FTP URLs via GDAL's vsicurl handler. Only
+# the following local filesystem schemes are supported.
+SCHEMES = ['gzip', 'zip', 'tar']
+
+
+def parse_path(path, vfs=None):
+    """Parse a file path or Apache VFS URL into its parts."""
+    archive = scheme = None
+    if vfs:
+        parts = vfs.split("://")
+        scheme = parts.pop(0) if parts else None
+        archive = parts.pop(0) if parts else None
+    else:
+        parts = path.split("://")
+        path = parts.pop() if parts else None
+        scheme = parts.pop() if parts else None
+        if scheme in SCHEMES:
+            parts = path.split('!')
+            path = parts.pop() if parts else None
+            archive = parts.pop() if parts else None
+        elif scheme in (None, 'file'):
+            pass
+        else:
+            raise ValueError("VFS scheme {0} is unknown".format(scheme))
+    return path, archive, scheme
+
+
+def vsi_path(path, archive=None, scheme=None):
+    """Convert a parsed path to a GDAL VSI path."""
+    # If a VSF and archive file are specified, we convert the path to
+    # a GDAL VSI path (see cpl_vsi.h).
+    if scheme and scheme != 'file':
+        path = path.strip(os.path.sep)
+        result = os.path.sep.join(['/vsi{0}'.format(scheme), archive, path])
+    else:
+        result = path
+    return result
diff --git a/rasterio/warnings.py b/rasterio/warnings.py
new file mode 100644
index 0000000..bdc22b4
--- /dev/null
+++ b/rasterio/warnings.py
@@ -0,0 +1,9 @@
+"""Rasterio warnings."""
+
+
+class NodataShadowWarning(Warning):
+    """Warn that a dataset's nodata attribute is shadowing its alpha band"""
+    def __str__(self):
+        return ("The dataset's nodata attribute is shadowing "
+                "the alpha band. All masks will be determined "
+                "by the nodata attribute")
diff --git a/requirements-dev.txt b/requirements-dev.txt
index b9e2dd1..16048fe 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,12 +1,11 @@
 affine
 cligj
-coveralls>=0.4
-cython>=0.23.1
+cython>=0.23.4
 delocate
 enum34
 numpy>=1.8
 snuggs>=1.2
-pytest
-pytest-cov
+pytest>=2.8.2
+pytest-cov>=2.2.0
 setuptools>=0.9.8
 wheel
diff --git a/setup.cfg b/setup.cfg
index 5eea52d..5ee6477 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,8 +1,2 @@
-[nosetests]
-tests=rasterio/tests
-nocapture=True
-verbosity=3
-logging-filter=rasterio
-logging-level=DEBUG
-with-coverage=1
-cover-package=rasterio
+[pytest]
+testpaths = tests
diff --git a/setup.py b/setup.py
index 0cefec9..60789fd 100755
--- a/setup.py
+++ b/setup.py
@@ -19,12 +19,10 @@ import sys
 from setuptools import setup
 from setuptools.extension import Extension
 
+
 logging.basicConfig()
 log = logging.getLogger()
 
-# python -W all setup.py ...
-if 'all' in sys.warnoptions:
-    log.level = logging.DEBUG
 
 def check_output(cmd):
     # since subprocess.check_output doesn't exist in 2.6
@@ -39,6 +37,7 @@ def check_output(cmd):
         out, err = p.communicate()
         return out
 
+
 def copy_data_tree(datadir, destdir):
     try:
         shutil.rmtree(destdir)
@@ -46,6 +45,11 @@ def copy_data_tree(datadir, destdir):
         pass
     shutil.copytree(datadir, destdir)
 
+
+# python -W all setup.py ...
+if 'all' in sys.warnoptions:
+    log.level = logging.DEBUG
+
 # Parse the version from the rasterio module.
 with open('rasterio/__init__.py') as f:
     for line in f:
@@ -135,6 +139,13 @@ if not os.name == "nt":
     ext_options['extra_compile_args'] = ['-Wno-unused-parameter',
                                          '-Wno-unused-function']
 
+cythonize_options = {}
+if os.environ.get('CYTHON_COVERAGE'):
+    cythonize_options['compiler_directives'] = {'linetrace': True}
+    cythonize_options['annotate'] = True
+    ext_options['define_macros'] = [('CYTHON_TRACE', '1'),
+                                    ('CYTHON_TRACE_NOGIL', '1')]
+
 log.debug('ext_options:\n%s', pprint.pformat(ext_options))
 
 # When building from a repo, Cython is required.
@@ -164,7 +175,7 @@ if os.path.exists("MANIFEST.in") and "clean" not in sys.argv:
             'rasterio._err', ['rasterio/_err.pyx'], **ext_options),
         Extension(
             'rasterio._example', ['rasterio/_example.pyx'], **ext_options),
-        ], quiet=True)
+        ], quiet=True, **cythonize_options)
 
 # If there's no manifest template, as in an sdist, we just specify .c files.
 else:
diff --git a/tests/conftest.py b/tests/conftest.py
index 90aed75..16f39e7 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -7,6 +7,10 @@ import sys
 from click.testing import CliRunner
 import py
 import pytest
+import numpy
+
+
+DEFAULT_SHAPE = (10, 10)
 
 
 if sys.version_info > (3,):
@@ -38,3 +42,163 @@ def data():
     for filename in test_files:
         shutil.copy(filename, str(tmpdir))
     return tmpdir
+
+
+ at pytest.fixture
+def basic_geometry():
+    """
+    Returns
+    -------
+
+    dict: GeoJSON-style geometry object.
+        Coordinates are in grid coordinates (Affine.identity()).
+    """
+
+    return {
+        'type': 'Polygon',
+        'coordinates': [[(2, 2), (2, 4.25), (4.25, 4.25), (4.25, 2), (2, 2)]]
+    }
+
+
+ at pytest.fixture
+def basic_feature(basic_geometry):
+    """
+    Returns
+    -------
+
+    dict: GeoJSON object.
+        Coordinates are in grid coordinates (Affine.identity()).
+    """
+
+    return {
+        'geometry': basic_geometry,
+        'properties': {
+            'val': 15
+        },
+        'type': 'Feature'
+    }
+
+
+ at pytest.fixture
+def basic_featurecollection(basic_feature):
+    """
+    Returns
+    -------
+
+    dict: GeoJSON FeatureCollection object.
+        Coordinates are in grid coordinates (Affine.identity()).
+    """
+
+    return {
+        'features': [basic_feature],
+        'type': 'FeatureCollection'
+    }
+
+
+ at pytest.fixture
+def basic_image():
+    """
+    A basic 10x10 array for testing sieve and shapes functions.
+    Contains a square feature 3x3 (size 9).
+    Equivalent to results of rasterizing basic_geometry with all_touched=True.
+
+    Returns
+    -------
+
+    numpy ndarray
+    """
+
+    image = numpy.zeros(DEFAULT_SHAPE, dtype=numpy.uint8)
+    image[2:5, 2:5] = 1
+
+    return image
+
+
+ at pytest.fixture
+def basic_image_2x2():
+    """
+    A basic 10x10 array for testing sieve and shapes functions.
+    Contains a square feature 2x2 (size 4).
+    Equivalent to results of rasterizing basic_geometry with all_touched=False.
+
+    Returns
+    -------
+
+    numpy ndarray
+    """
+
+    image = numpy.zeros(DEFAULT_SHAPE, dtype=numpy.uint8)
+    image[2:4, 2:4] = 1
+
+    return image
+
+
+ at pytest.fixture
+def pixelated_image(basic_image):
+    """
+    A basic 10x10 array for testing sieve functions.  Contains a square feature
+    3x3 (size 9), with 2 isolated pixels.
+
+    Returns
+    -------
+
+    numpy ndarray
+    """
+
+    image = basic_image.copy()
+    image [0, 0] = 1
+    image [8, 8] = 1
+
+    return image
+
+
+ at pytest.fixture
+def diagonal_image():
+    """
+    A 10x10 array for testing sieve functions, with only one diagonal filled.
+
+    Returns
+    -------
+
+    numpy ndarray
+    """
+
+    image = numpy.zeros(DEFAULT_SHAPE, dtype=numpy.uint8)
+    numpy.fill_diagonal(image, 1)
+    return image
+
+
+ at pytest.fixture()
+def pixelated_image_file(tmpdir, pixelated_image):
+    """
+    A basic raster file with a 10x10 array for testing sieve functions.
+    Contains data from pixelated_image.
+
+    Returns
+    -------
+
+    string
+        Filename of test raster file
+    """
+
+    from affine import Affine
+    import rasterio
+
+    image = pixelated_image
+
+    outfilename = str(tmpdir.join('pixelated_image.tif'))
+    kwargs = {
+        "crs": {'init': 'epsg:4326'},
+        "transform": Affine.identity(),
+        "count": 1,
+        "dtype": rasterio.uint8,
+        "driver": "GTiff",
+        "width": image.shape[1],
+        "height": image.shape[0],
+        "nodata": 255
+    }
+    with rasterio.drivers():
+        with rasterio.open(outfilename, 'w', **kwargs) as out:
+            out.write_band(1, image)
+
+    return outfilename
diff --git a/tests/data/files.zip b/tests/data/files.zip
new file mode 100644
index 0000000..64a8474
Binary files /dev/null and b/tests/data/files.zip differ
diff --git a/tests/test_band_masks.py b/tests/test_band_masks.py
index 36dd382..0624b10 100644
--- a/tests/test_band_masks.py
+++ b/tests/test_band_masks.py
@@ -2,14 +2,91 @@ import logging
 import sys
 
 import numpy as np
-from pytest import fixture
+import pytest
 
 import rasterio
+from rasterio.enums import MaskFlags
+from rasterio.warnings import NodataShadowWarning
 
 
 logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
 
 
+
+
+ at pytest.fixture(scope='function')
+def tiffs(tmpdir):
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        profile = src.profile
+
+        shadowed_profile = profile.copy()
+        shadowed_profile['count'] = 4
+        with rasterio.open(
+                str(tmpdir.join('shadowed.tif')), 'w',
+                **shadowed_profile) as dst:
+
+            for i, band in enumerate(src.read(masked=False), 1):
+                dst.write(band, i)
+            dst.write(band, 4)
+
+        del profile['nodata']
+        with rasterio.open(
+                str(tmpdir.join('no-nodata.tif')), 'w',
+                **profile) as dst:
+            dst.write(src.read(masked=False))
+
+        with rasterio.open(
+                str(tmpdir.join('sidecar-masked.tif')), 'w',
+                **profile) as dst:
+            dst.write(src.read(masked=False))
+            mask = np.zeros(src.shape, dtype='uint8')
+            dst.write_mask(mask)
+
+    return tmpdir
+
+def test_mask_flags():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        for flags in src.mask_flags:
+            assert flags & MaskFlags.nodata
+            assert not flags & MaskFlags.per_dataset
+            assert not flags & MaskFlags.alpha
+
+
+def test_mask_flags_sidecar(tiffs):
+    filename = str(tiffs.join('sidecar-masked.tif'))
+    with rasterio.open(filename) as src:
+        for flags in src.mask_flags:
+            assert not flags & MaskFlags.nodata
+            assert not flags & MaskFlags.alpha
+            assert flags & MaskFlags.per_dataset
+
+
+def test_mask_flags_shadow(tiffs):
+    filename = str(tiffs.join('shadowed.tif'))
+    with rasterio.open(filename) as src:
+        for flags in src.mask_flags:
+            assert flags & MaskFlags.nodata
+            assert not flags & MaskFlags.alpha
+            assert not flags & MaskFlags.per_dataset
+
+
+def test_warning_no():
+    """No shadow warning is raised"""
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        try:
+            rm, gm, bm = src.read_masks()
+        except NodataShadowWarning:
+            pytest.fail("Unexpected NodataShadowWarning raised")
+
+
+def test_warning_shadow(tiffs):
+    """Shadow warning is raised"""
+    filename = str(tiffs.join('shadowed.tif'))
+    with rasterio.open(filename) as src:
+        with pytest.warns(NodataShadowWarning):
+            _ = src.read_masks()
+
+
 def test_masks():
     with rasterio.open('tests/data/RGB.byte.tif') as src:
         rm, gm, bm = src.read_masks()
@@ -37,27 +114,6 @@ def test_masked_none():
         assert (b.mask==~bm.astype('bool')).all()
 
 
- at fixture(scope='function')
-def tiffs(tmpdir):
-    with rasterio.open('tests/data/RGB.byte.tif') as src:
-        kwds = src.meta
-        
-        del kwds['nodata']
-        with rasterio.open(
-                str(tmpdir.join('no-nodata.tif')), 'w',
-                **kwds) as dst:
-            dst.write(src.read(masked=False))
-
-        with rasterio.open(
-                str(tmpdir.join('sidecar-masked.tif')), 'w',
-                **kwds) as dst:
-            dst.write(src.read(masked=False))
-            mask = np.zeros(src.shape, dtype='uint8')
-            dst.write_mask(mask)
-
-    return tmpdir
-
-
 def test_masking_no_nodata(tiffs):
     # if the dataset has no defined nodata values, all data is
     # considered valid data. The GDAL masks bands are arrays of
@@ -65,6 +121,10 @@ def test_masking_no_nodata(tiffs):
     # is False.
     filename = str(tiffs.join('no-nodata.tif'))
     with rasterio.open(filename) as src:
+        for flags in src.mask_flags:
+            assert flags & MaskFlags.all_valid
+            assert not flags & MaskFlags.alpha
+            assert not flags & MaskFlags.nodata
 
         rgb = src.read(masked=False)
         assert not hasattr(rgb, 'mask')
@@ -92,6 +152,10 @@ def test_masking_sidecar_mask(tiffs):
     # If the dataset has a .msk sidecar mask band file, all masks will
     # be derived from that file.
     with rasterio.open(str(tiffs.join('sidecar-masked.tif'))) as src:
+        for flags in src.mask_flags:
+            assert flags & MaskFlags.per_dataset
+            assert not flags & MaskFlags.alpha
+            assert not flags & MaskFlags.nodata
         rgb = src.read(masked=True)
         assert rgb.mask.all()
         r = src.read(1, masked=True)
diff --git a/tests/test_blocks.py b/tests/test_blocks.py
index 912baf5..4ca99be 100644
--- a/tests/test_blocks.py
+++ b/tests/test_blocks.py
@@ -1,15 +1,16 @@
 import logging
 import os.path
-import unittest
 import shutil
 import subprocess
 import sys
 import tempfile
+import unittest
 
 import numpy
 
 import rasterio
 
+
 logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
 
 class WindowTest(unittest.TestCase):
@@ -120,4 +121,3 @@ class WindowWriteTest(unittest.TestCase):
             "Minimum=0.000, Maximum=127.000, "
             "Mean=31.750, StdDev=54.993" in info.decode('utf-8'),
             info)
-
diff --git a/tests/test_colormap.py b/tests/test_colormap.py
index 5fc7778..2008fa0 100644
--- a/tests/test_colormap.py
+++ b/tests/test_colormap.py
@@ -1,8 +1,6 @@
 import logging
-import pytest
 import subprocess
 import sys
-import warnings
 
 import rasterio
 
diff --git a/tests/test_dtypes.py b/tests/test_dtypes.py
index 1826ec5..9c0310b 100644
--- a/tests/test_dtypes.py
+++ b/tests/test_dtypes.py
@@ -1,23 +1,58 @@
 import numpy as np
 
-from rasterio import dtypes, ubyte
+from rasterio import (
+    ubyte, uint8, uint16, uint32, int16, int32, float32, float64)
+from rasterio.dtypes import (
+    _gdal_typename, is_ndarray, check_dtype, get_minimum_dtype, can_cast_dtype,
+    validate_dtype
+)
 
 
 def test_is_ndarray():
-    assert dtypes.is_ndarray(np.zeros((1,)))
-    assert dtypes.is_ndarray([0]) == False
-    assert dtypes.is_ndarray((0,)) == False
+    assert is_ndarray(np.zeros((1,)))
+    assert is_ndarray([0]) == False
+    assert is_ndarray((0,)) == False
 
 
 def test_np_dt_uint8():
-    assert dtypes.check_dtype(np.uint8)
+    assert check_dtype(np.uint8)
 
 
 def test_dt_ubyte():
-    assert dtypes.check_dtype(ubyte)
+    assert check_dtype(ubyte)
+
+
+def test_check_dtype_invalid():
+    assert check_dtype('foo') == False
 
 
 def test_gdal_name():
-    assert dtypes._gdal_typename(ubyte) == 'Byte'
-    assert dtypes._gdal_typename(np.uint8) == 'Byte'
-    assert dtypes._gdal_typename(np.uint16) == 'UInt16'
+    assert _gdal_typename(ubyte) == 'Byte'
+    assert _gdal_typename(np.uint8) == 'Byte'
+    assert _gdal_typename(np.uint16) == 'UInt16'
+
+
+def test_get_minimum_dtype():
+    assert get_minimum_dtype([0, 1]) == uint8
+    assert get_minimum_dtype([0, 1000]) == uint16
+    assert get_minimum_dtype([0, 100000]) == uint32
+    assert get_minimum_dtype([-1, 0, 1]) == int16
+    assert get_minimum_dtype([-1, 0, 100000]) == int32
+    assert get_minimum_dtype([-1.5, 0, 1.5]) == float32
+    assert get_minimum_dtype([-1.5e+100, 0, 1.5e+100]) == float64
+
+
+def test_can_cast_dtype():
+    assert can_cast_dtype((1, 2, 3), np.uint8) == True
+    assert can_cast_dtype(np.array([1, 2, 3]), np.uint8) == True
+    assert can_cast_dtype(np.array([1, 2, 3], dtype=np.uint8), np.uint8) == True
+    assert can_cast_dtype(np.array([1, 2, 3]), np.float32) == True
+    assert can_cast_dtype(np.array([1.4, 2.1, 3.65]), np.float32) == True
+    assert can_cast_dtype(np.array([1.4, 2.1, 3.65]), np.uint8) == False
+
+
+def test_validate_dtype():
+    assert validate_dtype([1, 2, 3], ('uint8', 'uint16')) == True
+    assert validate_dtype(np.array([1, 2, 3]), ('uint8', 'uint16')) == True
+    assert validate_dtype(np.array([1.4, 2.1, 3.65]), ('float32',)) == True
+    assert validate_dtype(np.array([1.4, 2.1, 3.65]),('uint8',)) == False
diff --git a/tests/test_features.py b/tests/test_features.py
new file mode 100644
index 0000000..bdbdab2
--- /dev/null
+++ b/tests/test_features.py
@@ -0,0 +1,728 @@
+import logging
+import sys
+import numpy
+import pytest
+
+from affine import Affine
+import rasterio
+from rasterio.features import bounds, geometry_mask, rasterize, sieve, shapes
+
+
+DEFAULT_SHAPE = (10, 10)
+
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+
+def test_bounds_point():
+    g = {'type': 'Point', 'coordinates': [10, 10]}
+    assert bounds(g) == (10, 10, 10, 10)
+
+
+def test_bounds_line():
+    g = {'type': 'LineString', 'coordinates': [[0, 0], [10, 10]]}
+    assert bounds(g) == (0, 0, 10, 10)
+
+
+def test_bounds_polygon():
+    g = {'type': 'Polygon', 'coordinates': [[[0, 0], [10, 10], [10, 0]]]}
+    assert bounds(g) == (0, 0, 10, 10)
+
+
+def test_bounds_z():
+    g = {'type': 'Point', 'coordinates': [10, 10, 10]}
+    assert bounds(g) == (10, 10, 10, 10)
+
+
+def test_bounds_invalid_obj():
+    with pytest.raises(KeyError):
+        bounds({'type': 'bogus', 'not_coordinates': []})
+
+
+def test_feature_collection(basic_featurecollection):
+    fc = basic_featurecollection
+    assert bounds(fc) == bounds(fc['features'][0]) == (2, 2, 4.25, 4.25)
+
+
+def test_bounds_existing_bbox(basic_featurecollection):
+    """
+    Test with existing bbox in geojson, similar to that produced by
+    rasterio.  Values specifically modified here for testing, bboxes are not
+    valid as written.
+    """
+
+    fc = basic_featurecollection
+    fc['bbox'] = [0, 10, 10, 20]
+    fc['features'][0]['bbox'] = [0, 100, 10, 200]
+
+    assert bounds(fc['features'][0]) == (0, 100, 10, 200)
+    assert bounds(fc) == (0, 10, 10, 20)
+
+
+def test_geometry_mask(basic_geometry, basic_image_2x2):
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            basic_image_2x2 == 0,
+            geometry_mask(
+                [basic_geometry],
+                out_shape=DEFAULT_SHAPE,
+                transform=Affine.identity()
+            )
+        )
+
+
+def test_geometry_mask_invert(basic_geometry, basic_image_2x2):
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            basic_image_2x2,
+            geometry_mask(
+                [basic_geometry],
+                out_shape=DEFAULT_SHAPE,
+                transform=Affine.identity(),
+                invert=True
+            )
+        )
+
+
+def test_rasterize(basic_geometry, basic_image_2x2):
+    """ Rasterize operation should succeed for both an out_shape and out """
+
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            basic_image_2x2,
+            rasterize([basic_geometry], out_shape=DEFAULT_SHAPE)
+        )
+
+        out = numpy.zeros(DEFAULT_SHAPE)
+        rasterize([basic_geometry], out=out)
+        assert numpy.array_equal(basic_image_2x2, out)
+
+
+def test_rasterize_invalid_out_dtype(basic_geometry):
+    """ A non-supported data type for out should raise an exception """
+
+    out = numpy.zeros(DEFAULT_SHAPE, dtype=numpy.int64)
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            rasterize([basic_geometry], out=out)
+
+
+def test_rasterize_shapes_out_dtype_mismatch(basic_geometry):
+    """ Shape values must be able to fit in data type for out """
+
+    out = numpy.zeros(DEFAULT_SHAPE, dtype=numpy.uint8)
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            rasterize([(basic_geometry, 10000000)], out=out)
+
+
+def test_rasterize_missing_out(basic_geometry):
+    """ If both out and out_shape are missing, should raise exception """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            rasterize([basic_geometry], out=None, out_shape=None)
+
+
+def test_rasterize_missing_shapes():
+    """ Shapes are required for this operation """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError) as ex:
+            rasterize([], out_shape=DEFAULT_SHAPE)
+
+        assert 'No valid geometry objects' in str(ex.value)
+
+
+def test_rasterize_invalid_shapes():
+    """ Invalid shapes should raise an exception rather than be skipped """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError) as ex:
+            rasterize([{'foo': 'bar'}], out_shape=DEFAULT_SHAPE)
+
+        assert 'Invalid geometry object' in str(ex.value)
+
+
+def test_rasterize_default_value(basic_geometry, basic_image_2x2):
+    """ All shapes should rasterize to the default value """
+
+    default_value = 2
+    truth = basic_image_2x2 * default_value
+
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            truth,
+            rasterize(
+                [basic_geometry], out_shape=DEFAULT_SHAPE,
+                default_value=default_value
+            )
+        )
+
+
+def test_rasterize_invalid_default_value(basic_geometry):
+    """ A default value that requires an int64 should raise an exception """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            rasterize(
+                [basic_geometry], out_shape=DEFAULT_SHAPE,
+                default_value=1000000000000
+            )
+
+
+def test_rasterize_fill_value(basic_geometry, basic_image_2x2):
+    """ All pixels not covered by shapes should be given fill value """
+
+    default_value = 2
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            basic_image_2x2 + 1,
+            rasterize(
+                [basic_geometry], out_shape=DEFAULT_SHAPE, fill=1,
+                default_value=default_value
+            )
+        )
+
+
+def test_rasterize_invalid_fill_value(basic_geometry):
+    """ A fill value that requires an int64 should raise an exception """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            rasterize(
+                [basic_geometry], out_shape=DEFAULT_SHAPE, fill=1000000000000,
+                default_value=2
+            )
+
+
+def test_rasterize_fill_value_dtype_mismatch(basic_geometry):
+    """ A fill value that doesn't match dtype should fail """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            rasterize(
+                [basic_geometry], out_shape=DEFAULT_SHAPE, fill=1000000,
+                default_value=2, dtype=numpy.uint8
+            )
+
+
+def test_rasterize_all_touched(basic_geometry, basic_image):
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            basic_image,
+            rasterize(
+                [basic_geometry], out_shape=DEFAULT_SHAPE, all_touched=True
+            )
+        )
+
+
+def test_rasterize_value(basic_geometry, basic_image_2x2):
+    """
+    All shapes should rasterize to the value passed in a tuple alongside
+    each shape
+    """
+
+    value = 5
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            basic_image_2x2 * value,
+            rasterize(
+                [(basic_geometry, value)], out_shape=DEFAULT_SHAPE
+            )
+        )
+
+
+def test_rasterize_invalid_value(basic_geometry):
+    """ A shape value that requires an int64 should raise an exception """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError) as ex:
+            rasterize(
+                [(basic_geometry, 1000000000000)], out_shape=DEFAULT_SHAPE
+            )
+
+        assert 'dtype must be one of' in str(ex.value)
+
+
+def test_rasterize_supported_dtype(basic_geometry):
+    """ Supported data types should return valid results """
+
+    with rasterio.drivers():
+        supported_types = (
+            ('int16', -32768),
+            ('int32', -2147483648),
+            ('uint8', 255),
+            ('uint16', 65535),
+            ('uint32', 4294967295),
+            ('float32', 1.434532),
+            ('float64', -98332.133422114)
+        )
+
+        for dtype, default_value in supported_types:
+            truth = numpy.zeros(DEFAULT_SHAPE, dtype=dtype)
+            truth[2:4, 2:4] = default_value
+
+            result = rasterize(
+                [basic_geometry],
+                out_shape=DEFAULT_SHAPE,
+                default_value=default_value,
+                dtype=dtype
+            )
+            assert numpy.array_equal(result, truth)
+            assert numpy.dtype(result.dtype) == numpy.dtype(truth.dtype)
+
+            result = rasterize(
+                [(basic_geometry, default_value)],
+                out_shape=DEFAULT_SHAPE
+            )
+            if numpy.dtype(dtype).kind == 'f':
+                assert numpy.allclose(result, truth)
+            else:
+                assert numpy.array_equal(result, truth)
+            # Since dtype is auto-detected, it may not match due to upcasting
+
+
+def test_rasterize_unsupported_dtype(basic_geometry):
+    """ Unsupported types should all raise exceptions """
+
+    with rasterio.drivers():
+        unsupported_types = (
+            ('int8', -127),
+            ('int64', 20439845334323),
+            ('float16', -9343.232)
+        )
+
+        for dtype, default_value in unsupported_types:
+            with pytest.raises(ValueError):
+                rasterize(
+                    [basic_geometry],
+                    out_shape=DEFAULT_SHAPE,
+                    default_value=default_value,
+                    dtype=dtype
+                )
+
+            with pytest.raises(ValueError):
+                rasterize(
+                    [(basic_geometry, default_value)],
+                    out_shape=DEFAULT_SHAPE,
+                    dtype=dtype
+                )
+
+
+def test_rasterize_mismatched_dtype(basic_geometry):
+    """ Mismatched values and dtypes should raise exceptions """
+
+    with rasterio.drivers():
+        mismatched_types = (('uint8', 3.2423), ('uint8', -2147483648))
+        for dtype, default_value in mismatched_types:
+            with pytest.raises(ValueError):
+                rasterize(
+                    [basic_geometry],
+                    out_shape=DEFAULT_SHAPE,
+                    default_value=default_value,
+                    dtype=dtype
+                )
+
+            with pytest.raises(ValueError):
+                rasterize(
+                    [(basic_geometry, default_value)],
+                    out_shape=DEFAULT_SHAPE,
+                    dtype=dtype
+                )
+
+
+def test_rasterize_geometries_symmetric():
+    """ Make sure that rasterize is symmetric with shapes """
+
+    transform = (1.0, 0.0, 0.0, 0.0, -1.0, 0.0)
+    truth = numpy.zeros(DEFAULT_SHAPE, dtype=rasterio.ubyte)
+    truth[2:5, 2:5] = 1
+    with rasterio.drivers():
+        s = shapes(truth, transform=transform)
+        result = rasterize(s, out_shape=DEFAULT_SHAPE, transform=transform)
+        assert numpy.array_equal(result, truth)
+
+
+def test_rasterize_internal_driver_manager(basic_geometry):
+    """ Rasterize should work without explicitly calling driver manager """
+
+    assert rasterize([basic_geometry], out_shape=DEFAULT_SHAPE).sum() == 4
+
+
+def test_shapes(basic_image):
+    """ Test creation of shapes from pixel values """
+
+    with rasterio.drivers():
+        results = list(shapes(basic_image))
+
+        assert len(results) == 2
+
+        shape, value = results[0]
+        assert shape == {
+            'coordinates': [
+                [(2, 2), (2, 5), (5, 5), (5, 2), (2, 2)]
+            ],
+            'type': 'Polygon'
+        }
+        assert value == 1
+
+        shape, value = results[1]
+        assert shape == {
+            'coordinates': [
+                [(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)],
+                [(2, 2), (5, 2), (5, 5), (2, 5), (2, 2)]
+            ],
+            'type': 'Polygon'
+        }
+        assert value == 0
+
+
+def test_shapes_band(pixelated_image, pixelated_image_file):
+    """ Shapes from a band should match shapes from an array """
+
+    with rasterio.drivers():
+        truth = list(shapes(pixelated_image))
+
+        with rasterio.open(pixelated_image_file) as src:
+            band = rasterio.band(src, 1)
+            assert truth == list(shapes(band))
+
+            # Mask band should function, but will mask out some results
+            assert truth[0] == list(shapes(band, mask=band))[0]
+
+
+def test_shapes_connectivity_rook(diagonal_image):
+    """
+    Diagonals are not connected, so there will be 1 feature per pixel plus
+    background.
+    """
+
+    with rasterio.drivers():
+        assert len(list(shapes(diagonal_image, connectivity=4))) == 12
+
+
+def test_shapes_connectivity_queen(diagonal_image):
+    """
+    Diagonals are connected, so there will be 1 feature for all pixels plus
+    background.
+    """
+
+    with rasterio.drivers():
+        assert len(list(shapes(diagonal_image, connectivity=8))) == 2
+
+
+def test_shapes_connectivity_invalid(diagonal_image):
+    """ Invalid connectivity should raise exception """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            assert next(shapes(diagonal_image, connectivity=12))
+
+
+def test_shapes_mask(basic_image):
+    """ Only pixels not masked out should be converted to features """
+
+    mask = numpy.ones(basic_image.shape, dtype=rasterio.bool_)
+    mask[4:5, 4:5] = False
+
+    with rasterio.drivers():
+        results = list(shapes(basic_image, mask=mask))
+
+        assert len(results) == 2
+
+        shape, value = results[0]
+        assert shape == {
+            'coordinates': [
+                [(2, 2), (2, 5), (4, 5), (4, 4), (5, 4), (5, 2), (2, 2)]
+            ],
+            'type': 'Polygon'
+        }
+        assert value == 1
+
+
+def test_shapes_blank_mask(basic_image):
+    """ Mask is blank so results should mask shapes without mask """
+
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            list(shapes(
+                basic_image,
+                mask=numpy.ones(basic_image.shape, dtype=rasterio.bool_))
+            ),
+            list(shapes(basic_image))
+        )
+
+
+def test_shapes_invalid_mask_shape(basic_image):
+    """ A mask that is the wrong shape should fail """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            next(shapes(
+                basic_image,
+                mask=numpy.ones(
+                    (basic_image.shape[0] + 10, basic_image.shape[1] + 10),
+                    dtype=rasterio.bool_
+                )
+            ))
+
+
+def test_shapes_invalid_mask_dtype(basic_image):
+    """ A mask that is the wrong dtype should fail """
+
+    with rasterio.drivers():
+        for dtype in ('int8', 'int16', 'int32'):
+            with pytest.raises(ValueError):
+                next(shapes(
+                    basic_image,
+                    mask=numpy.ones(basic_image.shape, dtype=dtype)
+                ))
+
+
+def test_shapes_supported_dtypes(basic_image):
+    """ Supported data types should return valid results """
+
+    supported_types = (
+        ('int16', -32768),
+        ('int32', -2147483648),
+        ('uint8', 255),
+        ('uint16', 65535),
+        ('float32', 1.434532)
+    )
+
+    with rasterio.drivers():
+        for dtype, test_value in supported_types:
+            shape, value = next(shapes(basic_image.astype(dtype) * test_value))
+            assert numpy.allclose(value, test_value)
+
+
+def test_shapes_unsupported_dtypes(basic_image):
+    """ Unsupported data types should raise exceptions """
+
+    unsupported_types = (
+        ('int8', -127),
+        ('uint32', 4294967295),
+        ('int64', 20439845334323),
+        ('float16', -9343.232),
+        ('float64', -98332.133422114)
+    )
+
+    with rasterio.drivers():
+        for dtype, test_value in unsupported_types:
+            with pytest.raises(ValueError):
+                next(shapes(basic_image.astype(dtype) * test_value))
+
+
+def test_shapes_internal_driver_manager(basic_image):
+    """ Shapes should work without explicitly calling driver manager """
+
+    assert next(shapes(basic_image))[0]['type'] == 'Polygon'
+
+
+def test_sieve_small(basic_image, pixelated_image):
+    """
+    Setting the size smaller than or equal to the size of the feature in the
+    image should not change the image.
+    """
+
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            basic_image,
+            sieve(pixelated_image, basic_image.sum())
+        )
+
+
+def test_sieve_large(basic_image):
+    """
+    Setting the size larger than size of feature should leave us an empty image.
+    """
+
+    with rasterio.drivers():
+        assert not numpy.any(sieve(basic_image, basic_image.sum() + 1))
+
+
+def test_sieve_invalid_size(basic_image):
+    with rasterio.drivers():
+        for invalid_size in (0, 45.1234, basic_image.size + 1):
+            with pytest.raises(ValueError):
+                sieve(basic_image, invalid_size)
+
+
+def test_sieve_connectivity_rook(diagonal_image):
+    """ Diagonals are not connected, so feature is removed """
+
+    assert not numpy.any(
+        sieve(diagonal_image, diagonal_image.sum(), connectivity=4)
+    )
+
+
+def test_sieve_connectivity_queen(diagonal_image):
+    """ Diagonals are connected, so feature is retained """
+
+    assert numpy.array_equal(
+        diagonal_image,
+        sieve(diagonal_image, diagonal_image.sum(), connectivity=8)
+    )
+
+
+def test_sieve_connectivity_invalid(basic_image):
+    with pytest.raises(ValueError):
+        sieve(basic_image, 54, connectivity=12)
+
+
+def test_sieve_out(basic_image):
+    """ Output array passed in should match the returned array """
+
+    with rasterio.drivers():
+        output = numpy.zeros_like(basic_image)
+        output[1:3, 1:3] = 5
+        sieved_image = sieve(basic_image, basic_image.sum(), out=output)
+        assert numpy.array_equal(basic_image, sieved_image)
+        assert numpy.array_equal(output, sieved_image)
+
+
+def test_sieve_invalid_out(basic_image):
+    """ Output with different dtype or shape should fail """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            sieve(
+                basic_image, basic_image.sum(),
+                out=numpy.zeros(basic_image.shape, dtype=rasterio.int32)
+            )
+
+        with pytest.raises(ValueError):
+            sieve(
+                basic_image, basic_image.sum(),
+                out=numpy.zeros(
+                    (basic_image.shape[0] + 10, basic_image.shape[1] + 10),
+                    dtype=rasterio.ubyte
+                )
+            )
+
+
+def test_sieve_mask(basic_image):
+    """
+    Only areas within the overlap of mask and input will be kept, so long
+    as mask is a bool or uint8 dtype.
+    """
+
+    mask = numpy.ones(basic_image.shape, dtype=rasterio.bool_)
+    mask[4:5, 4:5] = False
+    truth = basic_image * numpy.invert(mask)
+
+    with rasterio.drivers():
+        sieved_image = sieve(basic_image, basic_image.sum(), mask=mask)
+        assert sieved_image.sum() > 0
+
+        assert numpy.array_equal(
+            truth,
+            sieved_image
+        )
+
+        assert numpy.array_equal(
+            truth.astype(rasterio.uint8),
+            sieved_image
+        )
+
+
+def test_sieve_blank_mask(basic_image):
+    """ A blank mask should have no effect """
+
+    mask = numpy.ones(basic_image.shape, dtype=rasterio.bool_)
+    with rasterio.drivers():
+        assert numpy.array_equal(
+            basic_image,
+            sieve(basic_image, basic_image.sum(), mask=mask)
+        )
+
+
+def test_sieve_invalid_mask_shape(basic_image):
+    """ A mask that is the wrong shape should fail """
+
+    with rasterio.drivers():
+        with pytest.raises(ValueError):
+            sieve(
+                basic_image, basic_image.sum(),
+                mask=numpy.ones(
+                    (basic_image.shape[0] + 10, basic_image.shape[1] + 10),
+                    dtype=rasterio.bool_
+                )
+            )
+
+
+def test_sieve_invalid_mask_dtype(basic_image):
+    """ A mask that is the wrong dtype should fail """
+
+    with rasterio.drivers():
+        for dtype in ('int8', 'int16', 'int32'):
+            with pytest.raises(ValueError):
+                sieve(
+                    basic_image, basic_image.sum(),
+                    mask=numpy.ones(basic_image.shape, dtype=dtype)
+                )
+
+
+def test_sieve_supported_dtypes(basic_image):
+    """ Supported data types should return valid results """
+
+    supported_types = (
+        ('int16', -32768),
+        ('int32', -2147483648),
+        ('uint8', 255),
+        ('uint16', 65535)
+    )
+
+    with rasterio.drivers():
+        for dtype, test_value in supported_types:
+            truth = (basic_image).astype(dtype) * test_value
+            sieved_image = sieve(truth, basic_image.sum())
+            assert numpy.array_equal(truth, sieved_image)
+            assert numpy.dtype(sieved_image.dtype) == numpy.dtype(dtype)
+
+
+def test_sieve_unsupported_dtypes(basic_image):
+    """ Unsupported data types should raise exceptions """
+
+    unsupported_types = (
+        ('int8', -127),
+        ('uint32', 4294967295),
+        ('int64', 20439845334323),
+        ('float16', -9343.232),
+        ('float32', 1.434532),
+        ('float64', -98332.133422114)
+    )
+
+    with rasterio.drivers():
+        for dtype, test_value in unsupported_types:
+            with pytest.raises(ValueError):
+                sieve(
+                    (basic_image).astype(dtype) * test_value,
+                    basic_image.sum()
+                )
+
+
+def test_sieve_band(pixelated_image, pixelated_image_file):
+    """ Sieving a band from a raster file should match sieve of array """
+
+    with rasterio.drivers():
+        truth = sieve(pixelated_image, 9)
+
+        with rasterio.open(pixelated_image_file) as src:
+            band = rasterio.band(src, 1)
+            assert numpy.array_equal(truth, sieve(band, 9))
+
+            # Mask band should also work but will be a no-op
+            assert numpy.array_equal(
+                pixelated_image,
+                sieve(band, 9, mask=band)
+            )
+
+
+def test_sieve_internal_driver_manager(basic_image, pixelated_image):
+    """ Sieve should work without explicitly calling driver manager """
+
+    assert numpy.array_equal(
+        basic_image,
+        sieve(pixelated_image, basic_image.sum())
+    )
diff --git a/tests/test_features_bounds.py b/tests/test_features_bounds.py
deleted file mode 100644
index a7015cd..0000000
--- a/tests/test_features_bounds.py
+++ /dev/null
@@ -1,65 +0,0 @@
-from rasterio.features import bounds
-
-
-# Tests copied from Fiona 1.4.1
-
-def test_bounds_point():
-    g = {'type': 'Point', 'coordinates': [10, 10]}
-    assert bounds(g) == (10, 10, 10, 10)
-
-
-def test_bounds_line():
-    g = {'type': 'LineString', 'coordinates': [[0, 0], [10, 10]]}
-    assert bounds(g) == (0, 0, 10, 10)
-
-
-def test_bounds_polygon():
-    g = {'type': 'Polygon', 'coordinates': [[[0, 0], [10, 10], [10, 0]]]}
-    assert bounds(g) == (0, 0, 10, 10)
-
-
-def test_bounds_z():
-    g = {'type': 'Point', 'coordinates': [10, 10, 10]}
-    assert bounds(g) == (10, 10, 10, 10)
-
-
-# TODO: add these to Fiona with update to bounds
-def test_bounds_existing_bbox():
-    """ Test with existing bbox in geojson, similar to that produced by
-    rasterio.  Values specifically modified here for testing, bboxes are not
-    valid as written.
-    """
-
-    fc = {
-        'bbox': [-107, 40, -105, 41],
-        'features': [{
-            'bbox': [-107, 40, -104, 42],
-            'geometry': {
-                'coordinates': [
-                    [[-107, 40], [-106, 40], [-106, 41], [-107, 41], [-107, 40]]
-                ],
-                'type': 'Polygon'
-            },
-            'type': 'Feature'
-        }],
-        'type': 'FeatureCollection'
-    }
-    assert bounds(fc['features'][0]) == (-107, 40, -104, 42)
-    assert bounds(fc) == (-107, 40, -105, 41)
-
-
-def test_feature_collection():
-    fc = {
-        'features': [{
-            'geometry': {
-                'coordinates': [
-                    [[-107, 40], [-106, 40], [-106, 41], [-107, 41], [-107, 40]]
-                ],
-                'type': 'Polygon'
-            },
-            'type': 'Feature'
-        }],
-        'type': 'FeatureCollection'
-    }
-    assert bounds(fc['features'][0]) == (-107, 40, -106, 41)
-    assert bounds(fc) == (-107, 40, -106, 41)
diff --git a/tests/test_features_rasterize.py b/tests/test_features_rasterize.py
deleted file mode 100644
index 8848a8f..0000000
--- a/tests/test_features_rasterize.py
+++ /dev/null
@@ -1,181 +0,0 @@
-import logging
-import sys
-import numpy
-import pytest
-
-import rasterio
-from rasterio.features import shapes, rasterize, geometry_mask
-
-
-logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
-
-
-def test_rasterize_geometries():
-    """
-    Make sure that geometries are correctly rasterized according to parameters
-    """
-
-    rows = cols = 10
-    transform = (1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
-    geometry = {
-        'type': 'Polygon',
-        'coordinates': [[(2, 2), (2, 4.25), (4.25, 4.25), (4.25, 2), (2, 2)]]
-    }
-
-    with rasterio.drivers():
-        # we expect a subset of the pixels using default mode
-        result = rasterize([geometry], out_shape=(rows, cols))
-        truth = numpy.zeros((rows, cols))
-        truth[2:4, 2:4] = 1
-        assert numpy.array_equal(result, truth)
-
-        out = numpy.zeros((rows, cols))
-        result = rasterize([geometry], out=out, default_value=1)
-        assert numpy.array_equal(out, truth)
-
-        # we expect all touched pixels
-        result = rasterize(
-            [geometry], out_shape=(rows, cols), all_touched=True
-        )
-        truth = numpy.zeros((rows, cols))
-        truth[2:5, 2:5] = 1
-        assert numpy.array_equal(result, truth)
-
-        # we expect the pixel value to match the one we pass in
-        value = 5
-        result = rasterize([(geometry, value)], out_shape=(rows, cols))
-        truth = numpy.zeros((rows, cols))
-        truth[2:4, 2:4] = value
-        assert numpy.array_equal(result, truth)
-
-        # Check the fill and default transform.
-        # we expect the pixel value to match the one we pass in
-        value = 5
-        result = rasterize(
-            [(geometry, value)],
-            out_shape=(rows, cols),
-            fill=1
-        )
-        truth = numpy.ones((rows, cols))
-        truth[2:4, 2:4] = value
-        assert numpy.array_equal(result, truth)
-
-
-def test_rasterize_dtype():
-    """Make sure that data types are handled correctly"""
-
-    rows = cols = 10
-    transform = (1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
-    geometry = {
-        'type': 'Polygon',
-        'coordinates': [[(2, 2), (2, 4.25), (4.25, 4.25), (4.25, 2), (2, 2)]]
-    }
-
-    with rasterio.drivers():
-        # Supported types should all work properly
-        supported_types = (
-            ('int16', -32768),
-            ('int32', -2147483648),
-            ('uint8', 255),
-            ('uint16', 65535),
-            ('uint32', 4294967295),
-            ('float32', 1.434532),
-            ('float64', -98332.133422114)
-        )
-
-        for dtype, default_value in supported_types:
-            truth = numpy.zeros((rows, cols), dtype=dtype)
-            truth[2:4, 2:4] = default_value
-
-            result = rasterize(
-                [geometry],
-                out_shape=(rows, cols),
-                default_value=default_value,
-                dtype=dtype
-            )
-            assert numpy.array_equal(result, truth)
-            assert numpy.dtype(result.dtype) == numpy.dtype(truth.dtype)
-
-            result = rasterize(
-                [(geometry, default_value)],
-                out_shape=(rows, cols)
-            )
-            if numpy.dtype(dtype).kind == 'f':
-                assert numpy.allclose(result, truth)
-            else:
-                assert numpy.array_equal(result, truth)
-            # Since dtype is auto-detected, it may not match due to upcasting
-
-        # Unsupported types should all raise exceptions
-        unsupported_types = (
-            ('int8', -127),
-            ('int64', 20439845334323),
-            ('float16', -9343.232)
-        )
-
-        for dtype, default_value in unsupported_types:
-            with pytest.raises(ValueError):
-                rasterize(
-                    [geometry],
-                    out_shape=(rows, cols),
-                    default_value=default_value,
-                    dtype=dtype
-                )
-
-            with pytest.raises(ValueError):
-                rasterize(
-                    [(geometry, default_value)],
-                    out_shape=(rows, cols),
-                    dtype=dtype
-                )
-
-        # Mismatched values and dtypes should raise exceptions
-        mismatched_types = (('uint8', 3.2423), ('uint8', -2147483648))
-        for dtype, default_value in mismatched_types:
-            with pytest.raises(ValueError):
-                rasterize(
-                    [geometry],
-                    out_shape=(rows, cols),
-                    default_value=default_value,
-                    dtype=dtype
-                )
-
-            with pytest.raises(ValueError):
-                rasterize(
-                    [(geometry, default_value)],
-                    out_shape=(rows, cols),
-                    dtype=dtype
-                )
-
-
-def test_rasterize_geometries_symmetric():
-    """Make sure that rasterize is symmetric with shapes"""
-
-    rows = cols = 10
-    transform = (1.0, 0.0, 0.0, 0.0, -1.0, 0.0)
-    truth = numpy.zeros((rows, cols), dtype=rasterio.ubyte)
-    truth[2:5, 2:5] = 1
-    with rasterio.drivers():
-        s = shapes(truth, transform=transform)
-        result = rasterize(s, out_shape=(rows, cols), transform=transform)
-        assert numpy.array_equal(result, truth)
-
-
-def test_geometry_mask():
-    rows = cols = 10
-    transform = (1.0, 0.0, 0.0, 0.0, -1.0, 0.0)
-    truth = numpy.zeros((rows, cols), dtype=rasterio.bool_)
-    truth[2:5, 2:5] = True
-    with rasterio.drivers():
-        s = shapes((truth * 10).astype(rasterio.ubyte), transform=transform)
-        # Strip out values returned from shapes, and only keep first shape
-        geoms = [next(s)[0]]
-
-        # Regular mask should be the inverse of truth raster
-        mask = geometry_mask(geoms, out_shape=(rows, cols), transform=transform)
-        assert numpy.array_equal(mask, numpy.invert(truth))
-
-        # Inverted mask should be the same as the truth raster
-        mask = geometry_mask(geoms, out_shape=(rows, cols), transform=transform,
-                             invert=True)
-        assert numpy.array_equal(mask, truth)
\ No newline at end of file
diff --git a/tests/test_features_shapes.py b/tests/test_features_shapes.py
deleted file mode 100644
index 2419a15..0000000
--- a/tests/test_features_shapes.py
+++ /dev/null
@@ -1,144 +0,0 @@
-import logging
-import sys
-import numpy
-import pytest
-
-import rasterio
-import rasterio.features as ftrz
-
-
-logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
-
-
-def test_shapes():
-    """Test creation of shapes from pixel values"""
-
-    image = numpy.zeros((20, 20), dtype=rasterio.ubyte)
-    image[5:15, 5:15] = 127
-    with rasterio.drivers():
-        shapes = ftrz.shapes(image)
-        shape, val = next(shapes)
-        assert shape['type'] == 'Polygon'
-        assert len(shape['coordinates']) == 2  # exterior and hole
-        assert val == 0
-        shape, val = next(shapes)
-        assert shape['type'] == 'Polygon'
-        assert len(shape['coordinates']) == 1  # no hole
-        assert val == 127
-        try:
-            shape, val = next(shapes)
-        except StopIteration:
-            assert True
-        else:
-            assert False
-
-
-def test_shapes_band_shortcut():
-    """Test rasterio bands as input to shapes"""
-
-    with rasterio.drivers():
-        with rasterio.open('tests/data/shade.tif') as src:
-            shapes = ftrz.shapes(rasterio.band(src, 1))
-            shape, val = next(shapes)
-            assert shape['type'] == 'Polygon'
-            assert len(shape['coordinates']) == 1
-            assert val == 255
-
-
-def test_shapes_internal_driver_manager():
-    """Make sure this works if driver is managed outside this test"""
-
-    image = numpy.zeros((20, 20), dtype=rasterio.ubyte)
-    image[5:15, 5:15] = 127
-    shapes = ftrz.shapes(image)
-    shape, val = next(shapes)
-    assert shape['type'] == 'Polygon'
-
-
-def test_shapes_connectivity():
-    """Test connectivity options"""
-
-    image = numpy.zeros((20, 20), dtype=rasterio.ubyte)
-    image[5:11, 5:11] = 1
-    image[11, 11] = 1
-
-    shapes = ftrz.shapes(image, connectivity=4)
-    shape, val = next(shapes)
-    assert len(shape['coordinates'][0]) == 5
-
-    shapes = ftrz.shapes(image, connectivity=8)
-    shape, val = next(shapes)
-    assert len(shape['coordinates'][0]) == 9
-    # Note: geometry is not technically valid at this point, it has a self
-    # intersection at 11,11
-
-    # Test invalid connectivity
-    with pytest.raises(ValueError):
-        shapes = ftrz.shapes(image, connectivity=12)
-        next(shapes)
-
-
-def test_shapes_dtype():
-    """Test image data type handling"""
-
-    rows = cols = 10
-    with rasterio.drivers():
-        supported_types = (
-            ('int16', -32768),
-            ('int32', -2147483648),
-            ('uint8', 255),
-            ('uint16', 65535),
-            ('float32', 1.434532)
-        )
-
-        for dtype, test_value in supported_types:
-            image = numpy.zeros((rows, cols), dtype=dtype)
-            image[2:5, 2:5] = test_value
-
-            shapes = ftrz.shapes(image)
-            shape, value = next(shapes)
-            if dtype == 'float32':
-                assert round(value, 6) == round(test_value, 6)
-            else:
-                assert value == test_value
-
-        # Unsupported types should all raise exceptions
-        unsupported_types = (
-            ('int8', -127),
-            ('uint32', 4294967295),
-            ('int64', 20439845334323),
-            ('float16', -9343.232),
-            ('float64', -98332.133422114)
-        )
-
-        for dtype, test_value in unsupported_types:
-            with pytest.raises(ValueError):
-                image = numpy.zeros((rows, cols), dtype=dtype)
-                image[2:5, 2:5] = test_value
-                shapes = ftrz.shapes(image)
-                next(shapes)
-
-        # Test mask types
-        image = numpy.zeros((rows, cols), dtype='uint8')
-        image.fill(255)
-        supported_mask_types = (
-            ('bool', 1),
-            ('uint8', 255)
-        )
-        for dtype, mask_value in supported_mask_types:
-            mask = numpy.zeros((rows, cols), dtype=dtype)
-            mask[2:5, 2:5] = mask_value
-            shapes = ftrz.shapes(image, mask=mask)
-            shape, value = next(shapes)
-            assert value == 255
-
-        unsupported_mask_types = (
-            ('int8', -127),
-            ('int16', -32768)
-        )
-        for dtype, mask_value in unsupported_mask_types:
-            with pytest.raises(ValueError):
-                mask = numpy.zeros((rows, cols), dtype=dtype)
-                mask[2:5, 2:5] = mask_value
-                shapes = ftrz.shapes(image, mask=mask)
-                next(shapes)
\ No newline at end of file
diff --git a/tests/test_features_sieve.py b/tests/test_features_sieve.py
deleted file mode 100644
index e34cf34..0000000
--- a/tests/test_features_sieve.py
+++ /dev/null
@@ -1,185 +0,0 @@
-import logging
-import sys
-import numpy
-import pytest
-
-import rasterio
-import rasterio.features as ftrz
-
-
-logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
-
-
-def test_sieve():
-    """Test sieving a 10x10 feature from an ndarray."""
-
-    image = numpy.zeros((20, 20), dtype=rasterio.ubyte)
-    image[5:15, 5:15] = 1
-
-    # An attempt to sieve out features smaller than 100 should not change the
-    # image.
-    with rasterio.drivers():
-        sieved_image = ftrz.sieve(image, 100)
-        assert numpy.array_equal(sieved_image, image)
-
-    # Setting the size to 100 should leave us an empty, False image.
-    with rasterio.drivers():
-        sieved_image = ftrz.sieve(image, 101)
-        assert not sieved_image.any()
-
-    # Invalid size value should fail
-    for invalid_size in (0, 45.1234, image.size + 1):
-        with pytest.raises(ValueError):
-            sieved_image = ftrz.sieve(image, invalid_size)
-
-
-def test_sieve_connectivity():
-    """Test proper behavior of connectivity"""
-
-    image = numpy.zeros((20, 20), dtype=rasterio.ubyte)
-    image[5:15:2, 5:15] = 1
-    image[6, 4] = 1
-    image[8, 15] = 1
-    image[10, 4] = 1
-    image[12, 15] = 1
-
-    # Diagonals not connected, all become small features that will be removed
-    sieved_image = ftrz.sieve(image, 54, connectivity=4)
-    assert not sieved_image.any()
-
-    # Diagonals connected, everything is retained
-    sieved_image = ftrz.sieve(image, 54, connectivity=8)
-    assert numpy.array_equal(sieved_image, image)
-
-    # Invalid connectivity value should fail
-    with pytest.raises(ValueError):
-        ftrz.sieve(image, 54, connectivity=12)
-
-
-def test_sieve_output():
-    """Test proper behavior of output image, if passed into sieve"""
-
-    with rasterio.drivers():
-        shape = (20, 20)
-        image = numpy.zeros(shape, dtype=rasterio.ubyte)
-        image[5:15, 5:15] = 1
-
-        # Output should match returned array
-        output = numpy.zeros_like(image)
-        output[1:3, 1:3] = 5
-        sieved_image = ftrz.sieve(image, 100, output=output)
-        assert numpy.array_equal(output, sieved_image)
-
-        # Output of different dtype should fail
-        output = numpy.zeros(shape, dtype=rasterio.int32)
-        with pytest.raises(ValueError):
-            ftrz.sieve(image, 100, output)
-
-        # Output of a different shape should fail
-        output = numpy.zeros((21, 21), dtype=rasterio.ubyte)
-        with pytest.raises(ValueError):
-            ftrz.sieve(image, 100, output)
-
-
-def test_sieve_mask():
-    """Test proper behavior of mask image, if passed int sieve"""
-
-    with rasterio.drivers():
-        shape = (20, 20)
-        image = numpy.zeros(shape, dtype=rasterio.ubyte)
-        image[5:15, 5:15] = 1
-        image[1:3, 1:3] = 2
-
-        # Blank mask has no effect, only areas smaller than size will be
-        # removed
-        mask = numpy.ones(shape, dtype=rasterio.bool_)
-        sieved_image = ftrz.sieve(image, 100, mask=mask)
-        truth = numpy.zeros_like(image)
-        truth[5:15, 5:15] = 1
-        assert numpy.array_equal(sieved_image, truth)
-
-        # Only areas within the overlap of the mask and values will be kept
-        mask = numpy.ones(shape, dtype=rasterio.bool_)
-        mask[7:10, 7:10] = False
-        sieved_image = ftrz.sieve(image, 100, mask=mask)
-        truth = numpy.zeros_like(image)
-        truth[7:10, 7:10] = 1
-        assert numpy.array_equal(sieved_image, truth)
-
-        with pytest.raises(ValueError):
-            mask = numpy.ones((21, 21), dtype=rasterio.bool_)
-            ftrz.sieve(image, 100, mask=mask)
-
-
-def test_dtypes():
-    """Test data type support for sieve"""
-
-    rows = cols = 10
-    with rasterio.drivers():
-        supported_types = (
-            ('int16', -32768),
-            ('int32', -2147483648),
-            ('uint8', 255),
-            ('uint16', 65535)
-        )
-
-        for dtype, test_value in supported_types:
-            image = numpy.zeros((rows, cols), dtype=dtype)
-            image[2:5, 2:5] = test_value
-
-            # Sieve should return the original image
-            sieved_image = ftrz.sieve(image, 2)
-            assert numpy.array_equal(image, sieved_image)
-            assert numpy.dtype(sieved_image.dtype).name == dtype
-
-            # Sieve should return a blank image
-            sieved_image = ftrz.sieve(image, 10)
-            assert numpy.array_equal(numpy.zeros_like(image), sieved_image)
-            assert numpy.dtype(sieved_image.dtype).name == dtype
-
-        # Unsupported types should all raise exceptions
-        unsupported_types = (
-            ('int8', -127),
-            ('uint32', 4294967295),
-            ('int64', 20439845334323),
-            ('float16', -9343.232),
-            ('float32', 1.434532),
-            ('float64', -98332.133422114)
-        )
-
-        for dtype, test_value in unsupported_types:
-            with pytest.raises(ValueError):
-                image = numpy.zeros((rows, cols), dtype=dtype)
-                image[2:5, 2:5] = test_value
-                ftrz.sieve(image, 2)
-
-        # Test mask types
-        image = numpy.zeros((rows, cols), dtype='uint8')
-        image.fill(255)
-
-        supported_mask_types = (
-            ('bool', 1),
-            ('uint8', 255)
-        )
-        for dtype, mask_value in supported_mask_types:
-            mask = numpy.zeros((rows, cols), dtype=dtype)
-            mask[2:5, 2:5] = mask_value
-            sieved_image = ftrz.sieve(image, 2, mask=mask)
-            assert numpy.array_equal(image, sieved_image)
-
-        unsupported_mask_types = (
-            ('int8', -127),
-            ('int16', -32768)
-        )
-        for dtype, mask_value in unsupported_mask_types:
-            with pytest.raises(ValueError):
-                mask = numpy.zeros((rows, cols), dtype=dtype)
-                mask[2:5, 2:5] = mask_value
-                ftrz.sieve(image, 2, mask=mask)
-
-
-def test_sieve_shade():
-    with rasterio.drivers():
-        with rasterio.open('tests/data/shade.tif') as src:
-            sieved_image = ftrz.sieve(rasterio.band(src, 1), 42)
-            assert sieved_image.shape == (1024, 1024)
diff --git a/tests/test_indexing.py b/tests/test_indexing.py
index 257e7dc..d42ca16 100644
--- a/tests/test_indexing.py
+++ b/tests/test_indexing.py
@@ -1,4 +1,13 @@
+import numpy
+import pytest
+
 import rasterio
+from rasterio import (
+    get_data_window, window_intersection, window_union, windows_intersect
+)
+
+
+DATA_WINDOW = ((3, 5), (2, 6))
 
 
 def test_index():
@@ -67,3 +76,119 @@ def test_window_full_cover():
         win = src.window(*bounds)
         bounds_calc = list(src.window_bounds(win))
         assert bound_covers(bounds_calc, bounds)
+
+
+ at pytest.fixture
+def data():
+    data = numpy.zeros((10, 10), dtype='uint8')
+    data[slice(*DATA_WINDOW[0]), slice(*DATA_WINDOW[1])] = 1
+    return data
+
+
+def test_data_window_unmasked(data):
+    window = get_data_window(data)
+    assert window == ((0, data.shape[0]), (0, data.shape[1]))
+
+
+def test_data_window_masked(data):
+    data = numpy.ma.masked_array(data, data == 0)
+    window = get_data_window(data)
+    assert window == DATA_WINDOW
+
+
+def test_data_window_nodata(data):
+    window = get_data_window(data, nodata=0)
+    assert window == DATA_WINDOW
+
+    window = get_data_window(numpy.ones_like(data), nodata=0)
+    assert window == ((0, data.shape[0]), (0, data.shape[1]))
+
+
+def test_data_window_nodata_disjunct():
+    data = numpy.zeros((3, 10, 10), dtype='uint8')
+    data[0, :4, 1:4] = 1
+    data[1, 2:5, 2:8] = 1
+    data[2, 1:6, 1:6] = 1
+    window = get_data_window(data, nodata=0)
+    assert window == ((0, 6), (1, 8))
+
+
+def test_data_window_empty_result():
+    data = numpy.zeros((3, 10, 10), dtype='uint8')
+    window = get_data_window(data, nodata=0)
+    assert window == ((0, 0), (0, 0))
+
+
+def test_data_window_masked_file():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        window = get_data_window(src.read(1, masked=True))
+        assert window == ((3, 714), (13, 770))
+
+        window = get_data_window(src.read(masked=True))
+        assert window == ((3, 714), (13, 770))
+
+
+def test_window_union():
+    assert window_union([
+        ((0, 6), (3, 6)),
+        ((2, 4), (1, 5))
+    ]) == ((0, 6), (1, 6))
+
+
+def test_window_intersection():
+    assert window_intersection([
+        ((0, 6), (3, 6)),
+        ((2, 4), (1, 5))
+    ]) == ((2, 4), (3, 5))
+
+    assert window_intersection([
+        ((0, 6), (3, 6)),
+        ((6, 10), (1, 5))
+    ]) == ((6, 6), (3, 5))
+
+    assert window_intersection([
+        ((0, 6), (3, 6)),
+        ((2, 4), (1, 5)),
+        ((3, 6), (0, 6))
+    ]) == ((3, 4), (3, 5))
+
+
+def test_window_intersection_disjunct():
+    with pytest.raises(ValueError):
+        window_intersection([
+            ((0, 6), (3, 6)),
+            ((100, 200), (0, 12)),
+            ((7, 12), (7, 12))
+        ])
+
+
+def test_windows_intersect():
+    assert windows_intersect([
+        ((0, 6), (3, 6)),
+        ((2, 4), (1, 5))
+    ]) == True
+
+    assert windows_intersect([
+        ((0, 6), (3, 6)),
+        ((2, 4), (1, 5)),
+        ((3, 6), (0, 6))
+    ]) == True
+
+
+def test_windows_intersect_disjunct():
+    assert windows_intersect([
+        ((0, 6), (3, 6)),
+        ((10, 20), (0, 6))
+    ]) == False
+
+    assert windows_intersect([
+        ((0, 6), (3, 6)),
+        ((2, 4), (1, 5)),
+        ((5, 6), (0, 6))
+    ]) == False
+
+    assert windows_intersect([
+        ((0, 6), (3, 6)),
+        ((2, 4), (1, 3)),
+        ((3, 6), (4, 6))
+    ]) == False
\ No newline at end of file
diff --git a/tests/test_mask_creation.py b/tests/test_mask_creation.py
index a6f8c3f..77c6c60 100644
--- a/tests/test_mask_creation.py
+++ b/tests/test_mask_creation.py
@@ -5,6 +5,7 @@ See https://github.com/mapbox/rasterio/issues/293 for bug report.
 """
 
 import rasterio
+from rasterio.enums import MaskFlags
 
 
 def test_create_internal_mask(data):
@@ -22,6 +23,10 @@ def test_create_internal_mask(data):
     # Check that the mask was saved correctly.
     with rasterio.open(str(data.join('RGB.byte.tif'))) as src:
         assert (mask == src.read_mask()).all()
+        for flags in src.mask_flags:
+            assert flags & MaskFlags.per_dataset
+            assert not flags & MaskFlags.alpha
+            assert not flags & MaskFlags.nodata
 
 
 def test_create_sidecar_mask(data):
@@ -39,6 +44,10 @@ def test_create_sidecar_mask(data):
     # Check that the mask was saved correctly.
     with rasterio.open(str(data.join('RGB.byte.tif'))) as src:
         assert (mask == src.read_mask()).all()
+        for flags in src.mask_flags:
+            assert flags & MaskFlags.per_dataset
+            assert not flags & MaskFlags.alpha
+            assert not flags & MaskFlags.nodata
 
     # Check the .msk file, too.
     with rasterio.open(str(data.join('RGB.byte.tif.msk'))) as msk:
diff --git a/tests/test_profile.py b/tests/test_profile.py
index 07cc2cf..b85b885 100644
--- a/tests/test_profile.py
+++ b/tests/test_profile.py
@@ -60,14 +60,29 @@ def test_gtiff_profile_protected_driver():
 
 
 def test_open_with_profile(tmpdir):
-        tiffname = str(tmpdir.join('foo.tif'))
-        with rasterio.open(
-                tiffname,
-                'w',
-                **default_gtiff_profile(
-                    count=1, width=1, height=1)) as dst:
-            data = dst.read()
-            assert data.flatten().tolist() == [0]
+    tiffname = str(tmpdir.join('foo.tif'))
+    with rasterio.open(
+            tiffname,
+            'w',
+            **default_gtiff_profile(
+                count=1, width=256, height=256)) as dst:
+        data = dst.read()
+
+
+def test_blockxsize_guard(tmpdir):
+    """blockxsize can't be greater than image width."""
+    tiffname = str(tmpdir.join('foo.tif'))
+    with pytest.raises(ValueError):
+        _ = rasterio.open(tiffname, 'w', **default_gtiff_profile(
+                count=1, width=128, height=256))
+
+
+def test_blockysize_guard(tmpdir):
+    """blockysize can't be greater than image height."""
+    tiffname = str(tmpdir.join('foo.tif'))
+    with pytest.raises(ValueError):
+        _ = rasterio.open(tiffname, 'w', **default_gtiff_profile(
+                count=1, width=256, height=128))
 
 
 def test_profile_overlay():
diff --git a/tests/test_rio_features.py b/tests/test_rio_features.py
index a1a4b15..b1f20a0 100644
--- a/tests/test_rio_features.py
+++ b/tests/test_rio_features.py
@@ -3,559 +3,738 @@ import os
 import re
 import sys
 import numpy
+import json
+from affine import Affine
 
 import rasterio
 from rasterio.rio import features
 
 
+DEFAULT_SHAPE = (10, 10)
+
 logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
 
-TEST_FEATURES = """{
-    "geometry": {
-        "coordinates": [
-            [
-                [-110, 40],
-                [-100, 40],
-                [-100, 45],
-                [-105, 45],
-                [-110, 40]
-            ]
-        ],
-        "type": "Polygon"
-    },
-    "properties": {
-        "val": 15
-    },
-    "type": "Feature"
-}"""
-
-TEST_MERC_FEATURES = """{
-  "geometry": {
-      "coordinates": [
-          [
-              [-11858134, 4808920],
-              [-11868134, 4804143],
-              [-11853357, 4804143],
-              [-11840000, 4812000],
-              [-11858134, 4808920]
-          ]
-      ],
-      "type": "Polygon"
-  },
-  "properties": {
-      "val": 10
-  },
-  "type": "Feature"
-}"""
-
-# > rio shapes tests/data/shade.tif --mask --sampling 500 --projected --precision 0
-TEST_MERC_FEATURECOLLECTION = """{
-    "bbox": [-11858135.0, 4803914.0, -11848351.0, 4813698.0],
-    "features": [{
-        "bbox": [-11853357.504145855, 4808920.97837715,
-                 -11848580.189878704, 4813698.2926443005],
-        "geometry": {
-            "coordinates": [
-                [
-                    [-11853357.504145855, 4813698.2926443005],
-                    [-11853357.504145855, 4808920.97837715],
-                    [-11848580.189878704, 4808920.97837715],
-                    [-11848580.189878704, 4813698.2926443005],
-                    [-11853357.504145855, 4813698.2926443005]
-                ]
-            ],
-            "type": "Polygon"
-        },
-        "properties": {
-            "val": 2
-        },
-        "type": "Feature"
-    }, {
-        "bbox": [-11858134.818413004, 4804143.66411,
-                 -11853357.504145855, 4808920.97837715],
-        "geometry": {
-            "coordinates": [
-                [
-                    [-11858134.818413004, 4808920.97837715],
-                    [-11858134.818413004, 4804143.66411],
-                    [-11853357.504145855, 4804143.66411],
-                    [-11853357.504145855, 4808920.97837715],
-                    [-11858134.818413004, 4808920.97837715]
-                ]
-            ],
-            "type": "Polygon"
-        },
-        "properties": {
-            "val": 3
-        },
-        "type": "Feature"
-    }],
-    "type": "FeatureCollection"
-}"""
-
-
-def test_mask(runner, tmpdir):
+
+def test_mask(runner, tmpdir, basic_feature, basic_image_2x2,
+              pixelated_image_file):
+
     output = str(tmpdir.join('test.tif'))
 
-    with rasterio.open('tests/data/shade.tif') as src:
-        src_data = src.read(1, masked=True)
+    result = runner.invoke(
+        features.mask,
+        [pixelated_image_file, output, '--geojson-mask', '-'],
+        input=json.dumps(basic_feature)
+    )
 
-        result = runner.invoke(
-            features.mask,
-            ['tests/data/shade.tif', output, '--geojson-mask', '-'],
-            input=TEST_MERC_FEATURES
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+
+    with rasterio.open(output) as out:
+        assert numpy.array_equal(
+            basic_image_2x2,
+            out.read(1, masked=True).filled(0)
         )
-        assert result.exit_code == 0
-        assert os.path.exists(output)
-
-        masked_count = 0
-        with rasterio.open(output) as out:
-            assert out.count == src.count
-            assert out.shape == src.shape
-            out_data = out.read(1, masked=True)
-
-            # Make sure that pixels with 0 value were converted to mask
-            masked_count = (src_data == 0).sum() - (out_data == 0).sum()
-            assert masked_count == 79743
-            assert out_data.mask.sum() - src_data.mask.sum() == masked_count
-
-        # Test using --all-touched option
-        result = runner.invoke(
-            features.mask,
-            [
-                'tests/data/shade.tif', output,
-                '--all',
-                '--geojson-mask', '-'
-            ],
-            input=TEST_MERC_FEATURES
+
+
+def test_mask_all_touched(runner, tmpdir, basic_feature, basic_image,
+                          pixelated_image_file):
+
+    output = str(tmpdir.join('test.tif'))
+
+    result = runner.invoke(
+        features.mask,
+        [pixelated_image_file, output, '--all', '--geojson-mask', '-'],
+        input=json.dumps(basic_feature)
+    )
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+
+    with rasterio.open(output) as out:
+        assert numpy.array_equal(
+            basic_image,
+            out.read(1, masked=True).filled(0)
         )
-        assert result.exit_code == 0
-        with rasterio.open(output) as out:
-            out_data = out.read(1, masked=True)
-
-            # Make sure that more pixels with 0 value were converted to mask
-            masked_count2 = (src_data == 0).sum() - (out_data == 0).sum()
-            assert masked_count2 > 0 and masked_count > masked_count2
-            assert out_data.mask.sum() - src_data.mask.sum() == masked_count2
-
-        # Test using --invert option
-        result = runner.invoke(
-            features.mask,
-            [
-                'tests/data/shade.tif', output,
-                '--invert',
-                '--geojson-mask', '-'
-            ],
-            input=TEST_MERC_FEATURES
+
+
+def test_mask_invert(runner, tmpdir, basic_feature, pixelated_image,
+                     pixelated_image_file):
+
+    truth = pixelated_image
+    truth[2:4, 2:4] = 0
+
+    output = str(tmpdir.join('test.tif'))
+
+    result = runner.invoke(
+        features.mask,
+        [pixelated_image_file, output, '--invert', '--geojson-mask', '-'],
+        input=json.dumps(basic_feature)
+    )
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+
+    with rasterio.open(output) as out:
+        assert numpy.array_equal(
+            truth,
+            out.read(1, masked=True).filled(0)
         )
-        assert result.exit_code == 0
-        with rasterio.open(output) as out:
-            out_data = out.read(1, masked=True)
-            # Areas that were masked when not inverted should now be 0
-            assert (out_data == 0).sum() == masked_count
 
-    # Test with feature collection
+
+def test_mask_featurecollection(runner, tmpdir, basic_featurecollection,
+                                 basic_image_2x2, pixelated_image_file):
+
+    output = str(tmpdir.join('test.tif'))
+
     result = runner.invoke(
         features.mask,
-        ['tests/data/shade.tif', output, '--geojson-mask', '-'],
-        input=TEST_MERC_FEATURECOLLECTION
+        [pixelated_image_file, output, '--geojson-mask', '-'],
+        input=json.dumps(basic_featurecollection)
     )
     assert result.exit_code == 0
+    assert os.path.exists(output)
+
+    with rasterio.open(output) as out:
+        assert numpy.array_equal(
+            basic_image_2x2,
+            out.read(1, masked=True).filled(0)
+        )
+
+
+def test_mask_out_of_bounds(runner, tmpdir, basic_feature,
+                            pixelated_image_file):
+    """
+    A GeoJSON mask that is outside bounds of raster should result in a
+    blank image.
+    """
+
+    coords = numpy.array(basic_feature['geometry']['coordinates']) - 10
+    basic_feature['geometry']['coordinates'] = coords.tolist()
+
+    output = str(tmpdir.join('test.tif'))
 
-    # Missing GeoJSON should make copy of input to output
-    output2 = str(tmpdir.join('test2.tif'))
     result = runner.invoke(
         features.mask,
-        ['tests/data/shade.tif', output2]
+        [pixelated_image_file, output, '--geojson-mask', '-'],
+        input=json.dumps(basic_feature)
     )
     assert result.exit_code == 0
-    assert os.path.exists(output2)
-    with rasterio.open('tests/data/shade.tif') as src:
-        with rasterio.open(output2) as out:
-            src_data = src.read(1, masked=True)
-            out_data = out.read(1, masked=True)
-            assert numpy.array_equal(src_data, out_data)
+    assert 'outside bounds' in result.output
+    assert os.path.exists(output)
+
+    with rasterio.open(output) as out:
+        assert not numpy.any(out.read(1, masked=True).filled(0))
+
+
+def test_mask_no_geojson(runner, tmpdir, pixelated_image, pixelated_image_file):
+    """ Mask without geojson input should simply return same raster as input """
+
+    output = str(tmpdir.join('test.tif'))
+
+    result = runner.invoke(features.mask, [pixelated_image_file, output])
+    assert result.exit_code == 0
+    assert os.path.exists(output)
 
-    # Invalid JSON should fail
+    with rasterio.open(output) as out:
+        assert numpy.array_equal(
+            pixelated_image,
+            out.read(1, masked=True).filled(0)
+        )
+
+
+def test_mask_invalid_geojson(runner, tmpdir, pixelated_image_file):
+    """ Invalid GeoJSON should fail """
+
+    output = str(tmpdir.join('test.tif'))
+
+    # Using invalid JSON
     result = runner.invoke(
         features.mask,
-        ['tests/data/shade.tif', output, '--geojson-mask', '-'],
+        [pixelated_image_file, output, '--geojson-mask', '-'],
         input='{bogus: value}'
     )
     assert result.exit_code == 2
     assert 'GeoJSON could not be read' in result.output
 
+    # Using invalid GeoJSON
     result = runner.invoke(
         features.mask,
-        ['tests/data/shade.tif', output, '--geojson-mask', '-'],
+        [pixelated_image_file, output, '--geojson-mask', '-'],
         input='{"bogus": "value"}'
     )
     assert result.exit_code == 2
     assert 'Invalid GeoJSON' in result.output
 
 
-def test_mask_crop(runner, tmpdir):
+def test_mask_crop(runner, tmpdir, basic_feature, pixelated_image):
+    """
+    In order to test --crop option, we need to use a transform more similar to
+    a normal raster, with a negative y pixel size.
+    """
+
+    image = pixelated_image
+    outfilename = str(tmpdir.join('pixelated_image.tif'))
+    kwargs = {
+        "crs": {'init': 'epsg:4326'},
+        "transform": Affine(1, 0, 0, 0, -1, 0),
+        "count": 1,
+        "dtype": rasterio.uint8,
+        "driver": "GTiff",
+        "width": image.shape[1],
+        "height": image.shape[0],
+        "nodata": 255
+    }
+    with rasterio.drivers():
+        with rasterio.open(outfilename, 'w', **kwargs) as out:
+            out.write_band(1, image)
+
+
     output = str(tmpdir.join('test.tif'))
 
-    with rasterio.open('tests/data/shade.tif') as src:
+    truth = numpy.zeros((4, 3))
+    truth[1:3, 0:2] = 1
+
+    result = runner.invoke(
+        features.mask,
+        [outfilename, output, '--crop', '--geojson-mask', '-'],
+        input=json.dumps(basic_feature)
+    )
 
-        result = runner.invoke(
-            features.mask,
-            [
-                'tests/data/shade.tif', output,
-                '--crop',
-                '--geojson-mask', '-'
-            ],
-            input=TEST_MERC_FEATURES
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert numpy.array_equal(
+            truth,
+            out.read(1, masked=True).filled(0)
         )
-        assert result.exit_code == 0
-        assert os.path.exists(output)
-        with rasterio.open(output) as out:
-            assert out.shape[1] == src.shape[1]
-            assert out.shape[0] < src.shape[0]
-            assert out.shape[0] == 824
 
-    # Adding invert option after crop should be ignored
+
+def test_mask_crop_inverted_y(runner, tmpdir, basic_feature, pixelated_image_file):
+    """
+    --crop option should also work if raster has a positive y pixel size
+    (e.g., Affine.identity() ).
+    """
+
+    output = str(tmpdir.join('test.tif'))
+
+    truth = numpy.zeros((4, 3))
+    truth[1:3, 0:2] = 1
+
     result = runner.invoke(
         features.mask,
-        [
-            'tests/data/shade.tif', output,
-            '--crop',
-            '--invert',
-            '--geojson-mask', '-'
-        ],
-        input=TEST_MERC_FEATURES
+        [pixelated_image_file, output, '--crop', '--geojson-mask', '-'],
+        input=json.dumps(basic_feature)
     )
+
     assert result.exit_code == 0
-    assert 'Invert option ignored' in result.output
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert numpy.array_equal(
+            truth,
+            out.read(1, masked=True).filled(0)
+        )
+
 
+def test_mask_crop_out_of_bounds(runner, tmpdir, basic_feature,
+                                 pixelated_image_file):
+    """
+    A GeoJSON mask that is outside bounds of raster should fail with
+    --crop option.
+    """
+
+    coords = numpy.array(basic_feature['geometry']['coordinates']) - 10
+    basic_feature['geometry']['coordinates'] = coords.tolist()
 
-def test_mask_out_of_bounds(runner, tmpdir):
     output = str(tmpdir.join('test.tif'))
-    # Crop with out of bounds raster should
+
     result = runner.invoke(
         features.mask,
-        ['tests/data/shade.tif', output, '--geojson-mask', '-'],
-        input=TEST_FEATURES
+        [pixelated_image_file, output, '--crop', '--geojson-mask', '-'],
+        input=json.dumps(basic_feature)
     )
-    assert result.exit_code == 0
-    assert 'outside bounds' in result.output
+    assert result.exit_code == 2
+    assert 'not allowed' in result.output
+
+
+def test_mask_crop_and_invert(runner, tmpdir, basic_feature, pixelated_image,
+                              pixelated_image_file):
+    """ Adding crop and invert options should ignore invert option """
+
+    output = str(tmpdir.join('test.tif'))
 
-    # Crop with out of bounds raster should fail
     result = runner.invoke(
         features.mask,
         [
-            'tests/data/shade.tif', output,
+            pixelated_image_file, output,
             '--crop',
+            '--invert',
             '--geojson-mask', '-'
         ],
-        input=TEST_FEATURES
+        input=json.dumps(basic_feature)
     )
-    assert result.exit_code == 2
-    assert 'not allowed' in result.output
+    assert result.exit_code == 0
+    assert 'Invert option ignored' in result.output
 
 
-def test_shapes(runner):
-    result = runner.invoke(features.shapes, ['tests/data/shade.tif'])
+def test_shapes(runner, pixelated_image_file):
+    result = runner.invoke(features.shapes, [pixelated_image_file])
+
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
-    assert result.output.count('"Feature"') == 232
+    assert result.output.count('"Feature"') == 4
+    assert numpy.allclose(
+        json.loads(result.output)['features'][0]['geometry']['coordinates'],
+        [[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]]
+    )
+
+
+def test_shapes_invalid_bidx(runner, pixelated_image_file):
+    result = runner.invoke(features.shapes, [pixelated_image_file, '--bidx', 4])
 
-    # Invalid band index should fail
-    result = runner.invoke(
-        features.shapes, ['tests/data/shade.tif', '--bidx', '4'])
     assert result.exit_code == 1
+    # Underlying exception message trapped by shapes
+
 
+def test_shapes_sequence(runner, pixelated_image_file):
+    """
+    --sequence option should produce 4 features in series rather than
+    inside a feature collection.
+    """
+
+    result = runner.invoke(features.shapes, [pixelated_image_file, '--sequence'])
 
-def test_shapes_sequence(runner):
-    result = runner.invoke(features.shapes, ['tests/data/shade.tif', '--sequence'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 0
-    assert result.output.count('"Feature"') == 232
+    assert result.output.count('"Feature"') == 4
+    assert result.output.count('\n') == 4
+
 
+def test_shapes_sequence_rs(runner, pixelated_image_file):
+    """ --rs option should use the feature separator character. """
 
-def test_shapes_sequence_rs(runner):
     result = runner.invoke(
-        features.shapes, [
-            'tests/data/shade.tif',
-            '--sequence',
-            '--rs'])
+        features.shapes, [pixelated_image_file, '--sequence', '--rs']
+    )
+
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 0
-    assert result.output.count('"Feature"') == 232
-    assert result.output.count(u'\u001e') == 232
+    assert result.output.count('"Feature"') == 4
+    assert result.output.count(u'\u001e') == 4
+
+
 
+def test_shapes_with_nodata(runner, pixelated_image, pixelated_image_file):
+    """
+    An area of nodata should also be represented with a shape when using
+    --with-nodata option
+    """
 
-def test_shapes_with_nodata(runner):
-    result = runner.invoke(features.shapes, ['tests/data/shade.tif', '--with-nodata'])
+    pixelated_image[0:2, 8:10] = 255
+
+    with rasterio.open(pixelated_image_file, 'r+') as out:
+        out.write_band(1, pixelated_image)
+
+    result = runner.invoke(
+        features.shapes, [pixelated_image_file, '--with-nodata']
+    )
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
-    assert result.output.count('"Feature"') == 288
+    assert result.output.count('"Feature"') == 5
+
 
+def test_shapes_indent(runner, pixelated_image_file):
+    """
+    --indent option should produce lots of newlines and contiguous spaces
+    """
+
+    result = runner.invoke(
+        features.shapes, [pixelated_image_file, '--indent', 2]
+    )
 
-def test_shapes_indent(runner):
-    result = runner.invoke(features.shapes, ['tests/data/shade.tif', '--indent', '2'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
-    assert result.output.count('\n') == 70371
+    assert result.output.count('"Feature"') == 4
+    assert result.output.count('\n') == 231
+    assert result.output.count('        ') == 180
+
 
+def test_shapes_compact(runner, pixelated_image_file):
+    result = runner.invoke(features.shapes, [pixelated_image_file, '--compact'])
 
-def test_shapes_compact(runner):
-    result = runner.invoke(features.shapes, ['tests/data/shade.tif', '--compact'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
+    assert result.output.count('"Feature"') == 4
     assert result.output.count(', ') == 0
     assert result.output.count(': ') == 0
 
 
-def test_shapes_sampling(runner):
+def test_shapes_sampling(runner, pixelated_image_file):
+    """ --sampling option should remove the single pixel features """
     result = runner.invoke(
-        features.shapes, ['tests/data/shade.tif', '--sampling', '11'])
+        features.shapes, [pixelated_image_file, '--sampling', 2]
+    )
+
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
-    assert result.output.count('"Feature"') == 117
+    assert result.output.count('"Feature"') == 2
 
 
-def test_shapes_precision(runner):
+def test_shapes_precision(runner, pixelated_image_file):
+    """ Output numbers should have no more than 1 decimal place """
+
     result = runner.invoke(
-        features.shapes, ['tests/data/shade.tif', '--precision', '1'])
+        features.shapes, [pixelated_image_file, '--precision', 1]
+    )
+
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
-    # Find no numbers with 2+ decimal places.
+    assert result.output.count('"Feature"') == 4
     assert re.search(r'\d*\.\d{2,}', result.output) is None
 
 
-def test_shapes_mask(runner):
-    result = runner.invoke(features.shapes, ['tests/data/RGB.byte.tif', '--mask'])
-    assert result.exit_code == 0
-    assert result.output.count('"FeatureCollection"') == 1
-    assert result.output.count('"Feature"') == 7
+def test_shapes_mask(runner, pixelated_image, pixelated_image_file):
+    """ --mask should extract the nodata area of the image """
 
+    pixelated_image[0:5, 0:10] = 255
+    pixelated_image[0:10, 0:3] = 255
+    pixelated_image[8:10, 8:10] = 255
+
+    with rasterio.open(pixelated_image_file, 'r+') as out:
+        out.write_band(1, pixelated_image)
+
+    result = runner.invoke(features.shapes, [pixelated_image_file, '--mask'])
+
+    print(result.output)
+    print(result.exception)
 
-def test_shapes_mask_decimated(runner):
-    result = runner.invoke(
-        features.shapes, 
-        ['tests/data/RGB.byte.tif', '--mask', '--sampling', '10'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
     assert result.output.count('"Feature"') == 1
 
+    assert numpy.allclose(
+        json.loads(result.output)['features'][0]['geometry']['coordinates'],
+        [[[3, 5], [3, 10], [8, 10], [8, 8], [9, 8], [10, 8], [10, 5], [3, 5]]]
+    )
+
+
+def test_shapes_mask_sampling(runner, pixelated_image, pixelated_image_file):
+    """
+    using --sampling with the mask should snap coordinates to the nearest
+    factor of 5
+    """
+    pixelated_image[0:5, 0:10] = 255
+    pixelated_image[0:10, 0:3] = 255
+    pixelated_image[8:10, 8:10] = 255
+
+    with rasterio.open(pixelated_image_file, 'r+') as out:
+        out.write_band(1, pixelated_image)
+
+    result = runner.invoke(
+        features.shapes, [pixelated_image_file, '--mask', '--sampling', 5]
+    )
 
-def test_shapes_band1_as_mask(runner):
-    result = runner.invoke(features.shapes,
-        ['tests/data/RGB.byte.tif', '--band', '--bidx', '1', '--as-mask'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
-    assert result.output.count('"Feature"') == 9
+    assert result.output.count('"Feature"') == 1
 
+    assert numpy.allclose(
+        json.loads(result.output)['features'][0]['geometry']['coordinates'],
+        [[[5, 5], [5, 10], [10, 10], [10, 5], [5, 5]]]
+    )
 
-def test_rasterize_err(tmpdir, runner):
-    output = str(tmpdir.join('test.tif'))
-    # Test invalid stdin
-    result = runner.invoke(features.rasterize, [output], input='BOGUS')
-    assert result.exit_code == -1
 
-    # Test invalid GeoJSON
-    result = runner.invoke(features.rasterize, [output],
-                           input='{"foo": "bar"}')
-    assert result.exit_code == 2
+def test_shapes_band1_as_mask(runner, pixelated_image, pixelated_image_file):
+    """
+    When using --as-mask option, pixel value should not matter, only depends
+    on pixels being contiguous.
+    """
 
-    # Test invalid res
-    result = runner.invoke(features.rasterize, [output], input=TEST_FEATURES)
-    assert result.exit_code == 2
+    pixelated_image[2:3, 2:3] = 4
 
-    # Test invalid CRS for bounds
-    result = runner.invoke(features.rasterize, [output, '--res', 1],
-                           input=TEST_MERC_FEATURECOLLECTION)
-    assert result.exit_code == 2
+    with rasterio.open(pixelated_image_file, 'r+') as out:
+        out.write_band(1, pixelated_image)
 
-    # Test invalid CRS value
-    result = runner.invoke(features.rasterize, [output,
-                                                '--res', 1,
-                                                '--src-crs', 'BOGUS'],
-                           input=TEST_MERC_FEATURECOLLECTION)
-    assert result.exit_code == 2
+    result = runner.invoke(
+        features.shapes,
+        [pixelated_image_file, '--band', '--bidx', '1', '--as-mask']
+    )
 
+    assert result.exit_code == 0
+    assert result.output.count('"FeatureCollection"') == 1
+    assert result.output.count('"Feature"') == 3
+    assert numpy.allclose(
+        json.loads(result.output)['features'][1]['geometry']['coordinates'],
+        [[[2, 2], [2, 5], [5, 5], [5, 2], [2, 2]]]
+    )
 
-def test_rasterize(tmpdir, runner):
-    # Test dimensions
+
+def test_rasterize(tmpdir, runner, basic_feature):
     output = str(tmpdir.join('test.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output, '--dimensions', 20, 10],
-                           input=TEST_FEATURES)
+    result = runner.invoke(
+        features.rasterize,
+        [output, '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1]],
+        input=json.dumps(basic_feature)
+    )
+
     assert result.exit_code == 0
     assert os.path.exists(output)
     with rasterio.open(output) as out:
-        assert out.count == 1
-        assert out.meta['width'] == 20
-        assert out.meta['height'] == 10
+        assert numpy.allclose(out.bounds, (2, 2, 4.25, 4.25))
         data = out.read(1, masked=False)
-        assert (data == 0).sum() == 55
-        assert (data == 1).sum() == 145
-
-    # Test dimensions and bounds
-    output = str(tmpdir.join('test2.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output,
-                            '--dimensions', 40, 20,
-                            '--bounds', -120, 30, -90, 50
-                           ], input=TEST_FEATURES)
+        assert data.shape == DEFAULT_SHAPE
+        assert numpy.all(data)
+
+
+def test_rasterize_bounds(tmpdir, runner, basic_feature, basic_image_2x2):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        features.rasterize,
+        [
+            output,
+            '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1],
+            '--bounds', 0, 10, 10, 0
+        ],
+        input=json.dumps(basic_feature)
+    )
+
     assert result.exit_code == 0
     assert os.path.exists(output)
     with rasterio.open(output) as out:
-        assert out.count == 1
-        assert out.meta['width'] == 40
-        assert out.meta['height'] == 20
+        assert numpy.allclose(out.bounds, (0, 10, 10, 0))
         data = out.read(1, masked=False)
-        assert (data == 0).sum() == 748
-        assert (data == 1).sum() == 52
+        assert numpy.array_equal(basic_image_2x2, data)
+        assert data.shape == DEFAULT_SHAPE
+
+
+def test_rasterize_resolution(tmpdir, runner, basic_feature):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        features.rasterize,
+        [output, '--res', 0.15],
+        input=json.dumps(basic_feature)
+    )
 
-    # Test resolution
-    output = str(tmpdir.join('test3.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output, '--res', 0.5], input=TEST_FEATURES)
     assert result.exit_code == 0
     assert os.path.exists(output)
     with rasterio.open(output) as out:
-        assert out.count == 1
-        assert out.meta['width'] == 20
-        assert out.meta['height'] == 10
+        assert numpy.allclose(out.bounds, (2, 2, 4.25, 4.25))
         data = out.read(1, masked=False)
-        assert (data == 0).sum() == 55
-        assert (data == 1).sum() == 145
-
-    # Test that src-crs is written into new output
-    output = str(tmpdir.join('test4.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output,
-                            '--dimensions', 20, 10,
-                            '--src-crs', 'EPSG:3857'
-                           ],
-                           input=TEST_MERC_FEATURECOLLECTION)
+        assert data.shape == (15, 15)
+        assert numpy.all(data)
+
+
+def test_rasterize_src_crs(tmpdir, runner, basic_feature):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        features.rasterize,
+        [
+            output,
+            '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1],
+            '--src-crs', 'EPSG:3857'
+        ],
+        input=json.dumps(basic_feature)
+    )
+
     assert result.exit_code == 0
+    assert os.path.exists(output)
     with rasterio.open(output) as out:
         assert out.crs['init'].lower() == 'epsg:3857'
 
 
-def test_rasterize_existing_output(tmpdir, runner):
+def test_rasterize_mismatched_src_crs(tmpdir, runner, basic_feature):
+    """
+    A --src-crs that is geographic with coordinates that are outside
+    world bounds should fail.
+    """
+
+    coords = numpy.array(basic_feature['geometry']['coordinates']) * 100000
+    basic_feature['geometry']['coordinates'] = coords.tolist()
+
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        features.rasterize,
+        [
+            output,
+            '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1],
+            '--src-crs', 'EPSG:4326'
+        ],
+        input=json.dumps(basic_feature)
+    )
+
+    assert result.exit_code == 2
+    assert 'Bounds are beyond the valid extent for EPSG:4326' in result.output
+
+
+def test_rasterize_invalid_src_crs(tmpdir, runner, basic_feature):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        features.rasterize,
+        [
+            output,
+            '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1],
+            '--src-crs', 'foo:bar'
+        ],
+        input=json.dumps(basic_feature)
+    )
+
+    assert result.exit_code == 2
+    assert 'invalid CRS.  Must be an EPSG code.' in result.output
+
+
+def test_rasterize_existing_output(tmpdir, runner, basic_feature):
+    """
+    Create a rasterized output, then rasterize additional pixels into it.
+    The final result should include rasterized pixels from both features.
+    """
+
+    truth = numpy.zeros(DEFAULT_SHAPE)
+    truth[2:4, 2:4] = 1
+    truth[4:6, 4:6] = 1
+
     output = str(tmpdir.join('test.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output, '--res', 0.5], input=TEST_FEATURES)
+    result = runner.invoke(
+        features.rasterize,
+        [
+            output,
+            '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1],
+            '--bounds', 0, 10, 10, 0
+        ],
+        input=json.dumps(basic_feature)
+    )
+
     assert result.exit_code == 0
     assert os.path.exists(output)
 
-    geojson = """{
-        "geometry": {
-            "coordinates": [
-                [
-                    [-102, 40],
-                    [-98, 40],
-                    [-98, 45],
-                    [-100, 45],
-                    [-102, 40]
-                ]
-            ],
-            "type": "Polygon"
-        },
-        "type": "Feature"
-    }"""
-
-    result = runner.invoke(features.rasterize, [output, '--default-value', 2],
-                           input=geojson)
+    coords = numpy.array(basic_feature['geometry']['coordinates']) + 2
+    basic_feature['geometry']['coordinates'] = coords.tolist()
+
+    result = runner.invoke(
+        features.rasterize,
+        [output, '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1]],
+        input=json.dumps(basic_feature)
+    )
+
+    assert result.exit_code == 0
 
     with rasterio.open(output) as out:
-        assert out.count == 1
-        data = out.read(1, masked=False)
-        assert (data == 0).sum() == 55
-        assert (data == 1).sum() == 125
-        assert (data == 2).sum() == 20
-
-    # Confirm that a different src-crs is rejected, even if a geographic crs
-    result = runner.invoke(features.rasterize,
-                           [output,
-                            '--res', 0.5,
-                            '--src-crs', 'EPSG:4269'
-                            ], input=TEST_FEATURES)
-    assert result.exit_code == 2
+        assert numpy.array_equal(truth, out.read(1, masked=False))
+
 
+def test_rasterize_like_raster(tmpdir, runner, basic_feature, basic_image_2x2,
+                               pixelated_image_file):
 
-def test_rasterize_like(tmpdir, runner):
     output = str(tmpdir.join('test.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output, '--like', 'tests/data/shade.tif'],
-                           input=TEST_MERC_FEATURECOLLECTION)
+
+    result = runner.invoke(
+        features.rasterize,
+        [output, '--like', pixelated_image_file],
+        input=json.dumps(basic_feature)
+    )
+
     assert result.exit_code == 0
     assert os.path.exists(output)
     with rasterio.open(output) as out:
-        assert out.count == 1
-        data = out.read(1, masked=False)
-        assert (data == 0).sum() == 548576
-        assert (data == 1).sum() == 500000
+        assert numpy.array_equal(basic_image_2x2, out.read(1, masked=False))
+
+        with rasterio.open(pixelated_image_file) as src:
+            assert out.crs == src.crs
+            assert out.bounds == src.bounds
+            assert src.affine == src.affine
+
+
+def test_rasterize_invalid_like_raster(tmpdir, runner, basic_feature):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        features.rasterize,
+        [output, '--like', str(tmpdir.join('foo.tif'))],
+        input=json.dumps(basic_feature)
+    )
 
-    # Test invalid like raster
-    output = str(tmpdir.join('test2.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output, '--like', str(tmpdir.join('foo.tif'))], input=TEST_FEATURES)
     assert result.exit_code == 2
+    assert 'Invalid value for "--like":' in result.output
+
+
+def test_rasterize_like_raster_src_crs_mismatch(tmpdir, runner, basic_feature,
+                                                pixelated_image_file):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        features.rasterize,
+        [output, '--like', pixelated_image_file, '--src-crs', 'EPSG:3857'],
+        input=json.dumps(basic_feature)
+    )
 
-    # Test that src-crs different than --like raster crs breaks
-    output = str(tmpdir.join('test3.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output,
-                            '--like', 'tests/data/shade.tif',
-                            '--src-crs', 'EPSG:4326'],
-                           input=TEST_FEATURES)
     assert result.exit_code == 2
+    assert 'GeoJSON does not match crs of --like raster' in result.output
 
 
-def test_rasterize_property_value(tmpdir, runner):
-    # Test feature collection property values
+def test_rasterize_property_value(tmpdir, runner, basic_feature):
     output = str(tmpdir.join('test.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output,
-                            '--res', 1000,
-                            '--property', 'val',
-                            '--src-crs', 'EPSG:3857'
-                           ],
-                           input=TEST_MERC_FEATURECOLLECTION)
-    assert result.exit_code == 0
-    assert os.path.exists(output)
-    with rasterio.open(output) as out:
-        assert out.count == 1
-        data = out.read(1, masked=False)
-        assert (data == 0).sum() == 50
-        assert (data == 2).sum() == 25
-        assert (data == 3).sum() == 25
-
-    # Test feature property values
-    output = str(tmpdir.join('test2.tif'))
-    result = runner.invoke(features.rasterize,
-                           [output, '--res', 0.5, '--property', 'val'],
-                           input=TEST_FEATURES)
+    result = runner.invoke(
+        features.rasterize,
+        [
+            output,
+            '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1],
+            '--property', 'val'
+        ],
+        input=json.dumps(basic_feature)
+    )
+
     assert result.exit_code == 0
     assert os.path.exists(output)
     with rasterio.open(output) as out:
-        assert out.count == 1
+        assert numpy.allclose(out.bounds, (2, 2, 4.25, 4.25))
         data = out.read(1, masked=False)
-        assert (data == 0).sum() == 55
-        assert (data == 15).sum() == 145
+        assert data.shape == DEFAULT_SHAPE
+        assert numpy.all(data == basic_feature['properties']['val'])
+
 
+def test_rasterize_like_raster_outside_bounds(tmpdir, runner, basic_feature,
+                                              pixelated_image_file):
+    """
+    Rasterizing a feature outside bounds of --like raster should result
+    in a blank image
+    """
+
+    coords = numpy.array(basic_feature['geometry']['coordinates']) + 100
+    basic_feature['geometry']['coordinates'] = coords.tolist()
 
-def test_rasterize_out_of_bounds(tmpdir, runner):
     output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        features.rasterize,
+        [output, '--like', pixelated_image_file],
+        input=json.dumps(basic_feature)
+    )
 
-    # Test out of bounds of --like raster
-    result = runner.invoke(features.rasterize,
-                           [output, '--like', 'tests/data/shade.tif'],
-                           input=TEST_FEATURES)
     assert result.exit_code == 0
     assert 'outside bounds' in result.output
     assert os.path.exists(output)
     with rasterio.open(output) as out:
-        data = out.read_band(1, masked=False)
-        assert data.sum() == 0
+        assert not numpy.any(out.read_band(1, masked=False))
 
-    # Confirm that this does not fail when out of bounds for existing raster
-    result = runner.invoke(features.rasterize, [output], input=TEST_FEATURES)
-    assert result.exit_code == 0
-    assert 'outside bounds' in result.output
+
+def test_rasterize_invalid_stdin(tmpdir, runner):
+    """ Invalid value for stdin should fail with exception """
+
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(features.rasterize, [output], input='BOGUS')
+
+    assert result.exit_code == -1
+
+
+def test_rasterize_invalid_geojson(tmpdir, runner):
+    """ Invalid GeoJSON should fail with error  """
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(features.rasterize, [output], input='{"A": "B"}')
+
+    assert result.exit_code == 2
+    assert 'Invalid GeoJSON' in result.output
+
+
+def test_rasterize_missing_parameters(tmpdir, runner, basic_feature):
+    """ At least --res or --dimensions are required """
+
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        features.rasterize, [output], input=json.dumps(basic_feature)
+    )
+
+    assert result.exit_code == 2
+    assert 'pixel dimensions are required' in result.output
diff --git a/tests/test_rio_options.py b/tests/test_rio_options.py
new file mode 100644
index 0000000..edb7245
--- /dev/null
+++ b/tests/test_rio_options.py
@@ -0,0 +1,63 @@
+import os.path
+import uuid
+
+import click
+import pytest
+
+from rasterio.rio.options import file_in_handler
+
+
+class MockContext:
+
+    def __init__(self):
+        self.obj = {}
+
+
+class MockOption:
+
+    def __init__(self, name):
+        self.name = name
+
+
+def test_file_in_handler_no_vfs_nonexistent():
+    """file does not exist"""
+    ctx = MockContext()
+    with pytest.raises(click.BadParameter):
+        file_in_handler(ctx, 'INPUT', '{0}.tif'.format(uuid.uuid4()))
+
+
+def test_file_in_handler_no_vfs():
+    """file path is expanded to abspath"""
+    ctx = MockContext()
+    retval = file_in_handler(ctx, 'INPUT', 'tests/data/RGB.byte.tif')
+    assert retval == os.path.abspath('tests/data/RGB.byte.tif')
+
+
+def test_file_in_handler_bad_scheme():
+    """lolwut scheme is not supported"""
+    ctx = MockContext()
+    with pytest.raises(click.BadParameter):
+        _ = file_in_handler(ctx, 'INPUT', 'lolwut://bogus')
+
+
+def test_file_in_handler_with_vfs_nonexistent():
+    """archive does not exist"""
+    ctx = MockContext()
+    with pytest.raises(click.BadParameter):
+        _ = file_in_handler(
+            ctx, 'INPUT', 'zip://{0}/files.zip!/inputs/RGB.byte.tif'.format(uuid.uuid4()))
+
+
+def test_file_in_handler_with_vfs():
+    """vfs file path path is expanded"""
+    ctx = MockContext()
+    retval = file_in_handler(ctx, 'INPUT', 'zip://tests/data/files.zip!/inputs/RGB.byte.tif')
+    assert retval.startswith('zip:///')
+    assert retval.endswith('tests/data/files.zip!/inputs/RGB.byte.tif')
+
+
+def test_file_in_handler_with_vfs_file():
+    """vfs file path path is expanded"""
+    ctx = MockContext()
+    retval = file_in_handler(ctx, 'INPUT', 'file://tests/data/RGB.byte.tif')
+    assert retval.endswith('tests/data/RGB.byte.tif')
diff --git a/tests/test_vfs.py b/tests/test_vfs.py
new file mode 100644
index 0000000..0041a55
--- /dev/null
+++ b/tests/test_vfs.py
@@ -0,0 +1,101 @@
+import logging
+import sys
+
+import pytest
+
+import rasterio
+from rasterio.profiles import default_gtiff_profile
+from rasterio.vfs import parse_path, vsi_path
+
+
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+
+def test_parse_path_with_vfs():
+    """Correctly parse path with legacy vfs parameter"""
+    assert parse_path('foo.tif', vfs='zip://tests/data/files.zip') == (
+        'foo.tif', 'tests/data/files.zip', 'zip')
+
+
+def test_parse_path_zip():
+    """Correctly parse VFS scheme URL"""
+    assert parse_path('zip://tests/data/files.zip!foo.tif') == (
+        'foo.tif', 'tests/data/files.zip', 'zip')
+
+
+def test_parse_path_file_scheme():
+    """Correctly parse file:// URL"""
+    assert parse_path('file://foo.tif') == (
+        'foo.tif', None, 'file')
+
+
+def test_parse_path_file():
+    """Correctly parse an ordinary filesystem path"""
+    assert parse_path('/foo.tif') == (
+        '/foo.tif', None, None)
+
+
+def test_parse_unknown_scheme():
+    """Raise exception for unknown WFS scheme"""
+    with pytest.raises(ValueError):
+        parse_path('http://foo.tif')
+
+
+def test_vsi_path_scheme():
+    """Correctly make a vsi path"""
+    assert vsi_path(
+        'foo.tif', 'tests/data/files.zip', 'zip') == '/vsizip/tests/data/files.zip/foo.tif'
+
+
+def test_vsi_path_file():
+    """Correctly make a ordinary file path from a parsed file:// URL"""
+    assert vsi_path(
+        'foo.tif', None, 'file') == 'foo.tif'
+
+
+def test_vsi_path_file():
+    """Correctly make and ordinary file path from a file path"""
+    assert vsi_path(
+        'foo.tif', None, 'file') == 'foo.tif'
+
+
+def test_read_vfs_zip():
+    with rasterio.open(
+            'zip://tests/data/files.zip!/RGB.byte.tif') as src:
+        assert src.name == 'zip://tests/data/files.zip!/RGB.byte.tif'
+        assert src.count == 3
+
+
+def test_read_vfs_file():
+    with rasterio.open(
+            'file://tests/data/RGB.byte.tif') as src:
+        assert src.name == 'file://tests/data/RGB.byte.tif'
+        assert src.count == 3
+
+def test_read_vfs_zip_cmp_array():
+    with rasterio.open(
+            'zip://tests/data/files.zip!/RGB.byte.tif') as src:
+        zip_arr = src.read()
+
+    with rasterio.open(
+            'file://tests/data/RGB.byte.tif') as src:
+        file_arr = src.read()
+
+    assert zip_arr.dumps() == file_arr.dumps()
+
+
+def test_read_vfs_none():
+    with rasterio.open(
+            'tests/data/RGB.byte.tif') as src:
+        assert src.name == 'tests/data/RGB.byte.tif'
+        assert src.count == 3
+
+
+ at pytest.mark.parametrize('mode', ['r+', 'w'])
+def test_update_vfs(tmpdir, mode):
+    """VFS datasets can not be created or updated"""
+    with pytest.raises(TypeError):
+        _ = rasterio.open(
+            'zip://{0}'.format(tmpdir), mode,
+            **default_gtiff_profile(
+                count=1, width=1, height=1))
diff --git a/tests/test_warnings.py b/tests/test_warnings.py
new file mode 100644
index 0000000..c2e6b06
--- /dev/null
+++ b/tests/test_warnings.py
@@ -0,0 +1,8 @@
+from rasterio.warnings import NodataShadowWarning
+
+
+def test_nodata_shadow():
+    assert str(NodataShadowWarning()) == (
+        "The dataset's nodata attribute is shadowing "
+        "the alpha band. All masks will be determined "
+        "by the nodata attribute")

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/rasterio.git



More information about the Pkg-grass-devel mailing list