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

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

commit 498df05d2b5509e9d77d3ba48f99a19ff66843e4
Author: Johan Van de Wauw <johan.vandew...@gmail.com>
Date:   Tue Feb 10 21:37:23 2015 +0100

    Imported Upstream version 0.18
---
 CHANGES.txt                   |  13 +
 build-wheels.sh               |  18 -
 docs/cli.rst                  | 191 +++++++--
 docs/reproject.rst            |  50 ++-
 rasterio/__init__.py          |   2 +-
 rasterio/_base.pxd            |   1 +
 rasterio/_base.pyx            | 119 +++++-
 rasterio/_drivers.pyx         |   4 +-
 rasterio/_err.pyx             |   2 +-
 rasterio/_example.pyx         |   2 +
 rasterio/_features.pyx        |  64 ++-
 rasterio/_fill.pyx            |  84 ++++
 rasterio/_gdal.pxd            |   4 +-
 rasterio/_io.pyx              |  48 ++-
 rasterio/_warp.pyx            | 109 +-----
 rasterio/features.py          |  23 +-
 rasterio/fill.py              |  47 +++
 rasterio/rasterfill.cpp       | 886 ++++++++++++++++++++++++++++++++++++++++++
 rasterio/rio/features.py      | 236 ++++++++++-
 rasterio/rio/info.py          |  42 +-
 rasterio/rio/main.py          |   3 +-
 rasterio/rio/sample.py        |  94 +++++
 rasterio/warp.py              |   3 +-
 requirements-dev.txt          |   2 +-
 setup.py                      |  22 +-
 tests/conftest.py             |   8 +
 tests/test_features_bounds.py |  65 ++++
 tests/test_fillnodata.py      |  45 +++
 tests/test_rio_features.py    | 282 +++++++++++++-
 tests/test_rio_info.py        |  42 ++
 tests/test_rio_sample.py      |  81 ++++
 tests/test_sampling.py        |  17 +
 tests/test_tags.py            |   8 +
 tests/test_transform.py       |  14 +-
 34 files changed, 2398 insertions(+), 233 deletions(-)

diff --git a/CHANGES.txt b/CHANGES.txt
index fc51402..9dce1f6 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,6 +1,19 @@
 Changes
 =======
 
