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

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

commit a423ade4b57ab4e66353d69215da323d871a0582
Author: Johan Van de Wauw <jo...@vandewauw.be>
Date:   Tue Sep 29 08:47:53 2015 +0200

    Imported Upstream version 0.27.0
---
 .travis.yml                   |  29 +++++--
 CHANGES.txt                   |  16 ++++
 docs/cli.rst                  |  26 +++++++
 rasterio/__init__.py          |   2 +-
 rasterio/_base.pyx            | 110 +++++++++++++++++++++++---
 rasterio/_io.pyx              |  16 ++--
 rasterio/coords.py            |  12 +++
 rasterio/crs.py               |   2 +-
 rasterio/dtypes.py            |  10 +++
 rasterio/enums.py             |  17 ++++
 rasterio/rio/convert.py       |  92 +++++++++++++++++++++-
 rasterio/rio/features.py      |  28 ++++---
 rasterio/rio/info.py          |   2 +-
 rasterio/rio/merge.py         | 175 ++++++++----------------------------------
 rasterio/tool.py              |  71 +++++++++++++++++
 rasterio/tools/__init__.py    |   0
 rasterio/tools/merge.py       | 164 +++++++++++++++++++++++++++++++++++++++
 rasterio/transform.py         |   9 +++
 requirements-dev.txt          |   2 +-
 setup.py                      |   1 +
 tests/data/rgb1.tif           | Bin 0 -> 481148 bytes
 tests/data/rgb2.tif           | Bin 0 -> 471548 bytes
 tests/data/rgb3.tif           | Bin 0 -> 383844 bytes
 tests/data/rgb4.tif           | Bin 0 -> 376188 bytes
 tests/data/rgb_deflate.tif    | Bin 0 -> 912290 bytes
 tests/test_crs.py             |   9 +++
 tests/test_image_structure.py |  74 ++++++++++++++++++
 tests/test_indexing.py        |  24 +++++-
 tests/test_read_boundless.py  |  15 ++++
 tests/test_read_resample.py   |  17 ++--
 tests/test_rio_convert.py     |  65 +++++++++++++++-
 tests/test_rio_features.py    |   6 +-
 tests/test_rio_merge.py       |  92 +++++++++++++++++-----
 tests/test_tool.py            |  23 +++++-
 34 files changed, 885 insertions(+), 224 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index bbf5cc5..955a27d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,16 +1,33 @@
 language: python
+sudo: false
+cache:
+  directories:
+    - ~/.cache/pip
+env:
+  global:
+    - PIP_WHEEL_DIR=$HOME/.cache/pip/wheels
+    - PIP_FIND_LINKS=file://$HOME/.cache/pip/wheels
+addons:
+  apt:
+    packages:
+    - libgdal1h
+    - gdal-bin
+    - libgdal-dev
+    - libatlas-dev
+    - libatlas-base-dev
+    - gfortran
 python:
   - "2.7"
   - "3.3"
   - "3.4"
 before_install:
-  - sudo add-apt-repository -y ppa:ubuntugis/ppa
-  - sudo apt-get update -qq
-  - sudo apt-get install -y libgdal1h gdal-bin libgdal-dev
+  - pip install -U pip
+  - pip install wheel
 install:
-  - pip install --install-option="--no-cython-compile" cython==0.22
-  - pip install -r requirements-dev.txt
-  - pip install -e .
+  - "pip wheel -r requirements-dev.txt"
+  - "pip install -r requirements-dev.txt"
+  - "pip install coveralls"
+  - "pip install -e ."
 script: 
   - py.test --cov rasterio --cov-report term-missing
 after_success:
diff --git a/CHANGES.txt b/CHANGES.txt
index 96ee33a..01ea66e 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,6 +1,22 @@
 Changes
 =======
 
