Merge branch 'master' of github.com:mapbox/rasterio

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

rasterio/_gdal.pxd

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

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)