+0.18.0 (2015-02-10)
+-------------------
+- New rio-rasterize command (#187).
+- New window_transform method (#215).
+- New sample method and rio-sample command (#251, #275).
+- New fillnodata function based on GDAL's rasterfill.cpp (#253).
+- Speedups for _features and _warp modules (#259).
+- Enhancements for rio-info: 'res', 'lnglat', and 'stats' (#269, #270).
+
+0.17.1 (2015-01-20)
+-------------------
+- Properly handle metadata tags with values that contain "=" (#254).
+
 0.17.0 (2015-01-15)
 -------------------
 - Enhancements to rio-merge: relaxation of same-extent and same-resolution
diff --git a/build-wheels.sh b/build-wheels.sh
deleted file mode 100644
index 36e68dc..0000000
--- a/build-wheels.sh
+++ /dev/null
@@ -1,18 +0,0 @@
-#!/bin/bash
-
-# Automation of this is a TODO. For now, it depends on manually built libraries
-# as detailed in https://gist.github.com/sgillies/a8a2fb910a98a8566d0a.
-
-export MACOSX_DEPLOYMENT_TARGET=10.6
-export GDAL_CONFIG="/usr/local/bin/gdal-config"
-export PACKAGE_DATA=1
-
-VERSION=$1
-
-source $HOME/envs/riowhl27/bin/activate
-CFLAGS="`$GDAL_CONFIG --cflags`" LDFLAGS="`$GDAL_CONFIG --libs` `$GDAL_CONFIG 
--dep-libs`" python setup.py bdist_wheel -d wheels/$VERSION
-source $HOME/envs/riowhl34/bin/activate
-CFLAGS="`$GDAL_CONFIG --cflags`" LDFLAGS="`$GDAL_CONFIG --libs` `$GDAL_CONFIG 
--dep-libs`" python setup.py bdist_wheel -d wheels/$VERSION
-
-parallel delocate-wheel -w fixed_wheels/$VERSION --require-archs=intel -v {} 
::: wheels/$VERSION/*.whl
-parallel cp {} 
{.}.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl
 ::: fixed_wheels/$VERSION/*.whl
diff --git a/docs/cli.rst b/docs/cli.rst
index a744c4a..4e6adf2 100644
--- a/docs/cli.rst
+++ b/docs/cli.rst
@@ -13,19 +13,21 @@ Rasterio's new command line interface is a program named 
"rio".
     Options:
       -v, --verbose  Increase verbosity.
       -q, --quiet    Decrease verbosity.
+      --version      Show the version and exit.
       --help         Show this message and exit.
 
     Commands:
       bounds     Write bounding boxes to stdout as GeoJSON.
+      env        Print information about the rio environment.
       info       Print information about a data file.
       insp       Open a data file and start an interpreter.
       merge      Merge a stack of raster datasets.
+      rasterize  Rasterize features.
+      sample     Sample a dataset.
       shapes     Write the shapes of features.
       stack      Stack a number of bands into a multiband dataset.
-      transform  Transform coordinates.
-
-It is developed using the ``click`` package.
 
+It is developed using `Click <http://click.pocoo.org/3/>`__.
 
 bounds
 ------
@@ -84,26 +86,89 @@ Shoot the GeoJSON into a Leaflet map using geojsonio-cli by 
typing
 info
 ----
 
-Rio's info command intends to serve some of the same uses as gdalinfo.
+Rio's info command prints structured information about a dataset.
 
 .. code-block:: console
 
-    $ rio info tests/data/RGB.byte.tif
-    { 'affine': Affine(300.0379266750948, 0.0, 101985.0,
-           0.0, -300.041782729805, 2826915.0),
-      'count': 3,
-      'crs': { 'init': u'epsg:32618'},
-      'driver': u'GTiff',
-      'dtype': <type 'numpy.uint8'>,
-      'height': 718,
-      'nodata': 0.0,
-      'transform': ( 101985.0,
-                     300.0379266750948,
-                     0.0,
-                     2826915.0,
-                     0.0,
-                     -300.041782729805),
-      'width': 791}
+    $ rio info tests/data/RGB.byte.tif --indent 2
+    {
+      "count": 3,
+      "crs": "EPSG:32618",
+      "dtype": "uint8",
+      "driver": "GTiff",
+      "bounds": [
+        101985.0,
+        2611485.0,
+        339315.0,
+        2826915.0
+      ],
+      "lnglat": [
+        -77.75790625255473,
+        24.561583285327067
+      ],
+      "height": 718,
+      "width": 791,
+      "shape": [
+        718,
+        791
+      ],
+      "res": [
+        300.0379266750948,
+        300.041782729805
+      ],
+      "nodata": 0.0
+    }
+
+More information, such as band statistics, can be had using the `--verbose`
+option.
+
+.. code-block:: console
+
+    $ rio info tests/data/RGB.byte.tif --indent 2
+    {
+      "count": 3,
+      "crs": "EPSG:32618",
+      "stats": [
+        {
+          "max": 255.0,
+          "mean": 44.434478650699106,
+          "min": 1.0
+        },
+        {
+          "max": 255.0,
+          "mean": 66.02203484105824,
+          "min": 1.0
+        },
+        {
+          "max": 255.0,
+          "mean": 71.39316199120559,
+          "min": 1.0
+        }
+      ],
+      "dtype": "uint8",
+      "driver": "GTiff",
+      "bounds": [
+        101985.0,
+        2611485.0,
+        339315.0,
+        2826915.0
+      ],
+      "lnglat": [
+        -77.75790625255473,
+        24.561583285327067
+      ],
+      "height": 718,
+      "width": 791,
+      "shape": [
+        718,
+        791
+      ],
+      "res": [
+        300.0379266750948,
+        300.041782729805
+      ],
+      "nodata": 0.0
+    }
 
 insp
 ----
@@ -113,25 +178,12 @@ The insp command opens a dataset and an interpreter.
 .. code-block:: console
 
     $ rio insp tests/data/RGB.byte.tif
-    Rasterio 0.9 Interactive Inspector (Python 2.7.5)
+    Rasterio 0.18 Interactive Inspector (Python 2.7.9)
     Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
-    >>> import pprint
-    >>> pprint.pprint(src.meta)
-    {'affine': Affine(300.0379266750948, 0.0, 101985.0,
-           0.0, -300.041782729805, 2826915.0),
-     'count': 3,
-     'crs': {'init': u'epsg:32618'},
-     'driver': u'GTiff',
-     'dtype': <type 'numpy.uint8'>,
-     'height': 718,
-     'nodata': 0.0,
-     'transform': (101985.0,
-                   300.0379266750948,
-                   0.0,
-                   2826915.0,
-                   0.0,
-                   -300.041782729805),
-     'width': 791}
+    >>> print src.name
+    tests/data/RGB.byte.tif
+    >>> print src.bounds
+    BoundingBox(left=101985.0, bottom=2611485.0, right=339315.0, top=2826915.0)
 
 merge
 -----
@@ -143,6 +195,69 @@ datasets.
 
     $ rio merge rasterio/tests/data/R*.tif merged.tif
 
+rasterize
+---------
+
+New in 0.18.
+
+The rasterize command rasterizes GeoJSON features into a new or existing
+raster.
+
+.. code-block:: console
+
+    $ rio rasterize test.tif --res 0.0167 < input.geojson
+
+The resulting file will have an upper left coordinate determined by the bounds
+of the GeoJSON (in EPSG:4326, which is the default), with a
+pixel size of approximately 30 arc seconds.  Pixels whose center is within the
+polygon or that are selected by brezenhams line algorithm will be burned in
+with a default value of 1.
+
+It is possible to rasterize into an existing raster and use an alternative
+default value:
+
+.. code-block:: console
+
+    $ rio rasterize existing.tif --default_value 10 < input.geojson
+
+It is also possible to rasterize using a template raster, which will be used
+to determine the transform, dimensions, and coordinate reference system of the
+output raster:
+
+.. code-block:: console
+
+    $ rio rasterize test.tif --like tests/data/shade.tif < input.geojson
+
+GeoJSON features may be provided using stdin or specified directly as first
+argument, and dimensions may be provided in place of pixel resolution:
+
+.. code-block:: console
+
+    $ rio rasterize input.geojson test.tif --dimensions 1024 1024
+
+Other options are available, see:
+
+.. code-block:: console
+
+    $ rio rasterize --help
+
+sample
+------
+
+New in 0.18.
+
+The sample command reads ``x, y`` positions from stdin and writes the dataset
+values at that position to stdout.
+
+.. code-block:: console
+
+    $ cat << EOF | rio sample tests/data/RGB.byte.tif
+    > [220649.99999832606, 2719199.999999095]
+    > EOF
+    [18, 25, 14]
+
+The output of the transform command (see below) makes good input for sample.
+
 shapes
 ------
 
diff --git a/docs/reproject.rst b/docs/reproject.rst
index ad290b4..0b467f4 100644
--- a/docs/reproject.rst
+++ b/docs/reproject.rst
@@ -53,9 +53,57 @@ transform.
         assert destination.any()
         assert not destination.all()
 
+
 See `examples/reproject.py 
<https://github.com/mapbox/rasterio/blob/master/examples/reproject.py>`__ for 
code that writes the destination array to a GeoTIFF file. I've 
 uploaded the resulting file to a Mapbox map to demonstrate that the 
reprojection is
-correct: 
https://a.tiles.mapbox.com/v3/sgillies.hfek2oko/page.html?secure=1#6/0.000/0.033
+correct: 
https://a.tiles.mapbox.com/v3/sgillies.hfek2oko/page.html?secure=1#6/0.000/0.033.
+
+Reprojecting a GeoTIFF dataset
+------------------------------
+
+Here's a more real-world example of using ``reproject()`` to make an output 
dataset zoomed out by a factor of 2.
+Methods of the ``rasterio.Affine`` class help us generate the output dataset's 
transform matrix and, thereby, its
+spatial extent. 
+
+.. code-block:: python
+
+    import numpy
+    import rasterio
+    from rasterio import Affine as A
+    from rasterio.warp import reproject, RESAMPLING
+
+    with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
+        src_transform = src.affine
+
+        # Zoom out by a factor of 2 from the center of the source
+        # dataset. The destination transform is the product of the
+        # source transform, a translation down and to the right, and
+        # a scaling.
+        dst_transform = src_transform*A.translation(
+            -src.width/2.0, -src.height/2.0)*A.scale(2.0)
+
+        data = src.read()
+
+        kwargs = src.meta
+        kwargs['transform'] = dst_transform
+
+        with rasterio.open('/tmp/zoomed-out.tif', 'w', **kwargs) as dst:
+
+            for i, band in enumerate(data, 1):
+                dest = numpy.zeros_like(band)
+
+                reproject(
+                    band,
+                    dest,
+                    src_transform=src_transform,
+                    src_crs=src.crs,
+                    dst_transform=dst_transform,
+                    dst_crs=src.crs,
+                    resampling=RESAMPLING.nearest)
+
+                dst.write_band(i, dest)
+
+.. image:: https://farm8.staticflickr.com/7399/16390100651_54f01b8601_b_d.jpg)
 
 References
 ----------
diff --git a/rasterio/__init__.py b/rasterio/__init__.py
index 531dae2..0ea3dd2 100644
--- a/rasterio/__init__.py
+++ b/rasterio/__init__.py
@@ -18,7 +18,7 @@ from rasterio.transform import Affine, guard_transform
 
 __all__ = [
     'band', 'open', 'drivers', 'copy', 'pad']
-__version__ = "0.17"
+__version__ = "0.18"
 
 log = logging.getLogger('rasterio')
 class NullHandler(logging.Handler):
diff --git a/rasterio/_base.pxd b/rasterio/_base.pxd
index ad5440b..ca070b8 100644
--- a/rasterio/_base.pxd
+++ b/rasterio/_base.pxd
@@ -23,3 +23,4 @@ cdef class DatasetReader:
 
     cdef void *band(self, int bidx)
 
+cdef void *_osr_from_crs(object crs)
diff --git a/rasterio/_base.pyx b/rasterio/_base.pyx
index 660d112..96c7540 100644
--- a/rasterio/_base.pyx
+++ b/rasterio/_base.pyx
@@ -99,13 +99,13 @@ cdef class DatasetReader(object):
 
     def read_crs(self):
         cdef char *proj_c = NULL
-        cdef char *auth_key = NULL
-        cdef char *auth_val = NULL
+        cdef const char * auth_key = NULL
+        cdef const char * auth_val = NULL
         cdef void *osr = NULL
         if self._hds == NULL:
             raise ValueError("Null dataset")
         crs = {}
-        cdef char *wkt = _gdal.GDALGetProjectionRef(self._hds)
+        cdef const char * wkt = _gdal.GDALGetProjectionRef(self._hds)
         if wkt is NULL:
             raise ValueError("Unexpected NULL spatial reference")
         wkt_b = wkt
@@ -166,7 +166,7 @@ cdef class DatasetReader(object):
         cdef char *proj_c = NULL
         cdef char *key_c = NULL
         cdef void *osr = NULL
-        cdef char *wkt = NULL
+        cdef const char * wkt = NULL
         if self._hds == NULL:
             raise ValueError("Null dataset")
         wkt = _gdal.GDALGetProjectionRef(self._hds)
@@ -379,6 +379,11 @@ cdef class DatasetReader(object):
         lr = self.index(right, bottom)
         return tuple(zip(ul, lr))
 
+    def window_transform(self, window):
+        """Returns the affine transform for a dataset window."""
+        (r, _), (c, _) = window
+        return self.affine * Affine.translation(c or 0, r or 0)
+
     @property
     def meta(self):
         m = {
@@ -394,7 +399,14 @@ cdef class DatasetReader(object):
         self._read = True
         return m
 
-    
+    def lnglat(self):
+        w, s, e, n = self.bounds
+        cx = (w + e)/2.0
+        cy = (s + n)/2.0
+        lng, lat = _transform(
+                self.crs, {'init': 'epsg:4326'}, [cx], [cy], None)
+        return lng.pop(), lat.pop()
+
     def get_crs(self):
         # _read tells us that the CRS was read before and really is
         # None.
@@ -498,7 +510,7 @@ cdef class DatasetReader(object):
             item_c = papszStrList[i]
             item_b = item_c
             item = item_b.decode('utf-8')
-            key, value = item.split('=')
+            key, value = item.split('=', 1)
             retval[key] = value
         return retval
     
@@ -521,7 +533,7 @@ cdef class DatasetReader(object):
         cdef void *hBand
         cdef void *hTable
         cdef int i
-        cdef _gdal.GDALColorEntry *color
+        cdef const _gdal.GDALColorEntry * color
         if self._hds == NULL:
             raise ValueError("can't read closed raster file")
         if bidx not in self.indexes:
@@ -590,6 +602,7 @@ cpdef eval_window(object window, int height, int width):
             "invalid window: col range (%d, %d)" % (c_start, c_stop))
     return (r_start, r_stop), (c_start, c_stop)
 
+
 def window_shape(window, height=-1, width=-1):
     """Returns shape of a window.
 
@@ -599,8 +612,100 @@ def window_shape(window, height=-1, width=-1):
     (a, b), (c, d) = eval_window(window, height, width)
     return b-a, d-c
 
+
 def window_index(window):
     return tuple(slice(*w) for w in window)
 
+
 def tastes_like_gdal(t):
     return t[2] == t[4] == 0.0 and t[1] > 0 and t[5] < 0
+
+
+cdef void *_osr_from_crs(object crs):
+    cdef char *proj_c = NULL
+    cdef void *osr = _gdal.OSRNewSpatialReference(NULL)
+    params = []
+    # Normally, we expect a CRS dict.
+    if isinstance(crs, dict):
+        # EPSG is a special case.
+        init = crs.get('init')
+        if init:
+            auth, val = init.split(':')
+            if auth.upper() == 'EPSG':
+                _gdal.OSRImportFromEPSG(osr, int(val))
+        else:
+            crs['wktext'] = True
+            for k, v in crs.items():
+                if v is True or (k in ('no_defs', 'wktext') and v):
+                    params.append("+%s" % k)
+                else:
+                    params.append("+%s=%s" % (k, v))
+            proj = " ".join(params)
+            log.debug("PROJ.4 to be imported: %r", proj)
+            proj_b = proj.encode('utf-8')
+            proj_c = proj_b
+            _gdal.OSRImportFromProj4(osr, proj_c)
+    # Fall back for CRS strings like "EPSG:3857."
+    else:
+        proj_b = crs.encode('utf-8')
+        proj_c = proj_b
+        _gdal.OSRSetFromUserInput(osr, proj_c)
+    return osr
+
+
+def _transform(src_crs, dst_crs, xs, ys, zs):
+    cdef double *x = NULL
+    cdef double *y = NULL
+    cdef double *z = NULL
+    cdef char *proj_c = NULL
+    cdef void *src = NULL
+    cdef void *dst = NULL
+    cdef void *transform = NULL
+    cdef int i
+
+    assert len(xs) == len(ys)
+    assert zs is None or len(xs) == len(zs)
+
+    src = _osr_from_crs(src_crs)
+    dst = _osr_from_crs(dst_crs)
+
+    n = len(xs)
+    x = <double *>_gdal.CPLMalloc(n*sizeof(double))
+    y = <double *>_gdal.CPLMalloc(n*sizeof(double))
+    for i in range(n):
+        x[i] = xs[i]
+        y[i] = ys[i]
+
+    if zs is not None:
+        z = <double *>_gdal.CPLMalloc(n*sizeof(double))
+        for i in range(n):
+            z[i] = zs[i]
+
+    transform = _gdal.OCTNewCoordinateTransformation(src, dst)
+    res = _gdal.OCTTransform(transform, n, x, y, z)
+    #if res:
+    #    raise ValueError("Failed coordinate transformation")
+
+    res_xs = [0]*n
+    res_ys = [0]*n
+
+    for i in range(n):
+        res_xs[i] = x[i]
+        res_ys[i] = y[i]
+
+    if zs is not None:
+        res_zs = [0]*n
+        for i in range(n):
+            res_zs[i] = z[i]
+        _gdal.CPLFree(z)
+
+        retval = (res_xs, res_ys, res_zs)
+    else:
+        retval = (res_xs, res_ys)
+
+    _gdal.CPLFree(x)
+    _gdal.CPLFree(y)
+    _gdal.OCTDestroyCoordinateTransformation(transform)
+    _gdal.OSRDestroySpatialReference(src)
+    _gdal.OSRDestroySpatialReference(dst)
+    return retval
diff --git a/rasterio/_drivers.pyx b/rasterio/_drivers.pyx
index b8080da..40b8271 100644
--- a/rasterio/_drivers.pyx
+++ b/rasterio/_drivers.pyx
@@ -129,8 +129,8 @@ cdef class GDALEnv(object):
 
     def drivers(self):
         cdef void *drv = NULL
-        cdef char *key = NULL
-        cdef char *val = NULL
+        cdef const char *key = NULL
+        cdef const char *val = NULL
         cdef int i
         result = {}
         for i in range(GDALGetDriverCount()):
diff --git a/rasterio/_err.pyx b/rasterio/_err.pyx
index 44bb201..b59ce44 100644
--- a/rasterio/_err.pyx
+++ b/rasterio/_err.pyx
@@ -61,7 +61,7 @@ cdef class GDALErrCtxManager:
     def __exit__(self, exc_type=None, exc_val=None, exc_tb=None):
         cdef int err_type = CPLGetLastErrorType()
         cdef int err_no = CPLGetLastErrorNo()
-        cdef char *msg = CPLGetLastErrorMsg()
+        cdef const char *msg = CPLGetLastErrorMsg()
         # TODO: warn for err_type 2?
         if err_type >= 3:
             raise exception_map[err_no](msg)
diff --git a/rasterio/_example.pyx b/rasterio/_example.pyx
index c203847..e1b4396 100644
--- a/rasterio/_example.pyx
+++ b/rasterio/_example.pyx
@@ -1,3 +1,5 @@
+# cython: boundscheck=False
+
 import numpy
 cimport numpy
 
diff --git a/rasterio/_features.pyx b/rasterio/_features.pyx
index 3858af4..c7c63c6 100644
--- a/rasterio/_features.pyx
+++ b/rasterio/_features.pyx
@@ -1,4 +1,3 @@
-# cython: profile=True
 
 import logging
 import json
@@ -50,12 +49,12 @@ def _shapes(image, mask, connectivity, transform):
     """
 
     cdef int retval, rows, cols
-    cdef void *hband
-    cdef void *hmaskband
-    cdef void *hfdriver
-    cdef void *hfs
-    cdef void *hlayer
-    cdef void *fielddefn
+    cdef void *hband = NULL
+    cdef void *hmaskband = NULL
+    cdef void *hfdriver = NULL
+    cdef void *hfs = NULL
+    cdef void *hlayer = NULL
+    cdef void *fielddefn = NULL
     cdef _io.RasterReader rdr
     cdef _io.RasterReader mrdr
     cdef char **options = NULL
@@ -162,9 +161,9 @@ def _sieve(image, size, output, mask, connectivity):
     cdef InMemoryRaster in_mem_ds = None
     cdef InMemoryRaster out_mem_ds = None
     cdef InMemoryRaster mask_mem_ds = None
-    cdef void *in_band
-    cdef void *out_band
-    cdef void *mask_band
+    cdef void *in_band = NULL
+    cdef void *out_band = NULL
+    cdef void *mask_band = NULL
     cdef _io.RasterReader rdr
     cdef _io.RasterUpdater udr
     cdef _io.RasterReader mask_reader
@@ -291,6 +290,45 @@ def _rasterize(shapes, image, transform, all_touched):
             _gdal.CSLDestroy(options)
 
 
+def _explode(coords):
+    """Explode a GeoJSON geometry's coordinates object and yield
+    coordinate tuples. As long as the input is conforming, the type of
+    the geometry doesn't matter.  From Fiona 1.4.8"""
+    for e in coords:
+        if isinstance(e, (float, int)):
+            yield coords
+            break
+        else:
+            for f in _explode(e):
+                yield f
+
+
+def _bounds(geometry):
+    """Bounding box of a GeoJSON geometry.  
+    From Fiona 1.4.8 with updates here to handle feature collections.
+    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
+
+
 # Mapping of OGR integer geometry types to GeoJSON type names.
 GEOMETRY_TYPES = {
     0: 'Unknown',
@@ -416,12 +454,6 @@ cdef class GeomBuilder:
         return result
 
 
-cdef geometry(void *geom):
-    """Returns a GeoJSON object from an OGR geometry object."""
-
-    return GeomBuilder().build(geom)
-
-
 cdef class OGRGeomBuilder:
     """
     Builds an OGR geometry from GeoJSON geometry.
diff --git a/rasterio/_fill.pyx b/rasterio/_fill.pyx
new file mode 100644
index 0000000..8d7c998
--- /dev/null
+++ b/rasterio/_fill.pyx
@@ -0,0 +1,84 @@
+# distutils: language = c++
+# cython: profile=True
+#
+
+import numpy as np
+cimport numpy as np
+
+from rasterio import dtypes
+from rasterio._err import cpl_errs
+from rasterio cimport _gdal, _io
+
+from rasterio._io cimport InMemoryRaster
+
+
+def _fillnodata(image, mask, double max_search_distance=100.0,
+        int smoothing_iterations=0):
+    cdef void *memdriver = _gdal.GDALGetDriverByName("MEM")
+    cdef void *image_dataset
+    cdef void *image_band
+    cdef void *mask_dataset
+    cdef void *mask_band
+    cdef char **alg_options = NULL
+
+    if isinstance(image, np.ndarray):
+        # copy numpy ndarray into an in-memory dataset
+        image_dataset = _gdal.GDALCreate(
+            memdriver,
+            "image",
+            image.shape[1],
+            image.shape[0],
+            1,
+            <_gdal.GDALDataType>dtypes.dtype_rev[image.dtype.name],
+            NULL)
+        image_band = _gdal.GDALGetRasterBand(image_dataset, 1)
+        _io.io_auto(image, image_band, True)
+    elif isinstance(image, tuple):
+        # TODO
+        raise NotImplementedError()
+    else:
+        raise ValueError("Invalid source image")
+
+    if isinstance(mask, np.ndarray):
+        mask_cast = mask.astype('uint8')
+        mask_dataset = _gdal.GDALCreate(
+            memdriver,
+            "mask",
+            mask.shape[1],
+            mask.shape[0],
+            1,
+            <_gdal.GDALDataType>dtypes.dtype_rev['uint8'],
+            NULL)
+        mask_band = _gdal.GDALGetRasterBand(mask_dataset, 1)
+        _io.io_auto(mask_cast, mask_band, True)
+    elif isinstance(mask, tuple):
+        # TODO
+        raise NotImplementedError()
+    elif mask is None:
+        mask_band = NULL
+    else:
+        raise ValueError("Invalid source image mask")
+
+    with cpl_errs:
+        alg_options = _gdal.CSLSetNameValue(
+                alg_options, "TEMP_FILE_DRIVER", "MEM")
+
+        _gdal.GDALFillNodata(
+                image_band,
+                mask_band,
+                max_search_distance,
+                0,
+                smoothing_iterations,
+                alg_options,
+                NULL,
+                NULL)
+
+    # read the result into a numpy ndarray
+    result = np.empty(image.shape, dtype=image.dtype)
+    _io.io_auto(result, image_band, False)
+
+    _gdal.GDALClose(image_dataset)
+    _gdal.GDALClose(mask_dataset)
+    _gdal.CSLDestroy(alg_options)
+
+    return result
diff --git a/rasterio/_gdal.pxd b/rasterio/_gdal.pxd
index 1214678..0ae46f6 100644
--- a/rasterio/_gdal.pxd
+++ b/rasterio/_gdal.pxd
@@ -101,7 +101,7 @@ cdef extern from "gdal.h" nogil:
         short c3
         short c4
 
-    const GDALColorEntry *GDALGetColorEntry (void *hTable, int)
+    const GDALColorEntry * GDALGetColorEntry (void *hTable, int)
     void GDALSetColorEntry (void *hTable, int i, const GDALColorEntry *poEntry)
     int GDALSetRasterColorTable (void *hBand, void *hTable)
     void *GDALGetRasterColorTable (void *hBand)
@@ -211,5 +211,5 @@ cdef extern from "gdal_alg.h":
     int  GDALApproxTransform(void *pTransformArg, int bDstToSrc, int 
nPointCount, double *x, double *y, double *z, int *panSuccess )
     void GDALDestroyApproxTransformer( void * )
 
-
+    int GDALFillNodata(void *dst_band, void *mask_band, double 
max_search_distance, int deprecated, int smoothing_iterations, char **options, 
void *progress_func, void *progress_data)
 
diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx
index 98e18e1..e9339a5 100644
--- a/rasterio/_io.pyx
+++ b/rasterio/_io.pyx
@@ -559,8 +559,13 @@ cdef class RasterReader(_base.DatasetReader):
             a 2D array if it is a band index number.
 
         out: numpy ndarray, optional
-            An optional reference to an output array with the same
-            dimensions and shape.
+            As with Numpy ufuncs, this is an optional reference to an
+            output array with the same dimensions and shape into which
+            data will be placed.
+            
+            *Note*: the method's return value may be a view on this
+            array. In other words, `out` is likely to be an
+            incomplete representation of the method's results.
 
         window : a pair (tuple) of pairs of ints, optional
             The optional `window` argument is a 2 item tuple. The first
@@ -585,8 +590,11 @@ cdef class RasterReader(_base.DatasetReader):
 
         Returns
         -------
-        Numpy ndarray
+        Numpy ndarray or a view on a Numpy ndarray
 
+        Note: as with Numpy ufuncs, an object is returned even if you
+        use the optional `out` argument and the return value shall be
+        preferentially used by callers.
         """
 
         return2d = False
@@ -855,6 +863,29 @@ cdef class RasterReader(_base.DatasetReader):
             hmask, 0, xoff, yoff, width, height, out)
         return out
 
+    def sample(self, xy, indexes=None):
+        """Get the values of a dataset at certain positions
+
+        Parameters
+        ----------
+        xy : iterable, pairs of floats
+            A sequence or generator of (x, y) pairs.
+
+        indexes : list of ints or a single int, optional
+            If `indexes` is a list, the result is a 3D array, but is
+            a 2D array if it is a band index number.
+
+        Returns
+        -------
+        Iterable, yielding dataset values for the specified `indexes`
+        as an ndarray.
+        """
+        for x, y in xy:
+            r, c = self.index(x, y)
+            window = ((r, r+1), (c, c+1))
+            data = self.read(
+                    indexes, window=window, masked=False, boundless=True)
+            yield data[:,0,0]
 
 cdef class RasterUpdater(RasterReader):
     # Read-write access to raster data and metadata.
@@ -911,7 +942,8 @@ cdef class RasterUpdater(RasterReader):
     def start(self):
         cdef const char *drv_name = NULL
         cdef char **options = NULL
-        cdef char *key_c, *val_c = NULL
+        cdef char *key_c = NULL
+        cdef char *val_c = NULL
         cdef void *drv = NULL
         cdef void *hband = NULL
         cdef int success
@@ -1308,7 +1340,7 @@ cdef class RasterUpdater(RasterReader):
 
     def write_colormap(self, bidx, colormap):
         """Write a colormap for a band to the dataset."""
-        cdef void *hBand
+        cdef void *hBand = NULL
         cdef void *hTable
         cdef _gdal.GDALColorEntry color
         if self._hds == NULL:
@@ -1412,6 +1444,7 @@ cdef class InMemoryRaster:
         self.dataset = NULL
 
         cdef void *memdriver = _gdal.GDALGetDriverByName("MEM")
+        cdef int i = 0  # avoids Cython warning in for loop below
 
         # Several GDAL operations require the array of band IDs as input
         self.band_ids[0] = 1
@@ -1548,7 +1581,8 @@ cdef class IndirectRasterUpdater(RasterUpdater):
     def close(self):
         cdef const char *drv_name = NULL
         cdef char **options = NULL
-        cdef char *key_c, *val_c = NULL
+        cdef char *key_c = NULL
+        cdef char *val_c = NULL
         cdef void *drv = NULL
         cdef void *temp = NULL
         cdef int success
@@ -1598,7 +1632,7 @@ def writer(path, mode, **kwargs):
     # format driver's capabilities.
     cdef void *hds = NULL
     cdef void *drv = NULL
-    cdef char *drv_name = NULL
+    cdef const char *drv_name = NULL
     cdef const char *fname = NULL
 
     if mode == 'w' and 'driver' in kwargs:
diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx
index d9faaed..b84470a 100644
--- a/rasterio/_warp.pyx
+++ b/rasterio/_warp.pyx
@@ -1,6 +1,4 @@
 # distutils: language = c++
-# cython: profile=True
-#
 
 from collections import namedtuple
 import logging
@@ -8,7 +6,7 @@ import logging
 import numpy as np
 cimport numpy as np
 
-from rasterio cimport _gdal, _ogr, _io, _features
+from rasterio cimport _base, _gdal, _ogr, _io, _features
 from rasterio import dtypes
 
 
@@ -78,94 +76,6 @@ def tastes_like_gdal(t):
     return t[2] == t[4] == 0.0 and t[1] > 0 and t[5] < 0
 
 
-cdef void *_osr_from_crs(object crs):
-    cdef char *proj_c = NULL
-    cdef void *osr
-    osr = _gdal.OSRNewSpatialReference(NULL)
-    params = []
-    # Normally, we expect a CRS dict.
-    if isinstance(crs, dict):
-        # EPSG is a special case.
-        init = crs.get('init')
-        if init:
-            auth, val = init.split(':')
-            if auth.upper() == 'EPSG':
-                _gdal.OSRImportFromEPSG(osr, int(val))
-        else:
-            crs['wktext'] = True
-            for k, v in crs.items():
-                if v is True or (k in ('no_defs', 'wktext') and v):
-                    params.append("+%s" % k)
-                else:
-                    params.append("+%s=%s" % (k, v))
-            proj = " ".join(params)
-            log.debug("PROJ.4 to be imported: %r", proj)
-            proj_b = proj.encode('utf-8')
-            proj_c = proj_b
-            _gdal.OSRImportFromProj4(osr, proj_c)
-    # Fall back for CRS strings like "EPSG:3857."
-    else:
-        proj_b = crs.encode('utf-8')
-        proj_c = proj_b
-        _gdal.OSRSetFromUserInput(osr, proj_c)
-    return osr
-
-
-def _transform(src_crs, dst_crs, xs, ys, zs):
-    cdef double *x, *y, *z = NULL
-    cdef char *proj_c = NULL
-    cdef void *src, *dst
-    cdef void *transform
-    cdef int i
-
-    assert len(xs) == len(ys)
-    assert zs is None or len(xs) == len(zs)
-
-    src = _osr_from_crs(src_crs)
-    dst = _osr_from_crs(dst_crs)
-
-    n = len(xs)
-    x = <double *>_gdal.CPLMalloc(n*sizeof(double))
-    y = <double *>_gdal.CPLMalloc(n*sizeof(double))
-    for i in range(n):
-        x[i] = xs[i]
-        y[i] = ys[i]
-
-    if zs is not None:
-        z = <double *>_gdal.CPLMalloc(n*sizeof(double))
-        for i in range(n):
-            z[i] = zs[i]
-
-    transform = _gdal.OCTNewCoordinateTransformation(src, dst)
-    res = _gdal.OCTTransform(transform, n, x, y, z)
-    #if res:
-    #    raise ValueError("Failed coordinate transformation")
-
-    res_xs = [0]*n
-    res_ys = [0]*n
-
-    for i in range(n):
-        res_xs[i] = x[i]
-        res_ys[i] = y[i]
-
-    if zs is not None:
-        res_zs = [0]*n
-        for i in range(n):
-            res_zs[i] = z[i]
-        _gdal.CPLFree(z)
-
-        retval = (res_xs, res_ys, res_zs)
-    else:
-        retval = (res_xs, res_ys)
-
-    _gdal.CPLFree(x)
-    _gdal.CPLFree(y)
-    _gdal.OCTDestroyCoordinateTransformation(transform)
-    _gdal.OSRDestroySpatialReference(src)
-    _gdal.OSRDestroySpatialReference(dst)
-    return retval
-
-
 def _transform_geom(
         src_crs, dst_crs, geom, antimeridian_cutting, antimeridian_offset,
         precision):
@@ -174,15 +84,16 @@ def _transform_geom(
     cdef char *key_c = NULL
     cdef char *val_c = NULL
     cdef char **options = NULL
-    cdef void *src, *dst
-    cdef void *transform
-    cdef OGRGeometryFactory *factory
-    cdef void *src_ogr_geom
-    cdef void *dst_ogr_geom
+    cdef void *src = NULL
+    cdef void *dst = NULL
+    cdef void *transform = NULL
+    cdef OGRGeometryFactory *factory = NULL
+    cdef void *src_ogr_geom = NULL
+    cdef void *dst_ogr_geom = NULL
     cdef int i
 
-    src = _osr_from_crs(src_crs)
-    dst = _osr_from_crs(dst_crs)
+    src = _base._osr_from_crs(src_crs)
+    dst = _base._osr_from_crs(dst_crs)
     transform = _gdal.OCTNewCoordinateTransformation(src, dst)
 
     # Transform options.
@@ -259,7 +170,7 @@ def _reproject(
     bands of datasets on disk, the coordinate reference systems and
     transforms will be read from the appropriate datasets.
     """
-    cdef int retval, rows, cols
+    cdef int retval=0, rows, cols
     cdef void *hrdriver
     cdef void *hdsin
     cdef void *hdsout
diff --git a/rasterio/features.py b/rasterio/features.py
index ea5884a..3673c3a 100644
--- a/rasterio/features.py
+++ b/rasterio/features.py
@@ -8,7 +8,7 @@ import warnings
 import numpy as np
 
 import rasterio
-from rasterio._features import _shapes, _sieve, _rasterize
+from rasterio._features import _shapes, _sieve, _rasterize, _bounds
 from rasterio.transform import IDENTITY, guard_transform
 from rasterio.dtypes import get_minimum_int_dtype
 
@@ -321,3 +321,24 @@ def rasterize(
         _rasterize(valid_shapes, out, transform.to_gdal(), all_touched)
 
     return out
+
+
+def bounds(geometry):
+    """Returns a (minx, miny, maxx, maxy) bounding box.  From Fiona 1.4.8.
+    Modified to return bbox from geometry if available.
+
+    Parameters
+    ----------
+    geometry: GeoJSON-like feature, feature collection, or geometry.
+
+    Returns
+    -------
+    tuple
+        Bounding box: (minx, miny, maxx, maxy)
+    """
+
+    if 'bbox' in geometry:
+        return tuple(geometry['bbox'])
+
+    geom = geometry.get('geometry') or geometry
+    return _bounds(geom)
diff --git a/rasterio/fill.py b/rasterio/fill.py
new file mode 100644
index 0000000..37ec58e
--- /dev/null
+++ b/rasterio/fill.py
@@ -0,0 +1,47 @@
+import rasterio
+from rasterio._fill import _fillnodata
+
+def fillnodata(
+        image,
+        mask=None,
+        max_search_distance=100.0,
+        smoothing_iterations=0):
+    """
+    Fill selected raster regions by interpolation from the edges.
+
+    This algorithm will interpolate values for all designated nodata pixels
+    (marked by zeros in `mask`). For each pixel a four direction conic search
+    is done to find values to interpolate from (using inverse distance
+    weighting). Once all values are interpolated, zero or more smoothing
+    iterations (3x3 average filters on interpolated pixels) are applied to
+    smooth out artifacts.
+    
+    This algorithm is generally suitable for interpolating missing regions of
+    fairly continuously varying rasters (such as elevation models for
+    instance). It is also suitable for filling small holes and cracks in more
+    irregularly varying images (like aerial photos). It is generally not so
+    great for interpolating a raster from sparse point data.
+    
+    Parameters
+    ----------
+    image : numpy ndarray
+        The source  containing nodata holes.
+    mask : numpy ndarray
+        A mask band indicating which pixels to interpolate. Pixels to
+        interpolate into are indicated by the value 0. Values of 1 indicate
+        areas to use during interpolation. Must be same shape as image.
+    max_search_distance : float, optional
+        The maxmimum number of pixels to search in all directions to find
+        values to interpolate from. The default is 100.
+    smoothing_iterations : integer, optional
+        The number of 3x3 smoothing filter passes to run. The default is 0.
+    
+    Returns
+    -------
+    out : numpy ndarray
+        The filled raster array.
+    """
+    max_search_distance = float(max_search_distance)
+    smoothing_iterations = int(smoothing_iterations)
+    with rasterio.drivers():
+        return _fillnodata(image, mask, max_search_distance, 
smoothing_iterations)
diff --git a/rasterio/rasterfill.cpp b/rasterio/rasterfill.cpp
new file mode 100644
index 0000000..df2297b
--- /dev/null
+++ b/rasterio/rasterfill.cpp
@@ -0,0 +1,886 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  GDAL
+ * Purpose:  Interpolate in nodata areas.
+ * Author:   Frank Warmerdam, warmer...@pobox.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2008, Frank Warmerdam
+ * Copyright (c) 2015, Sean Gillies <s...@mapbox.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ***************************************************************************/
+
+#include "gdal_alg.h"
+#include "cpl_conv.h"
+#include "cpl_string.h"
+
+CPL_CVSID("$Id$");
+
+/************************************************************************/
+/*                           GDALFilterLine()                           */
+/*                                                                      */
+/*      Apply 3x3 filtering one one scanline with masking for which     */
+/*      pixels are to be interpolated (ThisFMask) and which window      */
+/*      pixels are valid to include in the interpolation (TMask).       */
+/************************************************************************/
+
+static void
+GDALFilterLine( float *pafLastLine, float *pafThisLine, float *pafNextLine,
+                float *pafOutLine,
+                GByte *pabyLastTMask, GByte *pabyThisTMask, 
GByte*pabyNextTMask,
+                GByte *pabyThisFMask, int nXSize )
+
+{
+    int iX;
+
+    for( iX = 0; iX < nXSize; iX++ )
+    {
+        if( !pabyThisFMask[iX] )
+        {
+            pafOutLine[iX] = pafThisLine[iX];
+            continue;
+        }
+
+        CPLAssert( pabyThisTMask[iX] );
+
+        double dfValSum = 0.0;
+        double dfWeightSum = 0.0;
+
+        // Previous line
+        if( pafLastLine != NULL )
+        {
+            if( iX > 0 && pabyLastTMask[iX-1] )
+            {
+                dfValSum += pafLastLine[iX-1];
+                dfWeightSum += 1.0;
+            }
+            if( pabyLastTMask[iX] )
+            {
+                dfValSum += pafLastLine[iX];
+                dfWeightSum += 1.0;
+            }
+            if( iX < nXSize-1 && pabyLastTMask[iX+1] )
+            {
+                dfValSum += pafLastLine[iX+1];
+                dfWeightSum += 1.0;
+            }
+        }
+
+        // Current Line
+        if( iX > 0 && pabyThisTMask[iX-1] )
+        {
+            dfValSum += pafThisLine[iX-1];
+            dfWeightSum += 1.0;
+        }
+        if( pabyThisTMask[iX] )
+        {
+            dfValSum += pafThisLine[iX];
+            dfWeightSum += 1.0;
+        }
+        if( iX < nXSize-1 && pabyThisTMask[iX+1] )
+        {
+            dfValSum += pafThisLine[iX+1];
+            dfWeightSum += 1.0;
+        }
+
+        // Next line
+        if( pafNextLine != NULL )
+        {
+            if( iX > 0 && pabyNextTMask[iX-1] )
+            {
+                dfValSum += pafNextLine[iX-1];
+                dfWeightSum += 1.0;
+            }
+            if( pabyNextTMask[iX] )
+            {
+                dfValSum += pafNextLine[iX];
+                dfWeightSum += 1.0;
+            }
+            if( iX < nXSize-1 && pabyNextTMask[iX+1] )
+            {
+                dfValSum += pafNextLine[iX+1];
+                dfWeightSum += 1.0;
+            }
+        }
+
+        pafOutLine[iX] = (float) (dfValSum / dfWeightSum);
+    }
+}
+
+/************************************************************************/
+/*                          GDALMultiFilter()                           */
+/*                                                                      */
+/*      Apply multiple iterations of a 3x3 smoothing filter over a      */
+/*      band with masking controlling what pixels should be             */
+/*      filtered (FiltMaskBand non zero) and which pixels can be        */
+/*      considered valid contributors to the filter                     */
+/*      (TargetMaskBand non zero).                                      */
+/*                                                                      */
+/*      This implementation attempts to apply many iterations in        */
+/*      one IO pass by managing the filtering over a rolling buffer     */
+/*      of nIternations+2 scanlines.  While possibly clever this        */
+/*      makes the algorithm implementation largely                      */
+/*      incomprehensible.                                               */
+/************************************************************************/
+
+static CPLErr
+GDALMultiFilter( GDALRasterBandH hTargetBand, 
+                 GDALRasterBandH hTargetMaskBand, 
+                 GDALRasterBandH hFiltMaskBand,
+                 int nIterations,
+                 GDALProgressFunc pfnProgress, 
+                 void * pProgressArg )
+
+{
+    float *paf3PassLineBuf;
+    GByte *pabyTMaskBuf;
+    GByte *pabyFMaskBuf;
+    float *pafThisPass, *pafLastPass, *pafSLastPass;
+
+    int   nBufLines = nIterations + 2;
+    int   iPassCounter = 0;
+    int   nNewLine; // the line being loaded this time (zero based scanline)
+    int   nXSize = GDALGetRasterBandXSize( hTargetBand );
+    int   nYSize = GDALGetRasterBandYSize( hTargetBand );
+    CPLErr eErr = CE_None;
+
+/* -------------------------------------------------------------------- */
+/*      Report starting progress value.                                 */
+/* -------------------------------------------------------------------- */
+    if( !pfnProgress( 0.0, "Smoothing Filter...", pProgressArg ) )
+    {
+        CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+        return CE_Failure;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Allocate rotating buffers.                                      */
+/* -------------------------------------------------------------------- */
+    pabyTMaskBuf = (GByte *) VSIMalloc2(nXSize, nBufLines);
+    pabyFMaskBuf = (GByte *) VSIMalloc2(nXSize, nBufLines);
+
+    paf3PassLineBuf = (float *) VSIMalloc3(nXSize, nBufLines, 3 * 
sizeof(float));
+    if (pabyTMaskBuf == NULL || pabyFMaskBuf == NULL || paf3PassLineBuf == 
NULL)
+    {
+        CPLError(CE_Failure, CPLE_OutOfMemory,
+                 "Could not allocate enough memory for temporary buffers");
+        eErr = CE_Failure;
+        goto end;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Process rotating buffers.                                       */
+/* -------------------------------------------------------------------- */
+    for( nNewLine = 0; 
+         eErr == CE_None && nNewLine < nYSize+nIterations; 
+         nNewLine++ )
+    {
+/* -------------------------------------------------------------------- */
+/*      Rotate pass buffers.                                            */
+/* -------------------------------------------------------------------- */
+        iPassCounter = (iPassCounter + 1) % 3;
+
+        pafSLastPass = paf3PassLineBuf 
+            + ((iPassCounter+0)%3) * nXSize*nBufLines;
+        pafLastPass = paf3PassLineBuf 
+            + ((iPassCounter+1)%3) * nXSize*nBufLines;
+        pafThisPass = paf3PassLineBuf 
+            + ((iPassCounter+2)%3) * nXSize*nBufLines;
+
+/* -------------------------------------------------------------------- */
+/*      Where does the new line go in the rotating buffer?              */
+/* -------------------------------------------------------------------- */
+        int iBufOffset = nNewLine % nBufLines;
+
+/* -------------------------------------------------------------------- */
+/*      Read the new data line if it is't off the bottom of the         */
+/*      image.                                                          */
+/* -------------------------------------------------------------------- */
+        if( nNewLine < nYSize )
+        {
+            eErr = 
+                GDALRasterIO( hTargetMaskBand, GF_Read, 
+                              0, nNewLine, nXSize, 1, 
+                              pabyTMaskBuf + nXSize * iBufOffset, nXSize, 1, 
+                              GDT_Byte, 0, 0 );
+            
+            if( eErr != CE_None )
+                break;
+
+            eErr = 
+                GDALRasterIO( hFiltMaskBand, GF_Read, 
+                              0, nNewLine, nXSize, 1, 
+                              pabyFMaskBuf + nXSize * iBufOffset, nXSize, 1, 
+                              GDT_Byte, 0, 0 );
+            
+            if( eErr != CE_None )
+                break;
+
+            eErr = 
+                GDALRasterIO( hTargetBand, GF_Read, 
+                              0, nNewLine, nXSize, 1, 
+                              pafThisPass + nXSize * iBufOffset, nXSize, 1, 
+                              GDT_Float32, 0, 0 );
+            
+            if( eErr != CE_None )
+                break;
+        }
+
+/* -------------------------------------------------------------------- */
+/*      Loop over the loaded data, applying the filter to all loaded    */
+/*      lines with neighbours.                                          */
+/* -------------------------------------------------------------------- */
+        int iFLine;
+
+        for( iFLine = nNewLine-1;
+             eErr == CE_None && iFLine >= nNewLine-nIterations;
+             iFLine-- )
+        {
+            int iLastOffset, iThisOffset, iNextOffset;
+
+            iLastOffset = (iFLine-1) % nBufLines; 
+            iThisOffset = (iFLine  ) % nBufLines;
+            iNextOffset = (iFLine+1) % nBufLines;
+
+            // default to preserving the old value.
+            if( iFLine >= 0 )
+                memcpy( pafThisPass + iThisOffset * nXSize, 
+                        pafLastPass + iThisOffset * nXSize, 
+                        sizeof(float) * nXSize );
+
+            // currently this skips the first and last line.  Eventually 
+            // we will enable these too.  TODO
+            if( iFLine < 1 || iFLine >= nYSize-1 )
+            {
+                continue;
+            }
+
+            GDALFilterLine( 
+                pafSLastPass + iLastOffset * nXSize,
+                pafLastPass  + iThisOffset * nXSize, 
+                pafThisPass  + iNextOffset * nXSize, 
+                pafThisPass  + iThisOffset * nXSize,
+                pabyTMaskBuf + iLastOffset * nXSize,
+                pabyTMaskBuf + iThisOffset * nXSize,
+                pabyTMaskBuf + iNextOffset * nXSize,
+                pabyFMaskBuf + iThisOffset * nXSize, 
+                nXSize );
+        }
+
+/* -------------------------------------------------------------------- */
+/*      Write out the top data line that will be rolling out of our     */
+/*      buffer.                                                         */
+/* -------------------------------------------------------------------- */
+        int iLineToSave = nNewLine - nIterations;
+
+        if( iLineToSave >= 0 && eErr == CE_None )
+        {
+            iBufOffset = iLineToSave % nBufLines;
+
+            eErr = 
+                GDALRasterIO( hTargetBand, GF_Write, 
+                              0, iLineToSave, nXSize, 1, 
+                              pafThisPass + nXSize * iBufOffset, nXSize, 1, 
+                              GDT_Float32, 0, 0 );
+        }
+
+/* -------------------------------------------------------------------- */
+/*      Report progress.                                                */
+/* -------------------------------------------------------------------- */
+        if( eErr == CE_None
+            && !pfnProgress( (nNewLine+1) / (double) (nYSize+nIterations), 
+                             "Smoothing Filter...", pProgressArg ) )
+        {
+            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+            eErr = CE_Failure;
+        }
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Cleanup                                                         */
+/* -------------------------------------------------------------------- */
+end:
+    CPLFree( pabyTMaskBuf );
+    CPLFree( pabyFMaskBuf );
+    CPLFree( paf3PassLineBuf );
+
+    return eErr;
+}
+ 
+/************************************************************************/
+/*                             QUAD_CHECK()                             */
+/*                                                                      */
+/*      macro for checking whether a point is nearer than the           */
+/*      existing closest point.                                         */
+/************************************************************************/
+#define QUAD_CHECK(quad_dist, quad_value,                              \
+target_x, target_y, origin_x, origin_y, target_value )                 \
+                                                                       \
+if( quad_value != nNoDataVal )                                                 
\
+{                                                                      \
+    double dfDx = (double)target_x - (double)origin_x;                 \
+    double dfDy = (double)target_y - (double)origin_y;                 \
+    double dfDistSq = dfDx * dfDx + dfDy * dfDy;                       \
+                                                                       \
+    if( dfDistSq < quad_dist*quad_dist )                               \
+    {                                                                  \
+       CPLAssert( dfDistSq > 0.0 );                                    \
+        quad_dist = sqrt(dfDistSq);                                    \
+        quad_value = target_value;                                     \
+    }                                                                  \
+}
+
+/************************************************************************/
+/*                           GDALFillNodata()                           */
+/************************************************************************/
+
+/**
+ * Fill selected raster regions by interpolation from the edges.
+ *
+ * This algorithm will interpolate values for all designated 
+ * nodata pixels (marked by zeros in hMaskBand).  For each pixel
+ * a four direction conic search is done to find values to interpolate
+ * from (using inverse distance weighting).  Once all values are
+ * interpolated, zero or more smoothing iterations (3x3 average
+ * filters on interpolated pixels) are applied to smooth out 
+ * artifacts. 
+ *
+ * This algorithm is generally suitable for interpolating missing
+ * regions of fairly continuously varying rasters (such as elevation
+ * models for instance).  It is also suitable for filling small holes
+ * and cracks in more irregularly varying images (like airphotos).  It
+ * is generally not so great for interpolating a raster from sparse 
+ * point data - see the algorithms defined in gdal_grid.h for that case.
+ *
+ * @param hTargetBand the raster band to be modified in place. 
+ * @param hMaskBand a mask band indicating pixels to be interpolated (zero 
valued
+ * @param dfMaxSearchDist the maximum number of pixels to search in all 
+ * directions to find values to interpolate from.
+ * @param bDeprecatedOption unused argument, should be zero.
+ * @param nSmoothingIterations the number of 3x3 smoothing filter passes to 
+ * run (0 or more).
+ * @param papszOptions additional name=value options in a string list (the
+ * temporary file driver can be specified like TEMP_FILE_DRIVER=MEM).
+ * @param pfnProgress the progress function to report completion.
+ * @param pProgressArg callback data for progress function.
+ * 
+ * @return CE_None on success or CE_Failure if something goes wrong. 
+ */
+
+CPLErr CPL_STDCALL
+GDALFillNodata( GDALRasterBandH hTargetBand,
+                GDALRasterBandH hMaskBand,
+                double dfMaxSearchDist,
+                int bDeprecatedOption,
+                int nSmoothingIterations,
+                char **papszOptions,
+                GDALProgressFunc pfnProgress,
+                void * pProgressArg )
+
+{
+    VALIDATE_POINTER1( hTargetBand, "GDALFillNodata", CE_Failure );
+
+    int nXSize = GDALGetRasterBandXSize( hTargetBand );
+    int nYSize = GDALGetRasterBandYSize( hTargetBand );
+    CPLErr eErr = CE_None;
+
+    // Special "x" pixel values identifying pixels as special.
+    GUInt32 nNoDataVal;
+    GDALDataType eType;
+
+    if( dfMaxSearchDist == 0.0 )
+        dfMaxSearchDist = MAX(nXSize,nYSize) + 1;
+
+    int nMaxSearchDist = (int) floor(dfMaxSearchDist);
+
+    if( nXSize > 65533 || nYSize > 65533 )
+    {
+        eType = GDT_UInt32;
+        nNoDataVal = 4000002;
+    }
+    else
+    {
+        eType = GDT_UInt16;
+        nNoDataVal = 65535;
+    }
+
+    if( hMaskBand == NULL )
+        hMaskBand = GDALGetMaskBand( hTargetBand );
+
+    /* If there are smoothing iterations, reserve 10% of the progress for them 
*/
+    double dfProgressRatio = (nSmoothingIterations > 0) ? 0.9 : 1.0;
+
+/* -------------------------------------------------------------------- */
+/*      Initialize progress counter.                                    */
+/* -------------------------------------------------------------------- */
+    if( pfnProgress == NULL )
+        pfnProgress = GDALDummyProgress;
+
+    if( !pfnProgress( 0.0, "Filling...", pProgressArg ) )
+    {
+        CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+        return CE_Failure;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Determine format driver for temp work files.                    */
+/* -------------------------------------------------------------------- */
+    CPLString osTmpFileDriver = CSLFetchNameValueDef(
+            papszOptions, "TEMP_FILE_DRIVER", "MEM");
+    GDALDriverH hDriver = GDALGetDriverByName((const char *) osTmpFileDriver);
+
+    if (hDriver == NULL)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "Given driver is not registered");
+        return CE_Failure;
+    }
+
+    if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, NULL) == NULL)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "Given driver is incapable of creating temp work files");
+        return CE_Failure;
+    }
+
+    char **papszWorkFileOptions = NULL;
+    if (osTmpFileDriver == "GTiff") {
+        papszWorkFileOptions = CSLSetNameValue(
+                papszWorkFileOptions, "COMPRESS", "LZW");
+        papszWorkFileOptions = CSLSetNameValue(
+                papszWorkFileOptions, "BIGTIFF", "IF_SAFER");
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Create a work file to hold the Y "last value" indices.          */
+/* -------------------------------------------------------------------- */
+    GDALDatasetH hYDS;
+    GDALRasterBandH hYBand;
+
+    CPLString osTmpFile = CPLGenerateTempFilename("");
+    CPLString osYTmpFile = osTmpFile + "fill_y_work.tif";
+
+    hYDS = GDALCreate( hDriver, osYTmpFile, nXSize, nYSize, 1,
+                       eType, (char **) papszWorkFileOptions );
+
+    if ( hYDS == NULL )
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+            "Could not create Y index work file. Check driver capabilities.");
+        return CE_Failure;
+    }
+
+    hYBand = GDALGetRasterBand( hYDS, 1 );
+
+/* -------------------------------------------------------------------- */
+/*      Create a work file to hold the pixel value associated with      */
+/*      the "last xy value" pixel.                                      */
+/* -------------------------------------------------------------------- */
+    GDALDatasetH hValDS;
+    GDALRasterBandH hValBand;
+    CPLString osValTmpFile = osTmpFile + "fill_val_work.tif";
+
+    hValDS = GDALCreate( hDriver, osValTmpFile, nXSize, nYSize, 1,
+                         GDALGetRasterDataType( hTargetBand ), 
+                         (char **) papszWorkFileOptions );
+
+    if ( hValDS == NULL )
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+            "Could not create XY value work file. Check driver capabilities.");
+        return CE_Failure;
+    }
+
+    hValBand = GDALGetRasterBand( hValDS, 1 );
+
+/* -------------------------------------------------------------------- */
+/*      Create a mask file to make it clear what pixels can be filtered */
+/*      on the filtering pass.                                          */
+/* -------------------------------------------------------------------- */
+    GDALDatasetH hFiltMaskDS;
+    GDALRasterBandH hFiltMaskBand;
+    CPLString osFiltMaskTmpFile = osTmpFile + "fill_filtmask_work.tif";
+
+    hFiltMaskDS = 
+        GDALCreate( hDriver, osFiltMaskTmpFile, nXSize, nYSize, 1,
+                    GDT_Byte, (char **) papszWorkFileOptions );
+
+    if ( hFiltMaskDS == NULL )
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+            "Could not create mask work file. Check driver capabilities.");
+        return CE_Failure;
+    }
+
+    hFiltMaskBand = GDALGetRasterBand( hFiltMaskDS, 1 );
+
+/* -------------------------------------------------------------------- */
+/*      Allocate buffers for last scanline and this scanline.           */
+/* -------------------------------------------------------------------- */
+    GUInt32 *panLastY, *panThisY, *panTopDownY;
+    float   *pafLastValue, *pafThisValue, *pafScanline, *pafTopDownValue;
+    GByte   *pabyMask, *pabyFiltMask;
+    int     iX;
+    int     iY;
+
+    panLastY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
+    panThisY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
+    panTopDownY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
+    pafLastValue = (float *) VSICalloc(nXSize,sizeof(float));
+    pafThisValue = (float *) VSICalloc(nXSize,sizeof(float));
+    pafTopDownValue = (float *) VSICalloc(nXSize,sizeof(float));
+    pafScanline = (float *) VSICalloc(nXSize,sizeof(float));
+    pabyMask = (GByte *) VSICalloc(nXSize,1);
+    pabyFiltMask = (GByte *) VSICalloc(nXSize,1);
+    if (panLastY == NULL || panThisY == NULL || panTopDownY == NULL ||
+        pafLastValue == NULL || pafThisValue == NULL || pafTopDownValue == 
NULL ||
+        pafScanline == NULL || pabyMask == NULL || pabyFiltMask == NULL)
+    {
+        CPLError(CE_Failure, CPLE_OutOfMemory,
+                 "Could not allocate enough memory for temporary buffers");
+
+        eErr = CE_Failure;
+        goto end;
+    }
+
+    for( iX = 0; iX < nXSize; iX++ )
+    {
+        panLastY[iX] = nNoDataVal;
+    }
+
+/* ==================================================================== */
+/*      Make first pass from top to bottom collecting the "last         */
+/*      known value" for each column and writing it out to the work     */
+/*      files.                                                          */
+/* ==================================================================== */
+    
+    for( iY = 0; iY < nYSize && eErr == CE_None; iY++ )
+    {
+/* -------------------------------------------------------------------- */
+/*      Read data and mask for this line.                               */
+/* -------------------------------------------------------------------- */
+        eErr = 
+            GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1, 
+                          pabyMask, nXSize, 1, GDT_Byte, 0, 0 );
+
+        if( eErr != CE_None )
+            break;
+
+        eErr = 
+            GDALRasterIO( hTargetBand, GF_Read, 0, iY, nXSize, 1, 
+                          pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
+        
+        if( eErr != CE_None )
+            break;
+        
+/* -------------------------------------------------------------------- */
+/*      Figure out the most recent pixel for each column.               */
+/* -------------------------------------------------------------------- */
+        
+        for( iX = 0; iX < nXSize; iX++ )
+        {
+            if( pabyMask[iX] )
+            {
+                pafThisValue[iX] = pafScanline[iX];
+                panThisY[iX] = iY;
+            }
+            else if( iY <= dfMaxSearchDist + panLastY[iX] )
+            {
+                pafThisValue[iX] = pafLastValue[iX];
+                panThisY[iX] = panLastY[iX];
+            }
+            else
+            {
+                panThisY[iX] = nNoDataVal;
+            }
+        }
+        
+/* -------------------------------------------------------------------- */
+/*      Write out best index/value to working files.                    */
+/* -------------------------------------------------------------------- */
+        eErr = GDALRasterIO( hYBand, GF_Write, 0, iY, nXSize, 1, 
+                             panThisY, nXSize, 1, GDT_UInt32, 0, 0 );
+        if( eErr != CE_None )
+            break;
+
+        eErr = GDALRasterIO( hValBand, GF_Write, 0, iY, nXSize, 1, 
+                             pafThisValue, nXSize, 1, GDT_Float32, 0, 0 );
+        if( eErr != CE_None )
+            break;
+
+/* -------------------------------------------------------------------- */
+/*      Flip this/last buffers.                                         */
+/* -------------------------------------------------------------------- */
+        {
+            float *pafTmp = pafThisValue;
+            pafThisValue = pafLastValue;
+            pafLastValue = pafTmp;
+
+            GUInt32 *panTmp = panThisY;
+            panThisY = panLastY;
+            panLastY = panTmp;
+        }
+
+/* -------------------------------------------------------------------- */
+/*      report progress.                                                */
+/* -------------------------------------------------------------------- */
+        if( eErr == CE_None
+            && !pfnProgress( dfProgressRatio * (0.5*(iY+1) / (double)nYSize), 
+                             "Filling...", pProgressArg ) )
+        {
+            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+            eErr = CE_Failure;
+        }
+    }
+
+/* ==================================================================== */
+/*      Now we will do collect similar this/last information from       */
+/*      bottom to top and use it in combination with the top to         */
+/*      bottom search info to interpolate.                              */
+/* ==================================================================== */
+    for( iY = nYSize-1; iY >= 0 && eErr == CE_None; iY-- )
+    {
+        eErr = 
+            GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1, 
+                          pabyMask, nXSize, 1, GDT_Byte, 0, 0 );
+
+        if( eErr != CE_None )
+            break;
+
+        eErr = 
+            GDALRasterIO( hTargetBand, GF_Read, 0, iY, nXSize, 1, 
+                          pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
+        
+        if( eErr != CE_None )
+            break;
+        
+/* -------------------------------------------------------------------- */
+/*      Figure out the most recent pixel for each column.               */
+/* -------------------------------------------------------------------- */
+        
+        for( iX = 0; iX < nXSize; iX++ )
+        {
+            if( pabyMask[iX] )
+            {
+                pafThisValue[iX] = pafScanline[iX];
+                panThisY[iX] = iY;
+            }
+            else if( panLastY[iX] - iY <= dfMaxSearchDist )
+            {
+                pafThisValue[iX] = pafLastValue[iX];
+                panThisY[iX] = panLastY[iX];
+            }
+            else
+            {
+                panThisY[iX] = nNoDataVal;
+            }
+        }
+        
+/* -------------------------------------------------------------------- */
+/*      Load the last y and corresponding value from the top down pass. */
+/* -------------------------------------------------------------------- */
+        eErr = 
+            GDALRasterIO( hYBand, GF_Read, 0, iY, nXSize, 1, 
+                          panTopDownY, nXSize, 1, GDT_UInt32, 0, 0 );
+
+        if( eErr != CE_None )
+            break;
+
+        eErr = 
+            GDALRasterIO( hValBand, GF_Read, 0, iY, nXSize, 1, 
+                          pafTopDownValue, nXSize, 1, GDT_Float32, 0, 0 );
+
+        if( eErr != CE_None )
+            break;
+
+/* -------------------------------------------------------------------- */
+/*      Attempt to interpolate any pixels that are nodata.              */
+/* -------------------------------------------------------------------- */
+        memset( pabyFiltMask, 0, nXSize );
+        for( iX = 0; iX < nXSize; iX++ )
+        {
+            int iStep, iQuad;
+            int nThisMaxSearchDist = nMaxSearchDist;
+
+            // If this was a valid target - no change.
+            if( pabyMask[iX] )
+                continue;
+
+            // Quadrants 0:topleft, 1:bottomleft, 2:topright, 3:bottomright
+            double adfQuadDist[4];
+            double adfQuadValue[4];
+
+            for( iQuad = 0; iQuad < 4; iQuad++ )
+            {
+                adfQuadDist[iQuad] = dfMaxSearchDist + 1.0;
+                adfQuadValue[iQuad] = 0.0;
+            }
+            
+            // Step left and right by one pixel searching for the closest 
+            // target value for each quadrant. 
+            for( iStep = 0; iStep < nThisMaxSearchDist; iStep++ )
+            {
+                int iLeftX = MAX(0,iX - iStep);
+                int iRightX = MIN(nXSize-1,iX + iStep);
+                
+                // top left includes current line 
+                QUAD_CHECK(adfQuadDist[0],adfQuadValue[0], 
+                           iLeftX, panTopDownY[iLeftX], iX, iY,
+                           pafTopDownValue[iLeftX] );
+
+                // bottom left 
+                QUAD_CHECK(adfQuadDist[1],adfQuadValue[1], 
+                           iLeftX, panLastY[iLeftX], iX, iY, 
+                           pafLastValue[iLeftX] );
+
+                // top right and bottom right do no include center pixel.
+                if( iStep == 0 )
+                     continue;
+                    
+                // top right includes current line 
+                QUAD_CHECK(adfQuadDist[2],adfQuadValue[2], 
+                           iRightX, panTopDownY[iRightX], iX, iY,
+                           pafTopDownValue[iRightX] );
+
+                // bottom right
+                QUAD_CHECK(adfQuadDist[3],adfQuadValue[3], 
+                           iRightX, panLastY[iRightX], iX, iY,
+                           pafLastValue[iRightX] );
+
+                // every four steps, recompute maximum distance.
+                if( (iStep & 0x3) == 0 )
+                    nThisMaxSearchDist = (int) floor(
+                        MAX(MAX(adfQuadDist[0],adfQuadDist[1]),
+                            MAX(adfQuadDist[2],adfQuadDist[3])) );
+            }
+
+            double dfWeightSum = 0.0;
+            double dfValueSum = 0.0;
+            
+            for( iQuad = 0; iQuad < 4; iQuad++ )
+            {
+                if( adfQuadDist[iQuad] <= dfMaxSearchDist )
+                {
+                    double dfWeight = 1.0 / adfQuadDist[iQuad];
+
+                    dfWeightSum += dfWeight;
+                    dfValueSum += adfQuadValue[iQuad] * dfWeight;
+                }
+            }
+
+            if( dfWeightSum > 0.0 )
+            {
+                pabyMask[iX] = 255;
+                pabyFiltMask[iX] = 255;
+                pafScanline[iX] = (float) (dfValueSum / dfWeightSum);
+            }
+
+        }
+
+/* -------------------------------------------------------------------- */
+/*      Write out the updated data and mask information.                */
+/* -------------------------------------------------------------------- */
+        eErr = 
+            GDALRasterIO( hTargetBand, GF_Write, 0, iY, nXSize, 1, 
+                          pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
+        
+        if( eErr != CE_None )
+            break;
+
+        eErr = 
+            GDALRasterIO( hFiltMaskBand, GF_Write, 0, iY, nXSize, 1, 
+                          pabyFiltMask, nXSize, 1, GDT_Byte, 0, 0 );
+        
+        if( eErr != CE_None )
+            break;
+
+/* -------------------------------------------------------------------- */
+/*      Flip this/last buffers.                                         */
+/* -------------------------------------------------------------------- */
+        {
+            float *pafTmp = pafThisValue;
+            pafThisValue = pafLastValue;
+            pafLastValue = pafTmp;
+            
+            GUInt32 *panTmp = panThisY;
+            panThisY = panLastY;
+            panLastY = panTmp;
+        }
+
+/* -------------------------------------------------------------------- */
+/*      report progress.                                                */
+/* -------------------------------------------------------------------- */
+        if( eErr == CE_None
+            && !pfnProgress( dfProgressRatio*(0.5+0.5*(nYSize-iY) / 
(double)nYSize), 
+                             "Filling...", pProgressArg ) )
+        {
+            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+            eErr = CE_Failure;
+        }
+    }        
+
+/* ==================================================================== */
+/*      Now we will do iterative average filters over the               */
+/*      interpolated values to smooth things out and make linear        */
+/*      artifacts less obvious.                                         */
+/* ==================================================================== */
+    if( eErr == CE_None && nSmoothingIterations > 0 )
+    {
+        // force masks to be to flushed and recomputed.
+        GDALFlushRasterCache( hMaskBand );
+
+        void *pScaledProgress;
+        pScaledProgress =
+            GDALCreateScaledProgress( dfProgressRatio, 1.0, pfnProgress, NULL 
);
+
+        eErr = GDALMultiFilter( hTargetBand, hMaskBand, hFiltMaskBand, 
+                                nSmoothingIterations,
+                                GDALScaledProgress, pScaledProgress );
+
+        GDALDestroyScaledProgress( pScaledProgress );
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Close and clean up temporary files. Free working buffers        */
+/* -------------------------------------------------------------------- */
+end:
+    CPLFree(panLastY);
+    CPLFree(panThisY);
+    CPLFree(panTopDownY);
+    CPLFree(pafLastValue);
+    CPLFree(pafThisValue);
+    CPLFree(pafTopDownValue);
+    CPLFree(pafScanline);
+    CPLFree(pabyMask);
+    CPLFree(pabyFiltMask);
+
+    GDALClose( hYDS );
+    GDALClose( hValDS );
+    GDALClose( hFiltMaskDS );
+
+    CSLDestroy(papszWorkFileOptions);
+
+    GDALDeleteDataset( hDriver, osYTmpFile );
+    GDALDeleteDataset( hDriver, osValTmpFile );
+    GDALDeleteDataset( hDriver, osFiltMaskTmpFile );
+
+    return eErr;
+}
diff --git a/rasterio/rio/features.py b/rasterio/rio/features.py
index f6d4c24..ade9b31 100644
--- a/rasterio/rio/features.py
+++ b/rasterio/rio/features.py
@@ -1,6 +1,8 @@
 import functools
 import json
 import logging
+from math import ceil, floor
+import os
 import sys
 import warnings
 
@@ -8,13 +10,15 @@ import click
 from cligj import (
     precision_opt, indent_opt, compact_opt, projection_geographic_opt,
     projection_projected_opt, sequence_opt, use_rs_opt,
-    geojson_type_feature_opt, geojson_type_bbox_opt)
+    geojson_type_feature_opt, geojson_type_bbox_opt, files_inout_arg,
+    format_opt)
 
 import rasterio
 from rasterio.transform import Affine
 from rasterio.rio.cli import cli, coords, write_features
 
 