+0.27.0 (2015-09-25)
+-------------------
+- Ensure local uniqueness of the rio-shapes feature ids (#479).
+- Surface compression and interleaving as dataset properties and in rio-info
+  (#481). In the module, these are enums (`enums.Compression` and 
+  `enums.Interleaving`); the values of the enums correspond to GDAL terms
+  (i.e, "DEFLATE") and the names are what surface in the CLI ("deflate").
+- Change get_window() and DatasetReader.window() to return a window guaranteed
+  to cover the input bounding box (#464, #485).
+- Bug fix for improperly computed transforms of output file in tools.merge and
+  rio-merge (#485).
+- More robust determination of dataset nodata values. In particular, the 
+  absence of a nodata value is much more clear: dataset.nodata should never
+  return an out of range value when there is no nodata value, it should always
+  return `None` (#485).
+
 0.26.0 (2015-08-11)
 -------------------
 - Add dependency on click-plugins, a new project that takes over the plugin
diff --git a/docs/cli.rst b/docs/cli.rst
index 56f560c..86e591e 100644
--- a/docs/cli.rst
+++ b/docs/cli.rst
@@ -19,6 +19,7 @@ Rasterio's command line interface is a program named "rio".
     Commands:
       bounds     Write bounding boxes to stdout as GeoJSON.
       calc       Raster data calculator.
+      clip       Clip a raster to given bounds.
       convert    Copy and convert raster dataset.
       edit-info  Edit dataset metadata.
       env        Print information about the rio environment.
@@ -157,6 +158,31 @@ efficiently in Python.
 Please see `calc.rst <calc.rst>`__ for more details.
 
 
+clip
+----
+
+Added in 0.27
+
+The ``clip`` command clips a raster using bounds input directly or from a
+template raster.
+
+.. code-block:: console
+
+    $ rio clip input.tif output.tif --bounds xmin ymin xmax ymax
+    $ rio clip input.tif output.tif --like template.tif
+
+If using ``--bounds``, values must be in coordinate reference system of input.
+If using ``--like``, bounds will automatically be transformed to match the
+coordinate reference system of the input.
+
+It can also be combined to read bounds of a feature dataset using Fiona:
+
+.. code-block:: console
+
+    $ rio clip input.tif output.tif --bounds $(fio info features.shp --bounds)
+
+
+
 convert
 -------
 
diff --git a/rasterio/__init__.py b/rasterio/__init__.py
index 7cf55cf..11934fd 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.26.0"
+__version__ = "0.27.0"
 
 log = logging.getLogger('rasterio')
 class NullHandler(logging.Handler):
diff --git a/rasterio/_base.pyx b/rasterio/_base.pyx
index caf5e85..422db77 100644
--- a/rasterio/_base.pyx
+++ b/rasterio/_base.pyx
@@ -15,7 +15,7 @@ from rasterio._err import cpl_errs
 from rasterio import dtypes
 from rasterio.coords import BoundingBox
 from rasterio.transform import Affine
-from rasterio.enums import ColorInterp
+from rasterio.enums import ColorInterp, Compression, Interleaving
 
 
 log = logging.getLogger('rasterio')
@@ -277,7 +277,7 @@ cdef class DatasetReader(object):
     def get_nodatavals(self):
         cdef void *hband = NULL
         cdef double nodataval
-        cdef int success
+        cdef int success = 0
 
         if not self._nodatavals:
             if self._hds == NULL:
@@ -286,11 +286,23 @@ cdef class DatasetReader(object):
                 hband = _gdal.GDALGetRasterBand(self._hds, i+1)
                 if hband == NULL:
                     raise ValueError("Null band")
+                dtype = dtypes.dtype_fwd[_gdal.GDALGetRasterDataType(hband)]
                 nodataval = _gdal.GDALGetRasterNoDataValue(hband, &success)
                 val = nodataval
-                if not success:
+                # GDALGetRasterNoDataValue() has two ways of telling you that
+                # there's no nodata value. The success flag might come back
+                # 0 (FALSE). Even if it comes back 1 (TRUE), you still need
+                # to check that the return value is within the range of the
+                # data type. If so, the band has a nodata value. If not,
+                # there's no nodata value.
+                if (success == 0 or
+                        val < dtypes.dtype_ranges[dtype][0] or
+                        val > dtypes.dtype_ranges[dtype][1]):
                     val = None
+                log.debug("Nodata success: %d", success)
+                log.debug("Nodata value: %f", nodataval)
                 self._nodatavals.append(val)
+
         return self._nodatavals
 
     property nodatavals:
@@ -377,16 +389,15 @@ cdef class DatasetReader(object):
             row += self.height
         return c+a*col, f+e*row
 
-    def index(self, x, y):
+    def index(self, x, y, op=math.floor):
         """Returns the (row, col) index of the pixel containing (x, y)."""
-        a, b, c, d, e, f, _, _, _ = self.affine
-        return int(math.floor((y-f)/e)), int(math.floor((x-c)/a))
+        return get_index(x, y, self.affine)
 
     def window(self, left, bottom, right, top, boundless=False):
         """Returns the window corresponding to the world bounding box.
         If boundless is False, window is limited to extent of this dataset."""
-        EPS = 1.0e-8
-        window = tuple(zip(self.index(left + EPS, top - EPS), self.index(right 
+ EPS, bottom - EPS)))
+
+        window = get_window(left, bottom, right, top, self.affine)
         if boundless:
             return window
         else:
@@ -421,6 +432,21 @@ cdef class DatasetReader(object):
         self._read = True
         return m
 
+    @property
+    def compression(self):
+        val = self.tags(ns='IMAGE_STRUCTURE').get('COMPRESSION')
+        if val:
+            return Compression(val)
+        else:
+            return None
+
+    @property
+    def interleaving(self):
+        val = self.tags(ns='IMAGE_STRUCTURE').get('INTERLEAVE')
+        if val:
+            return Interleaving(val)
+        else:
+            return None
 
     property profile:
         """Basic metadata and creation options of this dataset.
@@ -435,9 +461,12 @@ cdef class DatasetReader(object):
                 blockxsize=self.block_shapes[0][1],
                 blockysize=self.block_shapes[0][0],
                 tiled=self.block_shapes[0][1] != self.width)
+            if self.compression:
+                m['compress'] = self.compression.name
+            if self.interleaving:
+                m['interleave'] = self.interleaving.name
             return m
 
-
     def lnglat(self):
         w, s, e, n = self.bounds
         cx = (w + e)/2.0
@@ -699,6 +728,69 @@ cpdef eval_window(object window, int height, int width):
     return (r_start, r_stop), (c_start, c_stop)
 
 
+def get_index(x, y, affine, op=math.floor, precision=6):
+    """
+    Returns the (row, col) index of the pixel containing (x, y) given a
+    coordinate reference system.
+
+    Parameters
+    ----------
+    x : float
+        x value in coordinate reference system
+    y : float
+        y value in coordinate reference system
+    affine : tuple
+        Coefficients mapping pixel coordinates to coordinate reference system.
+    op : function
+        Function to convert fractional pixels to whole numbers (floor, ceiling,
+        round)
+    precision : int
+        Decimal places of precision in indexing, as in `round()`.
+
+    Returns
+    -------
+    row : int
+        row index
+    col : int
+        col index
+    """
+    # Use an epsilon, magnitude determined by the precision parameter
+    # and sign determined by the op function: positive for floor, negative
+    # for ceil.
+    eps = 10.0**-precision * (1.0 - 2.0*op(0.1))
+    row = int(op((y - eps - affine[5]) / affine[4]))
+    col = int(op((x + eps - affine[2]) / affine[0]))
+    return row, col
+
+
+def get_window(left, bottom, right, top, affine, precision=6):
+    """
+    Returns a window tuple given coordinate bounds and the coordinate reference
+    system.
+
+    Parameters
+    ----------
+    left : float
+        Left edge of window
+    bottom : float
+        Bottom edge of window
+    right : float
+        Right edge of window
+    top : float
+        top edge of window
+    affine : tuple
+        Coefficients mapping pixel coordinates to coordinate reference system.
+    precision : int
+        Decimal places of precision in indexing, as in `round()`.
+    """
+    window_start = get_index(
+        left, top, affine, op=math.floor, precision=precision)
+    window_stop = get_index(
+        right, bottom, affine, op=math.ceil, precision=precision)
+    window = tuple(zip(window_start, window_stop))
+    return window
+
+
 def window_shape(window, height=-1, width=-1):
     """Returns shape of a window.
 
diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx
index c317678..1d42712 100644
--- a/rasterio/_io.pyx
+++ b/rasterio/_io.pyx
@@ -816,10 +816,10 @@ cdef class RasterReader(_base.DatasetReader):
         else:
             # Compute the overlap between the dataset and the boundless window.
             overlap = ((
-                max(min(window[0][0] or 0, self.height), 0),
-                max(min(window[0][1] or self.height, self.height), 0)), (
-                max(min(window[1][0] or 0, self.width), 0),
-                max(min(window[1][1] or self.width, self.width), 0)))
+                max(min(window[0][0], self.height), 0),
+                max(min(window[0][1], self.height), 0)), (
+                max(min(window[1][0], self.width), 0),
+                max(min(window[1][1], self.width), 0)))
 
             if overlap != ((0, 0), (0, 0)):
                 # Prepare a buffer.
@@ -978,10 +978,10 @@ cdef class RasterReader(_base.DatasetReader):
         else:
             # Compute the overlap between the dataset and the boundless window.
             overlap = ((
-                max(min(window[0][0] or 0, self.height), 0),
-                max(min(window[0][1] or self.height, self.height), 0)), (
-                max(min(window[1][0] or 0, self.width), 0),
-                max(min(window[1][1] or self.width, self.width), 0)))
+                max(min(window[0][0], self.height), 0),
+                max(min(window[0][1], self.height), 0)), (
+                max(min(window[1][0], self.width), 0),
+                max(min(window[1][1], self.width), 0)))
 
             if overlap != ((0, 0), (0, 0)):
                 # Prepare a buffer.
diff --git a/rasterio/coords.py b/rasterio/coords.py
index 1d89156..50eb59d 100644
--- a/rasterio/coords.py
+++ b/rasterio/coords.py
@@ -3,3 +3,15 @@ from collections import namedtuple
 
 BoundingBox = namedtuple('BoundingBox', ('left', 'bottom', 'right', 'top'))
 
+
+def disjoint_bounds(bounds1, bounds2):
+    """Returns True if bounds do not overlap
+
+    Parameters
+    ----------
+    bounds1: rasterio bounds tuple (xmin, ymin, xmax, ymax)
+    bounds2: rasterio bounds tuple
+    """
+
+    return (bounds1[0] > bounds2[2] or bounds1[2] < bounds2[0] or
+            bounds1[1] > bounds2[3] or bounds1[3] < bounds2[1])
\ No newline at end of file
diff --git a/rasterio/crs.py b/rasterio/crs.py
index 775a54f..2ee4d45 100644
--- a/rasterio/crs.py
+++ b/rasterio/crs.py
@@ -58,7 +58,7 @@ def from_string(prjs):
         except ValueError:
             raise ValueError('crs appears to be JSON but is not valid')
 
-    if 'EPSG:' in prjs.upper():
+    if prjs.strip().upper().startswith('EPSG:'):
         return from_epsg(prjs.split(':')[1])
 
     parts = [o.lstrip('+') for o in prjs.strip().split()]