+logger = logging.getLogger('rio')
 warnings.simplefilter('default')
 
 
@@ -142,3 +146,233 @@ def shapes(
     except Exception:
         logger.exception("Failed. Exception caught")
         sys.exit(1)
+
+
+# Rasterize command.
+@cli.command(short_help='Rasterize features.')
+@files_inout_arg
+@format_opt
+@click.option('--like', type=click.Path(exists=True),
+              help='Raster dataset to use as a template for obtaining affine '
+              'transform (bounds and resolution), crs, data type, and driver '
+              'used to create the output.  Only a single band will be output.')
+@click.option('--bounds', nargs=4, type=float, default=None,
+              help='Output bounds: left, bottom, right, top.')
+@click.option('--dimensions', nargs=2, type=int, default=None,
+              help='Output dataset width, height in number of pixels.')
+@click.option('--res', multiple=True, type=float, default=None,
+              help='Output dataset resolution in units of coordinate '
+              'reference system. Pixels assumed to be square if this option is 
'
+              'used once, otherwise use: '
+              '--res pixel_width --res pixel_height')
+@click.option('--src_crs', default='EPSG:4326',
+              help='Source coordinate reference system.  Limited to EPSG '
+              'codes for now.  Used as output coordinate system if output does 
'
+              'not exist or --like option is not used. Default: EPSG:4326')
+@click.option('--all_touched', is_flag=True, default=False)
+@click.option('--default_value', type=float, default=1, help='Default value '
+              'for rasterized pixels')
+@click.option('--fill', type=float, default=0, help='Fill value for all pixels 
'
+              'not overlapping features.  Will be evaluated as NoData pixels '
+              'for output.  Default: 0')
+@click.option('--property', type=str, default=None, help='Property in '
+              'GeoJSON features to use for rasterized values.  Any features '
+              'that lack this property will be given --default_value instead.')
+@click.pass_context
+def rasterize(
+        ctx,
+        files,
+        driver,
+        like,
+        bounds,
+        dimensions,
+        res,
+        src_crs,
+        all_touched,
+        default_value,
+        fill,
+        property):
+
+    """Rasterize GeoJSON into a new or existing raster.
+
+    If the output raster exists, rio-rasterize will rasterize feature values
+    into all bands of that raster.  The GeoJSON is assumed to be in the same
+    coordinate reference system as the output.
+
+    --default_value or property values when using --property must be using a
+    data type valid for the data type of that raster.
+
+
+    If a template raster is provided using the --like option, the affine
+    transform, coordinate reference system, and data type from that raster will
+    be used to create the output.
+
+    --default_value or property values when using --property must be using a
+    data type valid for the data type of that raster.
+
+    --driver, --bounds, --dimensions, --res, and --src_crs are ignored when
+    output exists or --like raster is provided
+
+    The GeoJSON must be within the bounds of the existing output or --like
+    raster, or an error will be raised.
+
+
+    If the output does not exist and --like raster is not provided, the input
+    GeoJSON will be used to determine the bounds of the output unless
+    provided using --bounds.
+
+    --dimensions or --res are required in this case.
+
+    If --res is provided, the bottom and right coordinates of bounds are
+    ignored.
+
+
+    Note:
+    The GeoJSON is not projected to match the coordinate reference system
+    of the output or --like rasters at this time.  This functionality may be
+    added in the future.
+    """
+
+    import numpy as np
+    from rasterio.features import rasterize
+    from rasterio.features import bounds as calculate_bounds
+
+    verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+
+    files = list(files)
+    output = files.pop()
+    input = click.open_file(files.pop(0) if files else '-')
+
+    # If values are actually meant to be integers, we need to cast them
+    # as such or rasterize creates floating point outputs
+    if default_value == int(default_value):
+        default_value = int(default_value)
+    if fill == int(fill):
+        fill = int(fill)
+
+    with rasterio.drivers(CPL_DEBUG=verbosity > 2):
+
+        def feature_value(feature):
+            if property and 'properties' in feature:
+                return feature['properties'].get(property, default_value)
+            return default_value
+
+        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])
+
+        geojson = json.loads(input.read())
+        if 'features' in geojson:
+            geometries = []
+            for f in geojson['features']:
+                geometries.append((f['geometry'], feature_value(f)))
+        elif 'geometry' in geojson:
+            geometries = ((geojson['geometry'], feature_value(geojson)), )
+        else:
+            raise click.BadParameter('Invalid GeoJSON', param=input,
+                                     param_hint='input')
+
+        geojson_bounds = geojson.get('bbox', calculate_bounds(geojson))
+
+        if os.path.exists(output):
+            with rasterio.open(output, 'r+') as out:
+                if disjoint_bounds(geojson_bounds, out.bounds):
+                    raise click.BadParameter('GeoJSON outside bounds of '
+                                             'existing output raster',
+                                             param=input, param_hint='input')
+
+                meta = out.meta.copy()
+
+                result = rasterize(
+                    geometries,
+                    out_shape=(meta['height'], meta['width']),
+                    transform=meta.get('affine', meta['transform']),
+                    all_touched=all_touched,
+                    dtype=meta.get('dtype', None),
+                    default_value=default_value,
+                    fill = fill)
+
+                for bidx in range(1, meta['count'] + 1):
+                    data = out.read_band(bidx, masked=True)
+                    # Burn in any non-fill pixels, and update mask accordingly
+                    ne = result != fill
+                    data[ne] = result[ne]
+                    data.mask[ne] = False
+                    out.write_band(bidx, data)
+
+        else:
+            if like is not None:
+                template_ds = rasterio.open(like)
+                if disjoint_bounds(geojson_bounds, template_ds.bounds):
+                    raise click.BadParameter('GeoJSON outside bounds of '
+                                             '--like raster', param=input,
+                                             param_hint='input')
+
+                kwargs = template_ds.meta.copy()
+                kwargs['count'] = 1
+                template_ds.close()
+
+            else:
+                bounds = bounds or geojson_bounds
+
+                if src_crs == 'EPSG:4326':
+                    if (bounds[0] < -180 or bounds[2] > 180 or
+                        bounds[1] < -80 or bounds[3] > 80):
+                        raise click.BadParameter(
+                            'Bounds are beyond the valid extent for 
EPSG:4326.',
+                            ctx, param=bounds, param_hint='--bounds')
+
+                if dimensions:
+                    width, height = dimensions
+                    res = (
+                        (bounds[2] - bounds[0]) / float(width),
+                        (bounds[3] - bounds[1]) / float(height)
+                    )
+
+                else:
+                    if not res:
+                        raise click.BadParameter(
+                            'pixel dimensions are required',
+                            ctx, param=res, param_hint='--res')
+                    elif len(res) == 1:
+                        res = (res[0], res[0])
+
+                    width = max(int(ceil((bounds[2] - bounds[0]) /
+                                float(res[0]))), 1)
+                    height = max(int(ceil((bounds[3] - bounds[1]) /
+                                 float(res[1]))), 1)
+
+                src_crs = src_crs.upper()
+                if not src_crs.count('EPSG:'):
+                    raise click.BadParameter(
+                        'invalid CRS.  Must be an EPSG code.',
+                        ctx, param=src_crs, param_hint='--src_crs')
+
+                kwargs = {
+                    'count': 1,
+                    'crs': src_crs,
+                    'width': width,
+                    'height': height,
+                    'transform': Affine(res[0], 0, bounds[0], 0, -res[1],
+                                        bounds[3]),
+                    'driver': driver
+                }
+
+            result = rasterize(
+                geometries,
+                out_shape=(kwargs['height'], kwargs['width']),
+                transform=kwargs.get('affine', kwargs['transform']),
+                all_touched=all_touched,
+                dtype=kwargs.get('dtype', None),
+                default_value=default_value,
+                fill = fill)
+
+            if 'dtype' not in kwargs:
+                kwargs['dtype'] = result.dtype
+
+            kwargs['nodata'] = fill
+
+            with rasterio.open(output, 'w', **kwargs) as out:
+                out.write_band(1, result)
diff --git a/rasterio/rio/info.py b/rasterio/rio/info.py
index 9c34a9d..427f672 100644
--- a/rasterio/rio/info.py
+++ b/rasterio/rio/info.py
@@ -61,18 +61,29 @@ def env(ctx, key):
 @click.option('--bounds', 'meta_member', flag_value='bounds',
               help="Print the boundary coordinates "
                    "(left, bottom, right, top).")
+@click.option('--res', 'meta_member', flag_value='res',
+              help="Print pixel width and height.")
+@click.option('--lnglat', 'meta_member', flag_value='lnglat',
+              help="Print longitude and latitude at center.")
+@click.option('--stats', 'meta_member', flag_value='stats',
+              help="Print statistics (min, max, mean) of a single band "
+                   "(use --bidx).")
+@click.option('-v', '--tell-me-more', '--verbose', is_flag=True,
+              help="Output extra information.")
+@click.option('--bidx', type=int, default=1,
+              help="Input file band index (default: 1).")
 @click.pass_context
-def info(ctx, input, aspect, indent, namespace, meta_member):
+def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx):
     """Print metadata about the dataset as JSON.
 
     Optionally print a single metadata item as a string.
     """
     verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
     logger = logging.getLogger('rio')
-    stdout = click.get_text_stream('stdout')
+    mode = 'r' if (verbose or meta_member == 'stats') else 'r-'
     try:
         with rasterio.drivers(CPL_DEBUG=(verbosity > 2)):
-            with rasterio.open(input, 'r-') as src:
+            with rasterio.open(input, mode) as src:
                 info = src.meta
                 del info['affine']
                 del info['transform']
@@ -82,19 +93,30 @@ def info(ctx, input, aspect, indent, namespace, 
meta_member):
                 if proj4.startswith('+init=epsg'):
                     proj4 = proj4.split('=')[1].upper()
                 info['crs'] = proj4