diff --git a/rasterio/dtypes.py b/rasterio/dtypes.py
index 449ec5f..90097d8 100644
--- a/rasterio/dtypes.py
+++ b/rasterio/dtypes.py
@@ -58,6 +58,16 @@ typename_fwd = {
 
 typename_rev = dict((v, k) for k, v in typename_fwd.items())
 
+dtype_ranges = {
+    'uint8': (0, 255),
+    'uint16': (0, 65535),
+    'int16': (-32768, 32767),
+    'uint32': (0, 4294967295),
+    'int32': (-2147483648, 2147483647),
+    'float32': (-3.4028235e+38, 3.4028235e+38),
+    'float64': (-1.7976931348623157e+308, 1.7976931348623157e+308)}
+
+
 def _gdal_typename(dt):
     try:
         return typename_fwd[dtype_rev[dt]]
diff --git a/rasterio/enums.py b/rasterio/enums.py
index d33f71d..3e347a5 100644
--- a/rasterio/enums.py
+++ b/rasterio/enums.py
@@ -28,3 +28,20 @@ class Resampling(Enum):
     mode='MODE'
     average_magphase='AVERAGE_MAGPHASE'
     none='NONE'
+
+
+class Compression(Enum):
+    jpeg='JPEG'
+    lzw='LZW'
+    packbits='PACKBITS'
+    deflate='DEFLATE'
+    ccittrle='CCITTRLE'
+    ccittfax3='CCITTFAX3'
+    ccittfax4='CCITTFAX4'
+    lzma='LZMA'
+    none='NONE'
+
+
+class Interleaving(Enum):
+    pixel='PIXEL'
+    band='BAND'
diff --git a/rasterio/rio/convert.py b/rasterio/rio/convert.py
index 3443b51..11781f7 100644
--- a/rasterio/rio/convert.py
+++ b/rasterio/rio/convert.py
@@ -4,17 +4,107 @@ import logging
 import warnings
 
 import click
-from cligj import files_inout_arg, format_opt
+from cligj import format_opt
 import numpy as np
 
 from .helpers import resolve_inout
 from . import options
 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(
+    'files',
+    nargs=-1,
+    type=click.Path(resolve_path=True),
+    required=True,
+    metavar="INPUT OUTPUT")
+@options.output_opt
+@options.bounds_opt
+@click.option(
+    '--like',
+    type=click.Path(exists=True),
+    help='Raster dataset to use as a template for bounds')
+@format_opt
+@options.creation_options
+@click.pass_context
+def clip(
+        ctx,
+        files,
+        output,
+        bounds,
+        like,
+        driver,
+        creation_options):
+    """Clips a raster using bounds input directly or from a template raster.
+
+    \b
+      $ rio clip input.tif output.tif --bounds xmin ymin xmax ymax
+      $ rio clip input.tif output.tif --like template.tif
+
+    If using --bounds, values must be in coordinate reference system of input.
+    If using --like, bounds will automatically be transformed to match the
+    coordinate reference system of the input.
+
+    It can also be combined to read bounds of a feature dataset using Fiona:
+
+    \b
+      $ rio clip input.tif output.tif --bounds $(fio info features.shp 
--bounds)
+
+    """
+
+    from rasterio.warp import transform_bounds
+
+    verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+
+    with rasterio.drivers(CPL_DEBUG=verbosity > 2):
+
+        output, files = resolve_inout(files=files, output=output)
+        input = files[0]
+
+        with rasterio.open(input) as src:
+            if bounds:
+                if disjoint_bounds(bounds, src.bounds):
+                    raise click.BadParameter('must overlap the extent of '
+                                             'the input raster',
+                                             param='--bounds',
+                                             param_hint='--bounds')
+            elif like:
+                with rasterio.open(like) as template_ds:
+                    bounds = template_ds.bounds
+                    if template_ds.crs != src.crs:
+                        bounds = transform_bounds(template_ds.crs, src.crs,
+                                                  *bounds)
+
+                    if disjoint_bounds(bounds, src.bounds):
+                        raise click.BadParameter('must overlap the extent of '
+                                                 'the input raster',
+                                                 param='--like',
+                                                 param_hint='--like')
+
+            else:
+                raise click.UsageError('--bounds or --like required')
+
+            window = src.window(*bounds)
+
+            out_kwargs = src.meta.copy()
+            out_kwargs.update({
+                'driver': driver,
+                'height': window[0][1] - window[0][0],
+                'width': window[1][1] - window[1][0],
+                'transform': src.window_transform(window)
+            })
+            out_kwargs.update(**creation_options)
+
+            with rasterio.open(output, 'w', **out_kwargs) as out:
+                out.write(src.read(window=window))
+
+
 @click.command(short_help="Copy and convert raster dataset.")
 @click.argument(
     'files',
diff --git a/rasterio/rio/features.py b/rasterio/rio/features.py
index 31d84c2..0345bbb 100644
--- a/rasterio/rio/features.py
+++ b/rasterio/rio/features.py
@@ -3,6 +3,7 @@ import logging
 from math import ceil
 import os
 import shutil
+import re
 
 import click
 import cligj
@@ -16,7 +17,7 @@ from .helpers import coords, resolve_inout, write_features, 
to_lower
 from . import options
 import rasterio
 from rasterio.transform import Affine
-
+from rasterio.coords import disjoint_bounds
 
 logger = logging.getLogger('rio')
 
@@ -121,10 +122,10 @@ def mask(
         bounds = geojson.get('bbox', calculate_bounds(geojson))
 
         with rasterio.open(input) as src:
-            disjoint_bounds = _disjoint_bounds(bounds, src.bounds)
+            has_disjoint_bounds = disjoint_bounds(bounds, src.bounds)
 
             if crop:
-                if disjoint_bounds:
+                if has_disjoint_bounds:
                     raise click.BadParameter('not allowed for GeoJSON outside '
                                              'the extent of the input raster',
                                              param=crop, param_hint='--crop')
@@ -134,7 +135,7 @@ def mask(
                 (r1, r2), (c1, c2) = window
                 mask_shape = (r2 - r1, c2 - c1)
             else:
-                if disjoint_bounds:
+                if has_disjoint_bounds:
                     click.echo('GeoJSON outside bounds of existing output '
                                'raster. Are they in different coordinate '
                                'reference systems?',
@@ -329,8 +330,11 @@ def shapes(
                 if not with_nodata:
                     kwargs['mask'] = msk
 
+                src_basename = os.path.basename(src.name)
+
                 # Yield GeoJSON features.
-                for g, i in rasterio.features.shapes(img, **kwargs):
+                for i, (g, val) in enumerate(
+                        rasterio.features.shapes(img, **kwargs)):
                     if projection == 'geographic':
                         g = rasterio.warp.transform_geom(
                             src.crs, 'EPSG:4326', g,
@@ -338,9 +342,9 @@ def shapes(
                     xs, ys = zip(*coords(g))
                     yield {
                         'type': 'Feature',
-                        'id': str(i),
+                        'id': "{0}:{1}".format(src_basename, i),
                         'properties': {
-                            'val': i, 'filename': os.path.basename(src.name)
+                            'val': val, 'filename': src_basename
                         },
                         'bbox': [min(xs), min(ys), max(xs), max(ys)],
                         'geometry': g
@@ -487,7 +491,7 @@ def rasterize(
                                              'existing output raster',
                                              param='input', param_hint='input')
 
-                if _disjoint_bounds(geojson_bounds, out.bounds):
+                if disjoint_bounds(geojson_bounds, out.bounds):
                     click.echo("GeoJSON outside bounds of existing output "
                                "raster. Are they in different coordinate "
                                "reference systems?",
@@ -521,7 +525,7 @@ def rasterize(
                                              '--like raster',
                                              param='input', param_hint='input')
 
-                if _disjoint_bounds(geojson_bounds, template_ds.bounds):
+                if disjoint_bounds(geojson_bounds, template_ds.bounds):
                     click.echo("GeoJSON outside bounds of --like raster. "
                                "Are they in different coordinate reference "
                                "systems?",
@@ -702,9 +706,3 @@ def bounds(ctx, input, precision, indent, compact, 
projection, dst_crs,
     except Exception:
         logger.exception("Exception caught during processing")
         raise click.Abort()
-
-
-def _disjoint_bounds(bounds1, bounds2):
-    return (bounds1[0] > bounds2[2] or bounds1[2] < bounds2[0] or
-            bounds1[1] > bounds2[3] or bounds1[3] < bounds2[1])
-
diff --git a/rasterio/rio/info.py b/rasterio/rio/info.py
index d86d014..71f5a37 100644
--- a/rasterio/rio/info.py
+++ b/rasterio/rio/info.py
@@ -278,7 +278,7 @@ def info(ctx, input, aspect, indent, namespace, 
meta_member, verbose, bidx,
     try:
         with rasterio.drivers(CPL_DEBUG=(verbosity > 2)):
             with rasterio.open(input, mode) as src:
-                info = src.meta
+                info = src.profile
                 info['transform'] = info['affine'][:6]
                 del info['affine']
                 info['shape'] = info['height'], info['width']
diff --git a/rasterio/rio/merge.py b/rasterio/rio/merge.py
index 0d15e97..d816ac7 100644
--- a/rasterio/rio/merge.py
+++ b/rasterio/rio/merge.py
@@ -1,6 +1,5 @@
 # Merge command.
 
-
 import logging
 import math
 import os.path
@@ -23,9 +22,17 @@ from rasterio.transform import Affine
 @options.resolution_opt
 @click.option('--nodata', type=float, default=None,
               help="Override nodata values defined in input datasets")
+@click.option('--force-overwrite', '-f', 'force_overwrite', is_flag=True,
+              type=bool, default=False,
+              help="Do not prompt for confirmation before overwriting output "
+                   "file")
+@click.option('--precision', type=int, default=7,
+              help="Number of decimal places of precision in alignment of "
+                   "pixels")
 @options.creation_options
 @click.pass_context
-def merge(ctx, files, output, driver, bounds, res, nodata, creation_options):
+def merge(ctx, files, output, driver, bounds, res, nodata, force_overwrite,
+        precision, creation_options):
     """Copy valid pixels from input files to an output file.
 
     All files must have the same number of bands, data type, and
@@ -43,148 +50,30 @@ def merge(ctx, files, output, driver, bounds, res, nodata, 
creation_options):
       --res 0.1 0.1  => --res 0.1 (square)
       --res 0.1 0.2  => --res 0.1 --res 0.2  (rectangular)
     """
-    import numpy as np
+
+    from rasterio.tools.merge import merge as merge_tool
 
     verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
     logger = logging.getLogger('rio')
 
-    if len(res) == 1:
-        res = (res[0], res[0])
-
-    try:
-        with rasterio.drivers(CPL_DEBUG=verbosity>2):
-            output, files = resolve_inout(files=files, output=output)
-
-            with rasterio.open(files[0]) as first:
-                first_res = first.res
-                kwargs = first.meta
-                kwargs.update(**creation_options)
-                kwargs.pop('affine')
-                nodataval = first.nodatavals[0]
-                dtype = first.dtypes[0]
-
-            if os.path.exists(output):
-                # TODO: prompt user to update existing file (-i option) like:
-                # overwrite b.tif? (y/n [n]) n
-                # not overwritten
-                dst = rasterio.open(output, 'r+')
-                nodataval = dst.nodatavals[0]
-                dtype = dst.dtypes[0]
-                dest = np.zeros((dst.count,) + dst.shape, dtype=dtype)
-            else:
-                # Create new output file.
-                # Extent from option or extent of all inputs.
-                if not bounds:
-                    # scan input files.
-                    xs = []
-                    ys = []
-                    for f in files:
-                        with rasterio.open(f) as src:
-                            left, bottom, right, top = src.bounds
-                            xs.extend([left, right])
-                            ys.extend([bottom, top])
-                    bounds = min(xs), min(ys), max(xs), max(ys)
-                output_transform = Affine.translation(bounds[0], bounds[3])
-
-                # Resolution/pixel size.
-                if not res:
-                    res = first_res
-                output_transform *= Affine.scale(res[0], -res[1])
-
-                # Dataset shape.
-                output_width = int(math.ceil((bounds[2]-bounds[0])/res[0]))
-                output_height = int(math.ceil((bounds[3]-bounds[1])/res[1]))
-
-                kwargs['driver'] = driver
-                kwargs['transform'] = output_transform
-                kwargs['width'] = output_width
-                kwargs['height'] = output_height
-
-                logger.debug("Kwargs: %r", kwargs)
-                logger.debug("bounds: %r", bounds)
-                logger.debug("Res: %r", res)
-
-                dst = rasterio.open(output, 'w', **kwargs)
-                dest = np.zeros((first.count, output_height, output_width),
-                        dtype=dtype)
-
-                logger.debug("In merge, dest shape: %r", dest.shape)
-
-            if nodata is not None:
-                nodataval = nodata
-
-            if nodataval is not None:
-                # Only fill if the nodataval is within dtype's range.
-                inrange = False
-                if np.dtype(dtype).kind in ('i', 'u'):
-                    info = np.iinfo(dtype)
-                    inrange = (info.min <= nodataval <= info.max)
-                elif np.dtype(dtype).kind == 'f':
-                    info = np.finfo(dtype)
-                    inrange = (info.min <= nodataval <= info.max)
-                if inrange:
-                    dest.fill(nodataval)
-                else:
-                    warnings.warn(
-                        "Input file's nodata value, %s, is beyond the valid "
-                        "range of its data type, %s. Consider overriding it "
-                        "using the --nodata option for better results." % (
-                            nodataval, dtype))
-            else:
-                nodataval = 0
-
-            dst_w, dst_s, dst_e, dst_n = dst.bounds
-
-            for fname in reversed(files):
-                with rasterio.open(fname) as src:
-                    # Real World (tm) use of boundless reads.
-                    # This approach uses the maximum amount of memory to solve
-                    # the problem. Making it more efficient is a TODO.
-
-                    # 1. Compute spatial intersection of destination
-                    #    and source.
-                    src_w, src_s, src_e, src_n = src.bounds
-
-                    int_w = src_w if src_w > dst_w else dst_w
-                    int_s = src_s if src_s > dst_s else dst_s
-                    int_e = src_e if src_e < dst_e else dst_e
-                    int_n = src_n if src_n < dst_n else dst_n
-
-                    # 2. Compute the source window.
-                    src_window = src.window(int_w, int_s, int_e, int_n)
-
-                    # 3. Compute the destination window.
-                    dst_window = dst.window(int_w, int_s, int_e, int_n)
-
-                    # 4. Initialize temp array.
-                    temp = np.zeros(
-                            (first.count,) + tuple(b - a for a, b in 
dst_window),
-                            dtype=dtype)
-
-                    temp = src.read(
-                            out=temp,
-                            window=src_window,
-                            boundless=False,
-                            masked=True)
-
-                    # 5. Copy elements of temp into dest.
-                    roff, coff = dst.index(int_w, int_n)
-                    h, w = temp.shape[-2:]
-
-                    region = dest[:,roff:roff+h,coff:coff+w]
-                    np.copyto(region, temp,
-                        where=np.logical_and(
-                        region==nodataval, temp.mask==False))
-
-            if dst.mode == 'r+':
-                temp = dst.read(masked=True)
-                np.copyto(dest, temp,
-                    where=np.logical_and(
-                    dest==nodataval, temp.mask==False))
-
-            dst.write(dest)
-            dst.close()
-
-    except Exception:
-        logger.exception("Exception caught during processing")
-        raise click.Abort()
+    output, files = resolve_inout(files=files, output=output)
+
+    if os.path.exists(output) and not force_overwrite:
+        raise click.ClickException(
+            "Output exists and won't be overwritten without the "
+            "`-f` option")
+
+    sources = [rasterio.open(f) for f in files]
+    dest, output_transform = merge_tool(sources, bounds=bounds, res=res,
+                                        nodata=nodata, precision=precision)
+
+    profile = sources[0].profile
+    profile.pop('affine')
+    profile['transform'] = output_transform
+    profile['height'] = dest.shape[1]
+    profile['width'] = dest.shape[2]
+    profile['driver'] = driver
+    profile.update(**creation_options)
+
+    with rasterio.open(output, 'w', **profile) as dst:
+        dst.write(dest)
diff --git a/rasterio/tool.py b/rasterio/tool.py
index cb6364f..7804871 100644
--- a/rasterio/tool.py
+++ b/rasterio/tool.py
@@ -11,6 +11,7 @@ except ImportError:
 import numpy
 
 import rasterio
+from rasterio.five import zip_longest
 
 
 logger = logging.getLogger('rasterio')
@@ -48,6 +49,76 @@ def stats(source):
     return Stats(numpy.min(arr), numpy.max(arr), numpy.mean(arr))
 
 
+def show_hist(source, bins=10, masked=True, title='Histogram'):
+
+    """
+    Easily display a histogram with matplotlib.
+
+    Parameters
+    ----------
+    bins : int, optional
+        Compute histogram across N bins.
+    data : np.array or rasterio.Band or tuple(dataset, bidx)
+        Input data to display.  The first three arrays in multi-dimensional
+        arrays are plotted as red, green, and blue.
+    masked : bool, optional
+        When working with a `rasterio.Band()` object, specifies if the data
+        should be masked on read.
+    title : str, optional
+        Title for the figure.
+    """
+
+    if plt is None:
+        raise ImportError("Could not import matplotlib")
+
+    if isinstance(source, (tuple, rasterio.Band)):
+        arr = source[0].read(source[1], masked=masked)
+    else:
+        arr = source
+
+    # The histogram is computed individually for each 'band' in the array
+    # so we need the overall min/max to constrain the plot
+    rng = arr.min(), arr.max()
+
+    if len(arr.shape) is 2:
+        arr = [arr]
+        colors = ['gold']
+    else:
+        colors = ('red', 'green', 'blue', 'violet', 'gold', 'saddlebrown')
+
+    # If a rasterio.Band() is given make sure the proper index is displayed
+    # in the legend.
+    if isinstance(source, (tuple, rasterio.Band)):
+        labels = [str(source[1])]
+    else:
+        labels = (str(i + 1) for i in range(len(arr)))
+
+    # This loop should add a single plot each band in the input array,
+    # regardless of if the number of bands exceeds the number of colors.
+    # The colors slicing ensures that the number of iterations always
+    # matches the number of bands.
+    # The goal is to provide a curated set of colors for working with
+    # smaller datasets and let matplotlib define additional colors when
+    # working with larger datasets.
+    for bnd, color, label in zip_longest(arr, colors[:len(arr)], labels):
+
+        plt.hist(
+            bnd.flatten(),
+            bins=bins,
+            alpha=0.5,
+            color=color,
+            label=label,
+            range=rng
+        )
+
+    plt.legend(loc="upper right")
+    plt.title(title, fontweight='bold')
+    plt.grid(True)
+    plt.xlabel('DN')
+    plt.ylabel('Frequency')
+    plt.show()
+
+
 def main(banner, dataset, alt_interpreter=None):
     """ Main entry point for use with python interpreter """
     local = dict(funcs, src=dataset, np=numpy, rio=rasterio, plt=plt)
diff --git a/rasterio/tools/__init__.py b/rasterio/tools/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/rasterio/tools/merge.py b/rasterio/tools/merge.py
new file mode 100644
index 0000000..0fea083
--- /dev/null
+++ b/rasterio/tools/merge.py
@@ -0,0 +1,164 @@
+import logging
+import math
+import warnings
+
+import numpy as np
+
+import rasterio
+from rasterio._base import get_index, get_window
+from rasterio.transform import Affine
+
+
+logger = logging.getLogger('rasterio')
+
+
+def merge(sources, bounds=None, res=None, nodata=None, precision=7):
+    """Copy valid pixels from input files to an output file.
+
+    All files must have the same number of bands, data type, and
+    coordinate reference system.
+
+    Input files are merged in their listed order using the reverse
+    painter's algorithm. If the output file exists, its values will be
+    overwritten by input values.
+
+    Geospatial bounds and resolution of a new output file in the
+    units of the input file coordinate reference system may be provided
+    and are otherwise taken from the first input file.
+
+    Parameters
+    ----------
+    sources: list of source datasets 
+        Open rasterio RasterReader objects to be merged.
+    bounds: tuple, optional 
+        Bounds of the output image (left, bottom, right, top).
+        If not set, bounds are determined from bounds of input rasters. 
+    res: tuple, optional 
+        Output resolution in units of coordinate reference system. If not set,
+        the resolution of the first raster is used. If a single value is 
passed,
+        output pixels will be square.
+    nodata: float, optional
+        nodata value to use in output file. If not set, uses the nodata value
+        in the first input raster. 
+
+    Returns
+    -------
+    dest: numpy ndarray
+        Contents of all input rasters in single array.
+    out_transform: affine object
+        Information for mapping pixel coordinates in `dest` to another
+        coordinate system
+    """
+
+    first = sources[0]
+    first_res = first.res
+    nodataval = first.nodatavals[0]
+    dtype = first.dtypes[0]
+
+    # Extent from option or extent of all inputs.
+    if bounds:
+        dst_w, dst_s, dst_e, dst_n = bounds
+    else:
+        # scan input files.
+        xs = []
+        ys = []
+        for src in sources:
+           left, bottom, right, top = src.bounds
+           xs.extend([left, right])
+           ys.extend([bottom, top])
+        dst_w, dst_s, dst_e, dst_n = min(xs), min(ys), max(xs), max(ys)
+    
+    logger.debug("Output bounds: %r", (dst_w, dst_s, dst_e, dst_n))
+    output_transform = Affine.translation(dst_w, dst_n)
+    logger.debug("Output transform, before scaling: %r", output_transform)
+
+    # Resolution/pixel size.
+    if not res:
+        res = first_res
+    elif not np.iterable(res):
+        res = (res, res)
+    elif len(res) == 1:
+        res = (res[0], res[0])
+    output_transform *= Affine.scale(res[0], -res[1])
+    logger.debug("Output transform, after scaling: %r", output_transform)
+
+    # Compute output array shape. We guarantee it will cover the output
+    # bounds completely.
+    output_width = int(math.ceil((dst_e - dst_w) / res[0]))
+    output_height = int(math.ceil((dst_n - dst_s) / res[1]))
+
+    # Adjust bounds to fit.
+    dst_e, dst_s = output_transform * (output_width, output_height)
+    logger.debug("Output width: %d, height: %d", output_width, output_height)
+    logger.debug("Adjusted bounds: %r", (dst_w, dst_s, dst_e, dst_n))
+
+    # create destination array
+    dest = np.zeros((first.count, output_height, output_width), dtype=dtype)
+
+    if nodata is not None:
+        nodataval = nodata
+        logger.debug("Set nodataval: %r", nodataval)
+
+    if nodataval is not None:
+        # Only fill if the nodataval is within dtype's range.
+        inrange = False
+        if np.dtype(dtype).kind in ('i', 'u'):
+            info = np.iinfo(dtype)
+            inrange = (info.min <= nodataval <= info.max)
+        elif np.dtype(dtype).kind == 'f':
+            info = np.finfo(dtype)
+            inrange = (info.min <= nodataval <= info.max)
+        if inrange:
+            dest.fill(nodataval)
+        else:
+            warnings.warn(
+                "Input file's nodata value, %s, is beyond the valid "
+                "range of its data type, %s. Consider overriding it "
+                "using the --nodata option for better results." % (
+                    nodataval, dtype))
+    else:
+        nodataval = 0
+
+    for src in sources:
+        # Real World (tm) use of boundless reads.
+        # This approach uses the maximum amount of memory to solve the problem.
+        # Making it more efficient is a TODO.
+
+        # 1. Compute spatial intersection of destination and source.
+        src_w, src_s, src_e, src_n = src.bounds
+
+        int_w = src_w if src_w > dst_w else dst_w
+        int_s = src_s if src_s > dst_s else dst_s
+        int_e = src_e if src_e < dst_e else dst_e
+        int_n = src_n if src_n < dst_n else dst_n
+
+        # 2. Compute the source window.
+        src_window = get_window(
+            int_w, int_s, int_e, int_n, src.affine, precision=precision)
+        logger.debug("Src %s window: %r", src.name, src_window)
+
+        # 3. Compute the destination window.
+        dst_window = get_window(
+            int_w, int_s, int_e, int_n, output_transform, precision=precision)
+        logger.debug("Dst window: %r", dst_window)
+
+        # 4. Initialize temp array.
+        tcount = first.count
+        trows, tcols = tuple(b - a for a, b in dst_window)
+
+        temp_shape = (tcount, trows, tcols)
+        logger.debug("Temp shape: %r", temp_shape)
+
+        temp = np.zeros(temp_shape, dtype=dtype)
+        temp = src.read(out=temp, window=src_window, boundless=False,
+                        masked=True)
+
+        # 5. Copy elements of temp into dest.
+        roff, coff = dst_window[0][0], dst_window[1][0]
+
+        region = dest[:, roff:roff + trows, coff:coff + tcols]
+        np.copyto(
+            region, temp,
+            where=np.logical_and(region==nodataval, temp.mask==False))
+
+    return dest, output_transform
diff --git a/rasterio/transform.py b/rasterio/transform.py
index c02bc1c..b832711 100644
--- a/rasterio/transform.py
+++ b/rasterio/transform.py
@@ -2,6 +2,7 @@ import warnings
 
 from affine import Affine
 
+
 IDENTITY = Affine.identity()
 
 
@@ -38,3 +39,11 @@ def from_bounds(west, south, east, north, width, height):
     `height` in number of pixels."""
     return Affine.translation(west, north) * Affine.scale(
             (east - west)/width, (south - north)/height)
+
+
+def array_bounds(height, width, transform):
+    """Return the `west, south, east, north` bounds of an array given
+    its height, width, and an affine transform."""
+    w, n = transform.xoff, transform.yoff
+    e, s = transform * (width, height)
+    return w, s, e, n
diff --git a/requirements-dev.txt b/requirements-dev.txt
index a3e8d64..563cdc1 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,7 +1,7 @@
 affine
 cligj
 coveralls>=0.4
-cython==0.22
+cython>=0.23.1
 delocate
 enum34
 numpy>=1.8.0
diff --git a/setup.py b/setup.py
index c5b2c23..fe4cf4d 100755
--- a/setup.py
+++ b/setup.py
@@ -235,6 +235,7 @@ setup_args = dict(
         [rasterio.rio_commands]
         bounds=rasterio.rio.features:bounds
         calc=rasterio.rio.calc:calc
+        clip=rasterio.rio.convert:clip
         convert=rasterio.rio.convert:convert
         edit-info=rasterio.rio.info:edit
         env=rasterio.rio.info:env
diff --git a/tests/data/rgb1.tif b/tests/data/rgb1.tif
new file mode 100644
index 0000000..148b681
Binary files /dev/null and b/tests/data/rgb1.tif differ
diff --git a/tests/data/rgb2.tif b/tests/data/rgb2.tif
new file mode 100644
index 0000000..bc24905
Binary files /dev/null and b/tests/data/rgb2.tif differ
diff --git a/tests/data/rgb3.tif b/tests/data/rgb3.tif
new file mode 100644
index 0000000..ae89ad0
Binary files /dev/null and b/tests/data/rgb3.tif differ
diff --git a/tests/data/rgb4.tif b/tests/data/rgb4.tif
new file mode 100644
index 0000000..30b2c65
Binary files /dev/null and b/tests/data/rgb4.tif differ
diff --git a/tests/data/rgb_deflate.tif b/tests/data/rgb_deflate.tif
new file mode 100644
index 0000000..33eb8e5
Binary files /dev/null and b/tests/data/rgb_deflate.tif differ
diff --git a/tests/test_crs.py b/tests/test_crs.py
index 300b5cd..07b0e0c 100644
--- a/tests/test_crs.py
+++ b/tests/test_crs.py
@@ -89,6 +89,15 @@ def test_from_epsg_string():
         assert crs.from_string('epsg:xyz')
 
 
+def test_from_string():
+    wgs84_crs = crs.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 
+no_defs')
+    assert wgs84_crs == {'no_defs': True, 'ellps': 'WGS84', 'datum': 'WGS84', 
'proj': 'longlat'}
+
+    # Make sure this doesn't get handled using the from_epsg() even though 
'epsg' is in the string
+    epsg_init_crs = crs.from_string('+units=m +init=epsg:26911 +no_defs=True')
+    assert epsg_init_crs == {'units': 'm', 'init': 'epsg:26911', 'no_defs': 
True}
+
+
 def test_bare_parameters():
     """ Make sure that bare parameters (e.g., no_defs) are handled properly,
     even if they come in with key=True.  This covers interaction with pyproj,
diff --git a/tests/test_image_structure.py b/tests/test_image_structure.py
new file mode 100644
index 0000000..ed43973
--- /dev/null
+++ b/tests/test_image_structure.py
@@ -0,0 +1,74 @@
+import rasterio
+
+from rasterio.enums import Compression, Interleaving
+
+
+def test_enum_compression_JPEG():
+    assert Compression('JPEG').name == 'jpeg'
+
+
+def test_enum_compression_LZW():
+    assert Compression('LZW').name == 'lzw'
+
+
+def test_enum_compression_PACKBITS():
+    assert Compression('PACKBITS').name == 'packbits'
+
+
+def test_enum_compression_DEFLATE():
+    assert Compression('DEFLATE').name == 'deflate'
+
+
+def test_enum_compression_CCITTRLE():
+    assert Compression('CCITTRLE').name == 'ccittrle'
+
+
+def test_enum_compression_CCITTFAX3():
+    assert Compression('CCITTFAX3').name == 'ccittfax3'
+
+
+def test_enum_compression_CCITTFAX4():
+    assert Compression('CCITTFAX4').name == 'ccittfax4'
+
+
+def test_enum_compression_LZMA():
+    assert Compression('LZMA').name == 'lzma'
+
+
+def test_enum_compression_NONE():
+    assert Compression('NONE').name == 'none'
+
+
+def test_compression_none():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        assert src.compression is None
+        assert 'compress' not in src.profile
+
+
+def test_compression_deflate():
+    with rasterio.open('tests/data/rgb_deflate.tif') as src:
+        assert src.compression.name == 'deflate'
+        assert src.compression.value == 'DEFLATE'
+        assert src.profile['compress'] == 'deflate'
+
+
+def test_enum_interleaving_BAND():
+    assert Interleaving('BAND').name == 'band'
+
+
+def test_enum_interleaving_PIXEL():
+    assert Interleaving('PIXEL').name == 'pixel'
+
+
+def test_interleaving_pixel():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        assert src.interleaving.name == 'pixel'
+        assert src.interleaving.value == 'PIXEL'
+        assert src.profile['interleave'] == 'pixel'
+
+
+def test_interleaving_pixel():
+    with rasterio.open('tests/data/rgb_deflate.tif') as src:
+        assert src.interleaving.name == 'band'
+        assert src.interleaving.value == 'BAND'
+        assert src.profile['interleave'] == 'band'
diff --git a/tests/test_indexing.py b/tests/test_indexing.py
index c1c5f73..257e7dc 100644
--- a/tests/test_indexing.py
+++ b/tests/test_indexing.py
@@ -24,7 +24,7 @@ def test_window_no_exception():
                 (0, src.height), (-4, src.width))
 
 
-def test_index():
+def test_index_values():
     with rasterio.open('tests/data/RGB.byte.tif') as src:
         assert src.index(101985.0, 2826915.0) == (0, 0)
         assert src.index(101985.0+400.0, 2826915.0) == (0, 1)
@@ -41,11 +41,29 @@ def test_window():
                                                           (0, src.width))
         assert src.index(left+400, top-400) == (1, 1)
         assert src.index(left+dx+eps, top-dy-eps) == (1, 1)
-        assert src.window(left, top-400, left+400, top) == ((0, 1), (0, 1))
-        assert src.window(left, top-2*dy-eps, left+2*dx+eps, top) == ((0, 2), 
(0, 2))
+        assert src.window(left, top-400, left+400, top) == ((0, 2), (0, 2))
+        assert src.window(left, top-2*dy-eps, left+2*dx-eps, top) == ((0, 2), 
(0, 2))
 
 
 def test_window_bounds_roundtrip():
     with rasterio.open('tests/data/RGB.byte.tif') as src:
         assert ((100, 200), (100, 200)) == src.window(
             *src.window_bounds(((100, 200), (100, 200))))
+
+
+def test_window_full_cover():
+
+    def bound_covers(bounds1, bounds2):
+        """Does bounds1 cover bounds2?
+        """
+        return (bounds1[0] <= bounds2[0] and bounds1[1] <= bounds2[1] and
+                bounds1[2] >= bounds2[2] and bounds1[3] >= bounds2[3])
+
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        bounds = list(src.window_bounds(((100, 200), (100, 200))))
+        bounds[1] = bounds[1] - 10.0  # extend south
+        bounds[2] = bounds[2] + 10.0  # extend east
+
+        win = src.window(*bounds)
+        bounds_calc = list(src.window_bounds(win))
+        assert bound_covers(bounds_calc, bounds)
diff --git a/tests/test_read_boundless.py b/tests/test_read_boundless.py
index dc29c36..117479f 100644
--- a/tests/test_read_boundless.py
+++ b/tests/test_read_boundless.py
@@ -69,3 +69,18 @@ def test_read_boundless_masked_overlap():
         assert not data.mask.all()
         assert data.mask[0,399,399] == False
         assert data.mask[0,0,0] == True
+
+
+def test_read_boundless_zero_stop():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        data = src.read(
+            window=((-200, 0), (-200, 0)), boundless=True, masked=True)
+        assert data.shape == (3, 200, 200)
+        assert data.mask.all()
+
+
+def test_read_boundless_masks_zero_stop():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        data = src.read_masks(window=((-200, 0), (-200, 0)), boundless=True)
+        assert data.shape == (3, 200, 200)
+        assert data.min() == data.max() == src.nodata
diff --git a/tests/test_read_resample.py b/tests/test_read_resample.py
index d21c0a4..5593cf7 100644
--- a/tests/test_read_resample.py
+++ b/tests/test_read_resample.py
@@ -10,16 +10,17 @@ import rasterio
 
 def test_read_out_shape_resample_down():
     with rasterio.open('tests/data/RGB.byte.tif') as s:
-        out = numpy.zeros((7, 8), dtype=rasterio.ubyte)
+        out = numpy.zeros((8, 8), dtype=rasterio.ubyte)
         data = s.read(1, out=out)
         expected = numpy.array([
-            [  0,   8,   5,   7,   0,   0,   0,   0],
-            [  0,   6,  61,  15,  27,  15,  24, 128],
-            [  0,  20, 152,  23,  15,  19,  28,   0],
-            [  0,  17, 255,  25, 255,  22,  32,   0],
-            [  9,   7,  14,  16,  19,  18,  36,   0],
-            [  6,  27,  43, 207,  38,  31,  73,   0],
-            [  0,   0,   0,   0,  74,  23,   0,   0]], dtype=numpy.uint8)
+            [  0,   0,  20,  15,   0,   0,   0,   0],
+            [  0,   6, 193,   9, 255, 127,  23,  39],
+            [  0,   7,  27, 255, 193,  14,  28,  34],
+            [  0,  31,  29,  44,  14,  22,  43,   0],
+            [  0,   9,  69,  49,  17,  22, 255,   0],
+            [ 11,   7,  13,  25,  13,  29,  33,   0],
+            [  8,  10,  88,  27,  20,  33,  25,   0],
+            [  0,   0,   0,   0,  98,  23,   0,   0]], dtype=numpy.uint8)
         assert (data == expected).all() # all True.
 
 
diff --git a/tests/test_rio_convert.py b/tests/test_rio_convert.py
index 91590a2..e905180 100644
--- a/tests/test_rio_convert.py
+++ b/tests/test_rio_convert.py
@@ -1,15 +1,76 @@
 import sys
 import os
 import logging
-import click
+import numpy
 from click.testing import CliRunner
 
 import rasterio
-from rasterio.rio.convert import convert
+from rasterio.rio.main import main_group
+from rasterio.rio.convert import convert, clip
 
 
 logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
 
+TEST_BBOX = [-11850000, 4804000, -11840000, 4808000]
+
+
+def test_clip_bounds(runner, tmpdir):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        main_group,
+        ['clip', 'tests/data/shade.tif', output, '--bounds'] + TEST_BBOX
+    )
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+
+    with rasterio.open(output) as out:
+        assert out.shape == (420, 173)
+
+
+def test_clip_like(runner, tmpdir):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        clip,
+        ['tests/data/shade.tif', output, '--like', 'tests/data/shade.tif']
+    )
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+
+    with rasterio.open('tests/data/shade.tif') as template_ds:
+        with rasterio.open(output) as out:
+            assert out.shape == template_ds.shape
+            assert numpy.allclose(out.bounds, template_ds.bounds)
+
+
+def test_clip_missing_params(runner, tmpdir):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        clip,
+        ['tests/data/shade.tif', output]
+    )
+    assert result.exit_code == 2
+    assert '--bounds or --like required' in result.output
+
+
+def test_clip_bounds_disjunct(runner, tmpdir):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        clip,
+        ['tests/data/shade.tif', output, '--bounds'] + [0, 0, 10, 10]
+    )
+    assert result.exit_code == 2
+    assert '--bounds' in result.output
+
+
+def test_clip_like_disjunct(runner, tmpdir):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(
+        clip,
+        ['tests/data/shade.tif', output, '--like', 'tests/data/RGB.byte.tif']
+    )
+    assert result.exit_code == 2
+    assert '--like' in result.output
+
 
 # Tests: format and type conversion, --format and --dtype
 
diff --git a/tests/test_rio_features.py b/tests/test_rio_features.py
index dd3da8e..a1a4b15 100644
--- a/tests/test_rio_features.py
+++ b/tests/test_rio_features.py
@@ -213,7 +213,7 @@ def test_mask_crop(runner, tmpdir):
         with rasterio.open(output) as out:
             assert out.shape[1] == src.shape[1]
             assert out.shape[0] < src.shape[0]
-            assert out.shape[0] == 823
+            assert out.shape[0] == 824
 
     # Adding invert option after crop should be ignored
     result = runner.invoke(
@@ -310,10 +310,10 @@ def test_shapes_compact(runner):
 
 def test_shapes_sampling(runner):
     result = runner.invoke(
-        features.shapes, ['tests/data/shade.tif', '--sampling', '10'])
+        features.shapes, ['tests/data/shade.tif', '--sampling', '11'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
-    assert result.output.count('"Feature"') == 124
+    assert result.output.count('"Feature"') == 117
 
 
 def test_shapes_precision(runner):
diff --git a/tests/test_rio_merge.py b/tests/test_rio_merge.py
index d492373..db29e97 100644
--- a/tests/test_rio_merge.py
+++ b/tests/test_rio_merge.py
@@ -30,12 +30,12 @@ def test_data_dir_1(tmpdir):
 
     with rasterio.drivers():
 
-        with rasterio.open(str(tmpdir.join('a.tif')), 'w', **kwargs) as dst:
+        with rasterio.open(str(tmpdir.join('b.tif')), 'w', **kwargs) as dst:
             data = numpy.ones((10, 10), dtype=rasterio.uint8)
             data[0:6, 0:6] = 255
             dst.write_band(1, data)
 
-        with rasterio.open(str(tmpdir.join('b.tif')), 'w', **kwargs) as dst:
+        with rasterio.open(str(tmpdir.join('a.tif')), 'w', **kwargs) as dst:
             data = numpy.ones((10, 10), dtype=rasterio.uint8)
             data[4:8, 4:8] = 254
             dst.write_band(1, data)
@@ -58,12 +58,12 @@ def test_data_dir_2(tmpdir):
 
     with rasterio.drivers():
 
-        with rasterio.open(str(tmpdir.join('a.tif')), 'w', **kwargs) as dst:
+        with rasterio.open(str(tmpdir.join('b.tif')), 'w', **kwargs) as dst:
             data = numpy.zeros((10, 10), dtype=rasterio.uint8)
             data[0:6, 0:6] = 255
             dst.write_band(1, data)
 
-        with rasterio.open(str(tmpdir.join('b.tif')), 'w', **kwargs) as dst:
+        with rasterio.open(str(tmpdir.join('a.tif')), 'w', **kwargs) as dst:
             data = numpy.zeros((10, 10), dtype=rasterio.uint8)
             data[4:8, 4:8] = 254
             dst.write_band(1, data)
@@ -95,6 +95,7 @@ def test_merge_warn(test_data_dir_1):
     runner = CliRunner()
     result = runner.invoke(merge, inputs + [outputname] + ['--nodata', '-1'])
     assert result.exit_code == 0
+    assert os.path.exists(outputname)
     assert "using the --nodata option for better results" in result.output
 
 
@@ -130,12 +131,23 @@ def test_merge_output_exists(tmpdir):
         assert out.count == 3
 
 
-def test_merge_output_exists_without_nodata(test_data_dir_2):
+def test_merge_output_exists_without_nodata_fails(test_data_dir_2):
+    """Fails without -f or --force-overwrite"""
     runner = CliRunner()
     result = runner.invoke(
         merge,
         [str(test_data_dir_2.join('a.tif')),
             str(test_data_dir_2.join('b.tif'))])
+    assert result.exit_code == 1
+
+
+def test_merge_output_exists_without_nodata(test_data_dir_2):
+    """Succeeds with -f"""
+    runner = CliRunner()
+    result = runner.invoke(
+        merge,
+        ['-f', str(test_data_dir_2.join('a.tif')),
+            str(test_data_dir_2.join('b.tif'))])
     assert result.exit_code == 0
 
 
@@ -173,12 +185,12 @@ def test_data_dir_overlapping(tmpdir):
     }
 
     with rasterio.drivers():
-        with rasterio.open(str(tmpdir.join('nw.tif')), 'w', **kwargs) as dst:
+        with rasterio.open(str(tmpdir.join('se.tif')), 'w', **kwargs) as dst:
             data = numpy.ones((10, 10), dtype=rasterio.uint8)
             dst.write_band(1, data)
 
         kwargs['transform'] = (-113, 0.2, 0, 45, 0, -0.2)
-        with rasterio.open(str(tmpdir.join('se.tif')), 'w', **kwargs) as dst:
+        with rasterio.open(str(tmpdir.join('nw.tif')), 'w', **kwargs) as dst:
             data = numpy.ones((10, 10), dtype=rasterio.uint8) * 2
             dst.write_band(1, data)
 
@@ -219,12 +231,12 @@ def test_data_dir_float(tmpdir):
     }
 
     with rasterio.drivers():
-        with rasterio.open(str(tmpdir.join('one.tif')), 'w', **kwargs) as dst:
+        with rasterio.open(str(tmpdir.join('two.tif')), 'w', **kwargs) as dst:
             data = numpy.zeros((10, 10), dtype=rasterio.float64)
             data[0:6, 0:6] = 255
             dst.write_band(1, data)
 
-        with rasterio.open(str(tmpdir.join('two.tif')), 'w', **kwargs) as dst:
+        with rasterio.open(str(tmpdir.join('one.tif')), 'w', **kwargs) as dst:
             data = numpy.zeros((10, 10), dtype=rasterio.float64)
             data[4:8, 4:8] = 254
             dst.write_band(1, data)
@@ -295,15 +307,15 @@ def test_merge_tiny(tiffs):
 
     # Output should be
     #
-    # [[  0 120 120  90]
-    #  [  0 120 120  90]
+    # [[  0 120  90  90]
+    #  [  0 120  90  90]
     #  [  0  60   0   0]
     #  [ 40   0   0   0]]
 
     with rasterio.open(outputname) as src:
         data = src.read()
-        assert (data[0][0:2,1:3] == 120).all()
-        assert (data[0][0:2,3] == 90).all()
+        assert (data[0][0:2,1] == 120).all()
+        assert (data[0][0:2,2:4] == 90).all()
         assert data[0][2][1] == 60
         assert data[0][3][0] == 40
 
@@ -318,34 +330,72 @@ def test_merge_tiny_output_opt(tiffs):
 
     # Output should be
     #
-    # [[  0 120 120  90]
-    #  [  0 120 120  90]
+    # [[  0 120  90  90]
+    #  [  0 120  90  90]
     #  [  0  60   0   0]
     #  [ 40   0   0   0]]
 
     with rasterio.open(outputname) as src:
         data = src.read()
-        assert (data[0][0:2,1:3] == 120).all()
-        assert (data[0][0:2,3] == 90).all()
+        assert (data[0][0:2,1] == 120).all()
+        assert (data[0][0:2,2:4] == 90).all()
         assert data[0][2][1] == 60
         assert data[0][3][0] == 40
 
 
-def test_merge_tiny_res(tiffs):
+def test_merge_tiny_res_bounds(tiffs):
     outputname = str(tiffs.join('merged.tif'))
     inputs = [str(x) for x in tiffs.listdir()]
     inputs.sort()
     runner = CliRunner()
-    result = runner.invoke(merge, inputs + [outputname, '--res', 2])
+    result = runner.invoke(merge, inputs + [outputname, '--res', 2, 
'--bounds', 1, 0, 5, 4])
     assert result.exit_code == 0
 
     # Output should be
     # [[[120  90]
-    #   [  0   0]]]
+    #   [ 40   0]]]
 
     with rasterio.open(outputname) as src:
         data = src.read()
         print(data)
         assert data[0, 0, 0] == 120
         assert data[0, 0, 1] == 90
-        assert (data[0, 1, 0:1] == 0).all()
+        assert data[0, 1, 0] == 40
+        assert data[0, 1, 1] == 0
+
+
+def test_merge_tiny_res_high_precision(tiffs):
+    outputname = str(tiffs.join('merged.tif'))
+    inputs = [str(x) for x in tiffs.listdir()]
+    inputs.sort()
+    runner = CliRunner()
+    result = runner.invoke(merge, inputs + [outputname, '--res', 2, 
'--precision', 15])
+    assert result.exit_code == 0
+
+    # Output should be
+    # [[[120  90]
+    #   [ 40   0]]]
+
+    with rasterio.open(outputname) as src:
+        data = src.read()
+        print(data)
+        assert data[0, 0, 0] == 120
+        assert data[0, 0, 1] == 90
+        assert data[0, 1, 0] == 40
+        assert data[0, 1, 1] == 0
+
+
+def test_merge_rgb(tmpdir):
+    """Get back original image"""
+    outputname = str(tmpdir.join('merged.tif'))
+    inputs = [
+        'tests/data/rgb1.tif',
+        'tests/data/rgb2.tif',
+        'tests/data/rgb3.tif',
+        'tests/data/rgb4.tif']
+    runner = CliRunner()
+    result = runner.invoke(merge, inputs + [outputname])
+    assert result.exit_code == 0
+
+    with rasterio.open(outputname) as src:
+        assert [src.checksum(i) for i in src.indexes] == [25420, 29131, 37860]
diff --git a/tests/test_tool.py b/tests/test_tool.py
index f579c92..b0e5285 100644
--- a/tests/test_tool.py
+++ b/tests/test_tool.py
@@ -6,7 +6,7 @@ except ImportError:
     plt = None
 
 import rasterio
-from rasterio.tool import show, stats
+from rasterio.tool import show, show_hist, stats
 
 
 def test_stats():
@@ -42,3 +42,24 @@ def test_show():
                 show(src.read(1))
             except ImportError:
                 pass
+
+
+def test_show_hist():
+    """
+    This test only verifies that code up to the point of plotting with
+    matplotlib works correctly.  Tests do not exercise matplotlib.
+    """
+    if plt:
+        # Return because plotting causes the tests to block until the plot
+        # window is closed.
+        return
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        try:
+            show_hist((src, 1), bins=256)
+        except ImportError:
+            pass
+
+        try:
+            show_hist(src.read(), bins=256)
+        except ImportError:
+            pass

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

_______________________________________________
Pkg-grass-devel mailing list
Pkg-grass-devel@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel

Reply via email to