+                info['res'] = src.res
+                info['lnglat'] = src.lnglat()
+                if verbose:
+                    stats = [{'min': float(b.min()),
+                              'max': float(b.max()),
+                              'mean': float(b.mean())} for b in src.read()]
+                    info['stats'] = stats
                 if aspect == 'meta':
-                    if meta_member:
+                    if meta_member == 'stats':
+                        band = src.read(bidx)
+                        click.echo('%f %f %f' % (
+                            float(band.min()),
+                            float(band.max()),
+                            float(band.mean())))
+                    elif meta_member:
                         if isinstance(info[meta_member], (list, tuple)):
-                            print(" ".join(map(str, info[meta_member])))
+                            click.echo(" ".join(map(str, info[meta_member])))
                         else:
-                            print(info[meta_member])
+                            click.echo(info[meta_member])
                     else:
-                        stdout.write(json.dumps(info, indent=indent))
-                        stdout.write("\n")
+                        click.echo(json.dumps(info, indent=indent))
                 elif aspect == 'tags':
-                    stdout.write(json.dumps(src.tags(ns=namespace), 
+                    click.echo(json.dumps(src.tags(ns=namespace), 
                                             indent=indent))
-                    stdout.write("\n")
         sys.exit(0)
     except Exception:
         logger.exception("Failed. Exception caught")
diff --git a/rasterio/rio/main.py b/rasterio/rio/main.py
index 1bccb93..bd17c85 100644
--- a/rasterio/rio/main.py
+++ b/rasterio/rio/main.py
@@ -2,7 +2,8 @@
 
 from rasterio.rio.cli import cli
 from rasterio.rio.bands import stack
-from rasterio.rio.features import shapes
+from rasterio.rio.features import shapes, rasterize
 from rasterio.rio.info import env, info
 from rasterio.rio.merge import merge
 from rasterio.rio.rio import bounds, insp, transform
+from rasterio.rio.sample import sample
diff --git a/rasterio/rio/sample.py b/rasterio/rio/sample.py
new file mode 100644
index 0000000..fd4684a
--- /dev/null
+++ b/rasterio/rio/sample.py
@@ -0,0 +1,94 @@
+import json
+import logging
+import sys
+import warnings
+
+import click
+
+import rasterio
+from rasterio.rio.cli import cli
+
+
+warnings.simplefilter('default')
+
+
+@cli.command(short_help="Sample a dataset.")
+@click.argument('files', nargs=-1, required=True, metavar='FILE "[x, y]"')
+@click.option('--bidx', default=None, help="Indexes of input file bands.")
+@click.pass_context
+def sample(ctx, files, bidx):
+    """Sample a dataset at one or more points
+
+    Sampling points (x, y) encoded as JSON arrays, in the coordinate
+    reference system of the dataset, are read from the second
+    positional argument or stdin. Values of the dataset's bands
+    are also encoded as JSON arrays and are written to stdout.
+
+    Example:
+
+    \b
+        $ cat << EOF | rio sample tests/data/RGB.byte.tif
+        > [220650, 2719200]
+        > [219650, 2718200]
+        > EOF
+        [28, 29, 27]
+        [25, 29, 19]
+
+    By default, rio-sample will sample all bands. Optionally, bands
+    may be specified using a simple syntax:
+
+      --bidx N samples the Nth band (first band is 1).
+
+      --bidx M,N,0 samples bands M, N, and O.
+
+      --bidx M..O samples bands M-O, inclusive.
+
+      --bidx ..N samples all bands up to and including N.
+
+      --bidx N.. samples all bands from N to the end.
+
+    Example:
+
+    \b
+        $ cat << EOF | rio sample tests/data/RGB.byte.tif --bidx ..2
+        > [220650, 2719200]
+        > [219650, 2718200]
+        > EOF
+        [28, 29]
+        [25, 29]
+
+    """
+    verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+    logger = logging.getLogger('rio')
+
+    files = list(files)
+    source = files.pop(0)
+    input = files.pop(0) if files else '-'
+
+    # Handle the case of file, stream, or string input.
+    try:
+        points = click.open_file(input).readlines()
+    except IOError:
+        points = [input]
+
+    try:
+        with rasterio.drivers(CPL_DEBUG=verbosity>2):
+            with rasterio.open(source) as src:
+                if bidx is None:
+                    indexes = src.indexes
+                elif '..' in bidx:
+                    start, stop = map(
+                        lambda x: int(x) if x else None, bidx.split('..'))
+                    if start is None:
+                        start = 1
+                    indexes = src.indexes[slice(start-1, stop)]
+                else:
+                    indexes = list(map(int, bidx.split(',')))
+                for vals in src.sample(
+                            (json.loads(line) for line in points),
+                            indexes=indexes):
+                    click.echo(json.dumps(vals.tolist()))
+        sys.exit(0)
+    except Exception:
+        logger.exception("Failed. Exception caught")
+        sys.exit(1)
diff --git a/rasterio/warp.py b/rasterio/warp.py
index 25fd5a7..dd971d8 100644
--- a/rasterio/warp.py
+++ b/rasterio/warp.py
@@ -1,6 +1,7 @@
 """Raster warping and reprojection"""
 
-from rasterio._warp import _reproject, _transform, _transform_geom, RESAMPLING
+from rasterio._base import _transform
+from rasterio._warp import _transform_geom, _reproject, RESAMPLING
 from rasterio.transform import guard_transform
 
 
diff --git a/requirements-dev.txt b/requirements-dev.txt
index fb826ea..4d9ff53 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -6,5 +6,5 @@ delocate
 enum34
 numpy>=1.8.0
 pytest
-setuptools
+setuptools>=0.9.8
 wheel
diff --git a/setup.py b/setup.py
index cd491cf..1788e52 100755
--- a/setup.py
+++ b/setup.py
@@ -11,16 +11,21 @@
 
 import logging
 import os
+import pprint
 import shutil
 import subprocess
 import sys
-from setuptools import setup
 
-from distutils.extension import Extension
+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
+
 # Parse the version from the fiona module.
 with open('rasterio/__init__.py') as f:
     for line in f:
@@ -87,10 +92,9 @@ try:
 except Exception as e:
     log.warning("Failed to get options via gdal-config: %s", str(e))
 
-# Conditionally copy PROJ.4 data. Presumes PROJ.4 is installed locally
-# with --prefix=/usr/local.
+# Conditionally copy PROJ.4 data.
 if os.environ.get('PACKAGE_DATA'):
-    projdatadir = '/usr/local/share/proj'
+    projdatadir = os.environ.get('PROJ_LIB', '/usr/local/share/proj')
     if os.path.exists(projdatadir):
         try:
             shutil.rmtree('rasterio/proj_data')
@@ -104,6 +108,8 @@ ext_options = dict(
     libraries=libraries,
     extra_link_args=extra_link_args)
 
+log.debug('ext_options:\n%s', pprint.pformat(ext_options))
+
 # When building from a repo, Cython is required.
 if os.path.exists("MANIFEST.in") and "clean" not in sys.argv:
     log.info("MANIFEST.in found, presume a repo, cythonizing...")
@@ -126,10 +132,12 @@ if os.path.exists("MANIFEST.in") and "clean" not in 
sys.argv:
         Extension(
             'rasterio._warp', ['rasterio/_warp.pyx'], **ext_options),
         Extension(
+            'rasterio._fill', ['rasterio/_fill.pyx', 
'rasterio/rasterfill.cpp'], **ext_options),
+        Extension(
             'rasterio._err', ['rasterio/_err.pyx'], **ext_options),
         Extension(
             'rasterio._example', ['rasterio/_example.pyx'], **ext_options),
-            ])
+        ], quiet=True)
 
 # If there's no manifest template, as in an sdist, we just specify .c files.
 else:
@@ -147,6 +155,8 @@ else:
         Extension(
             'rasterio._warp', ['rasterio/_warp.cpp'], **ext_options),
         Extension(
+            'rasterio._fill', ['rasterio/_fill.cpp', 
'rasterio/rasterfill.cpp'], **ext_options),
+        Extension(
             'rasterio._err', ['rasterio/_err.c'], **ext_options),
         Extension(
             'rasterio._example', ['rasterio/_example.c'], **ext_options),
diff --git a/tests/conftest.py b/tests/conftest.py
index e55a25e..c5a9933 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -4,6 +4,8 @@ import os
 import sys
 
 import pytest
+from click.testing import CliRunner
+
 
 if sys.version_info > (3,):
     reduce = functools.reduce
@@ -11,6 +13,7 @@ if sys.version_info > (3,):
 test_files = [os.path.join(os.path.dirname(__file__), p) for p in [
     'data/RGB.byte.tif', 'data/float.tif', 'data/float_nan.tif', 
'data/shade.tif']]
 
+
 def pytest_cmdline_main(config):
     # Bail if the test raster data is not present. Test data is not 
     # distributed with sdists since 0.12.
@@ -19,3 +22,8 @@ def pytest_cmdline_main(config):
     else:
         print("Test data not present. See download directions in 
tests/README.txt")
         sys.exit(1)
+
+
+@pytest.fixture(scope='function')
+def runner():
+    return CliRunner()
diff --git a/tests/test_features_bounds.py b/tests/test_features_bounds.py
new file mode 100644
index 0000000..a7015cd
--- /dev/null
+++ b/tests/test_features_bounds.py
@@ -0,0 +1,65 @@
+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_fillnodata.py b/tests/test_fillnodata.py
new file mode 100644
index 0000000..2ccee1e
--- /dev/null
+++ b/tests/test_fillnodata.py
@@ -0,0 +1,45 @@
+import logging
+import sys
+import numpy
+import pytest
+
+import rasterio
+from rasterio.fill import fillnodata
+
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+def test_fillnodata():
+    """Test filling nodata values in an ndarray"""
+    # create a 5x5 array, with some missing data
+    a = numpy.ones([3, 3]) * 42
+    a[1][1] = 0
+    # find the missing data
+    mask = ~(a == 0)
+    # fill the missing data using interpolation from the edges
+    result = fillnodata(a, mask)
+    assert(numpy.all((numpy.ones([3, 3]) * 42) == result))
+
+def test_fillnodata_invalid_types():
+    a = numpy.ones([3, 3])
+    with pytest.raises(ValueError):
+        fillnodata(None, a)
+    with pytest.raises(ValueError):
+        fillnodata(a, 42)
+
+def test_fillnodata_mask_ones():
+    # when mask is all ones, image should be unmodified
+    a = numpy.ones([3, 3]) * 42
+    a[1][1] = 0
+    mask = numpy.ones([3, 3])
+    result = fillnodata(a, mask)
+    assert(numpy.all(a == result))
+
+'''
+def test_fillnodata_smooth():
+    a = numpy.array([[1,3,3,1],[2,0,0,2],[2,0,0,2],[1,3,3,1]], 
dtype=numpy.float64)
+    mask = ~(a == 0)
+    result = fillnodata(a, mask, max_search_distance=1, smoothing_iterations=0)
+    assert(result[1][1] == 3)
+    result = fillnodata(a, mask, max_search_distance=1, smoothing_iterations=1)
+    assert(round(result[1][1], 1) == 2.2)
+'''
diff --git a/tests/test_rio_features.py b/tests/test_rio_features.py
index c886515..e524cd7 100644
--- a/tests/test_rio_features.py
+++ b/tests/test_rio_features.py
@@ -1,14 +1,81 @@
 import logging
+import os
 import re
 import sys
 
 import click
 from click.testing import CliRunner
 
+import rasterio
 from rasterio.rio import features
 
 
-logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+# 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"
+}"""
+
+
+# > 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_err():
@@ -18,24 +85,21 @@ def test_err():
     assert result.exit_code == 1
 
 
-def test_shapes():
-    runner = CliRunner()
+def test_shapes(runner):
     result = runner.invoke(features.shapes, ['tests/data/shade.tif'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
     assert result.output.count('"Feature"') == 232
 
 
-def test_shapes_sequence():
-    runner = CliRunner()
+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
 
 
-def test_shapes_sequence_rs():
-    runner = CliRunner()
+def test_shapes_sequence_rs(runner):
     result = runner.invoke(
         features.shapes, [
             'tests/data/shade.tif',
@@ -47,24 +111,21 @@ def test_shapes_sequence_rs():
     assert result.output.count(u'\u001e') == 232
 
 
-def test_shapes_with_nodata():
-    runner = CliRunner()
+def test_shapes_with_nodata(runner):
     result = runner.invoke(features.shapes, ['tests/data/shade.tif', 
'--with-nodata'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
     assert result.output.count('"Feature"') == 288
 
 
-def test_shapes_indent():
-    runner = CliRunner()
+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') == 70139
 
 
-def test_shapes_compact():
-    runner = CliRunner()
+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
@@ -72,8 +133,7 @@ def test_shapes_compact():
     assert result.output.count(': ') == 0
 
 
-def test_shapes_sampling():
-    runner = CliRunner()
+def test_shapes_sampling(runner):
     result = runner.invoke(
         features.shapes, ['tests/data/shade.tif', '--sampling', '10'])
     assert result.exit_code == 0
@@ -81,8 +141,7 @@ def test_shapes_sampling():
     assert result.output.count('"Feature"') == 124
 
 
-def test_shapes_precision():
-    runner = CliRunner()
+def test_shapes_precision(runner):
     result = runner.invoke(
         features.shapes, ['tests/data/shade.tif', '--precision', '1'])
     assert result.exit_code == 0
@@ -91,9 +150,194 @@ def test_shapes_precision():
     assert re.search(r'\d*\.\d{2,}', result.output) is None
 
 
-def test_shapes_mask():
-    runner = CliRunner()
+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"') == 9
+
+
+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
+
+    # Test invalid res
+    result = runner.invoke(features.rasterize, [output], input=TEST_FEATURES)
+    assert result.exit_code == 2
+
+    # Test invalid CRS for bounds
+    result = runner.invoke(features.rasterize, [output, '--res', 1],
+                           input=TEST_MERC_FEATURECOLLECTION)
+    assert result.exit_code == 2
+
+    # Test invalid CRS value
+    result = runner.invoke(features.rasterize, [output,
+                                                '--res', 1,
+                                                '--src_crs', 'BOGUS'],
+                           input=TEST_MERC_FEATURECOLLECTION)
+    assert result.exit_code == 2
+
+
+def test_rasterize(tmpdir, runner):
+    # Test dimensions
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, '--dimensions', 20, 10],
+                           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
+        data = out.read_band(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 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
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 748
+        assert (data == 1).sum() == 52
+
+    # 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
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 55
+        assert (data == 1).sum() == 145
+
+
+def test_rasterize_existing_output(tmpdir, runner):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, '--res', 0.5], input=TEST_FEATURES)
+    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)
+
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 55
+        assert (data == 1).sum() == 125
+        assert (data == 2).sum() == 20
+
+
+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)
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 548576
+        assert (data == 1).sum() == 500000
+
+    # Test invalid like raster
+    output = str(tmpdir.join('test2.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, '--like', 'foo.tif'], input=TEST_FEATURES)
+    assert result.exit_code == 2
+
+
+def test_rasterize_property_value(tmpdir, runner):
+    # Test feature collection property values
+    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_band(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)
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 55
+        assert (data == 15).sum() == 145
+
+def test_rasterize_out_of_bounds(tmpdir, runner):
+    output = str(tmpdir.join('test.tif'))
+
+    # 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 == 2
+
+    # Test out of bounds of existing output raster (first have to create one)
+    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)
+
+    result = runner.invoke(features.rasterize, [output], input=TEST_FEATURES)
+    assert result.exit_code == 2
diff --git a/tests/test_rio_info.py b/tests/test_rio_info.py
index fce04eb..0dd5c65 100644
--- a/tests/test_rio_info.py
+++ b/tests/test_rio_info.py
@@ -71,3 +71,45 @@ def test_info_tags():
         ['tests/data/RGB.byte.tif', '--tags'])
     assert result.exit_code == 0
     assert result.output == '{"AREA_OR_POINT": "Area"}\n'
+
+
+def test_info_res():
+    runner = CliRunner()
+    result = runner.invoke(
+        info.info,
+        ['tests/data/RGB.byte.tif', '--res'])
+    assert result.exit_code == 0
+    assert result.output.startswith('300.037')
+
+
+def test_info_lnglat():
+    runner = CliRunner()
+    result = runner.invoke(
+        info.info,
+        ['tests/data/RGB.byte.tif', '--lnglat'])
+    assert result.exit_code == 0
+    assert result.output.startswith('-77.757')
+
+
+def test_mo_info():
+    runner = CliRunner()
+    result = runner.invoke(info.info, ['tests/data/RGB.byte.tif'])
+    assert result.exit_code == 0
+    assert '"res": [300.037' in result.output
+    assert '"lnglat": [-77.757' in result.output
+
+
+def test_info_stats():
+    runner = CliRunner()
+    result = runner.invoke(info.info, ['tests/data/RGB.byte.tif', 
'--tell-me-more'])
+    assert result.exit_code == 0
+    assert '"max": 255.0' in result.output
+    assert '"min": 1.0' in result.output
+    assert '"mean": 44.4344' in result.output
+
+
+def test_info_stats_only():
+    runner = CliRunner()
+    result = runner.invoke(info.info, ['tests/data/RGB.byte.tif', '--stats', 
'--bidx', '2'])
+    assert result.exit_code == 0
+    assert result.output.startswith('1.000000 255.000000 66.02')
diff --git a/tests/test_rio_sample.py b/tests/test_rio_sample.py
new file mode 100644
index 0000000..2c6a329
--- /dev/null
+++ b/tests/test_rio_sample.py
@@ -0,0 +1,81 @@
+import logging
+import sys
+
+import click
+from click.testing import CliRunner
+
+import rasterio
+from rasterio.rio import sample
+
+
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+
+def test_sample_err():
+    runner = CliRunner()
+    result = runner.invoke(
+        sample.sample,
+        ['bogus.tif'],
+        "[220650.0, 2719200.0]")
+    assert result.exit_code == 1
+
+
+def test_sample_stdin():
+    runner = CliRunner()
+    result = runner.invoke(
+        sample.sample,
+        ['tests/data/RGB.byte.tif'],
+        "[220650.0, 2719200.0]\n[220650.0, 2719200.0]",
+        catch_exceptions=False)
+    assert result.exit_code == 0
+    assert result.output.strip() == '[28, 29, 27]\n[28, 29, 27]'
+
+
+def test_sample_arg():
+    runner = CliRunner()
+    result = runner.invoke(
+        sample.sample,
+        ['tests/data/RGB.byte.tif', "[220650.0, 2719200.0]"],
+        catch_exceptions=False)
+    assert result.exit_code == 0
+    assert result.output.strip() == '[28, 29, 27]'
+
+
+def test_sample_bidx():
+    runner = CliRunner()
+    result = runner.invoke(
+        sample.sample,
+        ['tests/data/RGB.byte.tif', '--bidx', '1,2', "[220650.0, 2719200.0]"],
+        catch_exceptions=False)
+    assert result.exit_code == 0
+    assert result.output.strip() == '[28, 29]'
+
+
+def test_sample_bidx2():
+    runner = CliRunner()
+    result = runner.invoke(
+        sample.sample,
+        ['tests/data/RGB.byte.tif', '--bidx', '1..2', "[220650.0, 2719200.0]"],
+        catch_exceptions=False)
+    assert result.exit_code == 0
+    assert result.output.strip() == '[28, 29]'
+
+
+def test_sample_bidx3():
+    runner = CliRunner()
+    result = runner.invoke(
+        sample.sample,
+        ['tests/data/RGB.byte.tif', '--bidx', '..2', "[220650.0, 2719200.0]"],
+        catch_exceptions=False)
+    assert result.exit_code == 0
+    assert result.output.strip() == '[28, 29]'
+
+
+def test_sample_bidx4():
+    runner = CliRunner()
+    result = runner.invoke(
+        sample.sample,
+        ['tests/data/RGB.byte.tif', '--bidx', '3', "[220650.0, 2719200.0]"],
+        catch_exceptions=False)
+    assert result.exit_code == 0
+    assert result.output.strip() == '[27]'
diff --git a/tests/test_sampling.py b/tests/test_sampling.py
new file mode 100644
index 0000000..c4030cc
--- /dev/null
+++ b/tests/test_sampling.py
@@ -0,0 +1,17 @@
+import rasterio
+
+
+def test_sampling():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        data = next(src.sample([(220650.0, 2719200.0)]))
+        assert list(data) == [28, 29, 27]
+
+def test_sampling_beyond_bounds():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        data = next(src.sample([(-10, 2719200.0)]))
+        assert list(data) == [0, 0, 0]
+
+def test_sampling_indexes():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        data = next(src.sample([(220650.0, 2719200.0)], indexes=[2]))
+        assert list(data) == [29]
diff --git a/tests/test_tags.py b/tests/test_tags.py
index 4431055..514d068 100644
--- a/tests/test_tags.py
+++ b/tests/test_tags.py
@@ -54,3 +54,11 @@ def test_tags_update_twice():
         assert dst.tags() == {'a': '1', 'b': '2'}
         dst.update_tags(c=3)
         assert dst.tags() == {'a': '1', 'b': '2', 'c': '3'}
+
+
+def test_tags_eq():
+    with rasterio.open(
+            'test.tif', 'w', 
+            'GTiff', 3, 4, 1, dtype=rasterio.ubyte) as dst:
+        dst.update_tags(a="foo=bar")
+        assert dst.tags() == {'a': "foo=bar"}
diff --git a/tests/test_transform.py b/tests/test_transform.py
index b17feab..acd0d54 100644
--- a/tests/test_transform.py
+++ b/tests/test_transform.py
@@ -1,4 +1,3 @@
-
 import rasterio
 
 def test_window():
@@ -9,3 +8,16 @@ def test_window():
         assert src.window(left, top-src.res[1], left+src.res[0], top) == (
             (0, 1), (0, 1))
 
+
+def test_window_transform():
+    with rasterio.open('tests/data/RGB.byte.tif') as src:
+        assert src.window_transform(((0, None), (0, None))) == src.affine
+        assert src.window_transform(((None, None), (None, None))) == src.affine
+        assert src.window_transform(
+                ((1, None), (1, None))).c == src.bounds.left + src.res[0]
+        assert src.window_transform(
+                ((1, None), (1, None))).f == src.bounds.top - src.res[1]
+        assert src.window_transform(
+                ((-1, None), (-1, None))).c == src.bounds.left - src.res[0]
+        assert src.window_transform(
+                ((-1, None), (-1, None))).f == src.bounds.top + src.res[1]

-- 
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