CHANGES.txt

 Changes
 =======
 
+1.0.12 (2018-12-10)
+-------------------
+
+- Rasterio's lower level I/O functions now take Py_ssize_t index arrays and
+  will raise exceptions for type mismatches instead of swallowing the
+  exceptions (#1076).
+- The _boundless_vrt_doc function's background layer no longer needs an
+  in-memory dataset and the performance regression noted in #1499 has been
+  reversed. We've also found a performance improvement in masked boundless
+  reads in the case that the source data is entirely valid.
+- The signature of the private _boundless_vrt_doc function in rasterio.vrt has
+  changed. Its background keyword argument now takes an int or float, a new
+  masked keyword argument has been added, and the function returns a unicode
+  str instead of ascii-encoded bytes.
+- The copy and copyfiles functions of rasterio.shutil now raise an exception
+  when the source and destination paths identify the same dataset (#1569).
+
 1.0.11 (2019-11-30)
 -------------------
 

debian/changelog

+rasterio (1.0.12-1) unstable; urgency=medium
+
+  * Team upload.
+  * New upstream release.
+
+ -- Bas Couwenberg <sebastic@debian.org>  Tue, 11 Dec 2018 09:30:56 +0100
+
 rasterio (1.0.11-1) unstable; urgency=medium
 
   * Team upload.

rasterio/__init__.py

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

rasterio/_features.pyx

@@ -310,7 +310,7 @@ def _rasterize(shapes, image, transform, all_touched, merge_alg):
             exc_wrap_int(
                 GDALRasterizeGeometries(
                     mem.handle(), 1, mem.band_ids, num_geoms, geoms, NULL,
-                    mem.gdal_transform, pixel_values, options, NULL, NULL))
+                    NULL, pixel_values, options, NULL, NULL))
 
             # Read in-memory data back into image
             image = mem.read()

rasterio/_io.pxd

@@ -50,4 +50,4 @@ ctypedef np.float64_t DTYPE_FLOAT64_t
 
 cdef bint in_dtype_range(value, dtype)
 
-cdef int io_auto(image, GDALRasterBandH band, bint write, int resampling=*)
+cdef int io_auto(image, GDALRasterBandH band, bint write, int resampling=*) except -1

rasterio/_io.pyx

@@ -104,7 +104,7 @@ cdef bint in_dtype_range(value, dtype):
     return rng.min <= value <= rng.max
 
 
-cdef int io_auto(data, GDALRasterBandH band, bint write, int resampling=0):
+cdef int io_auto(data, GDALRasterBandH band, bint write, int resampling=0) except -1:
     """Convenience function to handle IO with a GDAL band.
 
     :param data: a numpy ndarray
@@ -121,7 +121,7 @@ cdef int io_auto(data, GDALRasterBandH band, bint write, int resampling=0):
         return io_band(band, write, 0.0, 0.0, width, height, data,
                        resampling=resampling)
     elif ndims == 3:
-        indexes = np.arange(1, data.shape[0] + 1)
+        indexes = np.arange(1, data.shape[0] + 1).astype('int')
         return io_multi_band(band, write, 0.0, 0.0, width, height, data,
                              indexes, resampling=resampling)
     else:
@@ -367,23 +367,15 @@ cdef class DatasetReaderBase(DatasetBase):
         else:
 
             if fill_value is not None:
-                dtype = self.dtypes[0]
-                bg_path = UnparsedPath('/vsimem/bg{}.tif'.format(uuid.uuid4()))
-                with DatasetWriterBase(
-                        bg_path, 'w',
-                        driver='GTiff', count=self.count, height=3, width=3,
-                        dtype=dtype, crs=None, transform=Affine.identity() * Affine.translation(1, 1)) as bg_dataset:
-                    bg_dataset.write(
-                        np.full((self.count, 3, 3), fill_value, dtype=dtype))
-                bg_dataset = DatasetReaderBase(bg_path)
+                nodataval = fill_value
             else:
-                bg_dataset = None
+                nodataval = ndv
 
             vrt_doc = _boundless_vrt_doc(
-                self, nodata=ndv, background=bg_dataset,
+                self, nodata=nodataval, background=nodataval,
                 width=max(self.width, window.width) + 1,
                 height=max(self.height, window.height) + 1,
-                transform=self.window_transform(window)).decode('ascii')
+                transform=self.window_transform(window))
 
             if not gdal_version().startswith('1'):
                 vrt_kwds = {'driver': 'VRT'}
@@ -396,58 +388,43 @@ cdef class DatasetReaderBase(DatasetBase):
                     indexes, out, Window(0, 0, window.width, window.height),
                     None, resampling=resampling)
 
-                if masked and all_valid:
-                    blank_path = UnparsedPath('/vsimem/blank-{}.tif'.format(uuid.uuid4()))
-                    transform = Affine.translation(self.transform.xoff, self.transform.yoff) * (Affine.scale(self.width / 3, self.height / 3) * (Affine.translation(-self.transform.xoff, -self.transform.yoff) * self.transform))
-                    with DatasetWriterBase(
-                            blank_path, 'w',
-                            driver='GTiff', count=self.count, height=3, width=3,
-                            dtype='uint8', crs=self.crs, transform=transform) as blank_dataset:
-                        blank_dataset.write(
-                            np.ones((self.count, 3, 3), dtype='uint8'))
+                if masked:
+
+                    # Below we use another VRT to compute the valid data mask
+                    # in this special case where all source pixels are valid.
+                    if all_valid:
 
-                    with DatasetReaderBase(blank_path) as blank_dataset:
                         mask_vrt_doc = _boundless_vrt_doc(
-                            blank_dataset, nodata=0,
+                            self, nodata=0,
                             width=max(self.width, window.width) + 1,
                             height=max(self.height, window.height) + 1,
-                            transform=self.window_transform(window)).decode('ascii')
+                            transform=self.window_transform(window),
+                            masked=True)
 
                         with DatasetReaderBase(UnparsedPath(mask_vrt_doc), **vrt_kwds) as mask_vrt:
                             mask = np.zeros(out.shape, 'uint8')
                             mask = ~mask_vrt._read(
                                 indexes, mask, Window(0, 0, window.width, window.height), None).astype('bool')
 
-                            kwds = {'mask': mask}
-                            # Set a fill value only if the read bands share a
-                            # single nodata value.
-                            if fill_value is not None:
-                                kwds['fill_value'] = fill_value
-                            elif len(set(nodatavals)) == 1:
-                                if nodatavals[0] is not None:
-                                    kwds['fill_value'] = nodatavals[0]
 
-                            out = np.ma.array(out, **kwds)
-
-                elif masked:
+                    else:
                         mask = np.zeros(out.shape, 'uint8')
                         mask = ~vrt._read(
                             indexes, mask, Window(0, 0, window.width, window.height), None, masks=True).astype('bool')
 
                     kwds = {'mask': mask}
+
                     # Set a fill value only if the read bands share a
                     # single nodata value.
                     if fill_value is not None:
                         kwds['fill_value'] = fill_value
+
                     elif len(set(nodatavals)) == 1:
                         if nodatavals[0] is not None:
                             kwds['fill_value'] = nodatavals[0]
 
                     out = np.ma.array(out, **kwds)
 
-            if bg_dataset is not None:
-                bg_dataset.close()
-
         if return2d:
             out.shape = out.shape[1:]
 
@@ -601,7 +578,7 @@ cdef class DatasetReaderBase(DatasetBase):
                         blank_dataset, nodata=0,
                         width=max(self.width, window.width) + 1,
                         height=max(self.height, window.height) + 1,
-                        transform=self.window_transform(window)).decode('ascii')
+                        transform=self.window_transform(window))
 
                     with DatasetReaderBase(UnparsedPath(mask_vrt_doc), **vrt_kwds) as mask_vrt:
                         out = np.zeros(out.shape, 'uint8')
@@ -612,7 +589,7 @@ cdef class DatasetReaderBase(DatasetBase):
                 vrt_doc = _boundless_vrt_doc(
                     self, width=max(self.width, window.width) + 1,
                     height=max(self.height, window.height) + 1,
-                    transform=self.window_transform(window)).decode('ascii')
+                    transform=self.window_transform(window))
 
                 with DatasetReaderBase(UnparsedPath(vrt_doc), **vrt_kwds) as vrt:
 
@@ -693,7 +670,7 @@ cdef class DatasetReaderBase(DatasetBase):
 
         # Call io_multi* functions with C type args so that they
         # can release the GIL.
-        indexes_arr = np.array(indexes, dtype=int)
+        indexes_arr = np.array(indexes).astype('int')
         indexes_count = <int>indexes_arr.shape[0]
 
         if masks:

rasterio/_shim.pxd

 include "gdal.pxi"
 
-cdef GDALDatasetH open_dataset(object filename, int mode,
-                               object allowed_drivers, object open_options,
-                               object siblings) except NULL
+cdef GDALDatasetH open_dataset(object filename, int mode, object allowed_drivers, object open_options, object siblings) except NULL
 cdef int delete_nodata_value(GDALRasterBandH hBand) except 3
-cdef int io_band(GDALRasterBandH band, int mode, float xoff, float yoff,
-                 float width, float height, object data, int resampling=*)
-cdef int io_multi_band(GDALDatasetH hds, int mode, float xoff, float yoff,
-                       float width, float height, object data, long[:] indexes,
-                       int resampling=*)
-cdef int io_multi_mask(GDALDatasetH hds, int mode, float xoff, float yoff,
-                       float width, float height, object data, long[:] indexes,
-                       int resampling=*)
+cdef int io_band(GDALRasterBandH band, int mode, float xoff, float yoff, float width, float height, object data, int resampling=*) except -1
+cdef int io_multi_band(GDALDatasetH hds, int mode, float xoff, float yoff, float width, float height, object data, Py_ssize_t[:] indexes, int resampling=*) except -1
+cdef int io_multi_mask(GDALDatasetH hds, int mode, float xoff, float yoff, float width, float height, object data, Py_ssize_t[:] indexes, int resampling=*) except -1

rasterio/_shim1.pyx

@@ -55,8 +55,9 @@ cdef int delete_nodata_value(GDALRasterBandH hBand) except 3:
         "GDAL versions < 2.1 do not support nodata deletion")
 
 
-cdef int io_band(GDALRasterBandH band, int mode, float x0, float y0,
-                 float width, float height, object data, int resampling=0):
+cdef int io_band(
+        GDALRasterBandH band, int mode, float x0, float y0,
+        float width, float height, object data, int resampling=0) except -1:
     """Read or write a region of data for the band.
 
     Implicit are
@@ -90,9 +91,9 @@ cdef int io_band(GDALRasterBandH band, int mode, float x0, float y0,
     return retval
 
 
-cdef int io_multi_band(GDALDatasetH hds, int mode, float x0, float y0,
-                       float width, float height, object data,
-                       long[:] indexes, int resampling=0):
+cdef int io_multi_band(
+        GDALDatasetH hds, int mode, float x0, float y0, float width,
+        float height, object data, Py_ssize_t[:] indexes, int resampling=0) except -1:
     """Read or write a region of data for multiple bands.
 
     Implicit are
@@ -123,7 +124,7 @@ cdef int io_multi_band(GDALDatasetH hds, int mode, float x0, float y0,
     with nogil:
         bandmap = <int *>CPLMalloc(count*sizeof(int))
         for i in range(count):
-            bandmap[i] = indexes[i]
+            bandmap[i] = <int>indexes[i]
         retval = GDALDatasetRasterIO(
             hds, mode, xoff, yoff, xsize, ysize, buf,
             bufxsize, bufysize, buftype, count, bandmap,
@@ -133,9 +134,9 @@ cdef int io_multi_band(GDALDatasetH hds, int mode, float x0, float y0,
     return retval
 
 
-cdef int io_multi_mask(GDALDatasetH hds, int mode, float x0, float y0,
-                       float width, float height, object data,
-                       long[:] indexes, int resampling=0):
+cdef int io_multi_mask(
+        GDALDatasetH hds, int mode, float x0, float y0, float width,
+        float height, object data, Py_ssize_t[:] indexes, int resampling=0) except -1:
     """Read or write a region of data for multiple band masks.
 
     Implicit are
@@ -165,7 +166,7 @@ cdef int io_multi_mask(GDALDatasetH hds, int mode, float x0, float y0,
     cdef int ysize = <int>height
 
     for i in range(count):
-        j = indexes[i]
+        j = <int>indexes[i]
         band = GDALGetRasterBand(hds, j)
         if band == NULL:
             raise ValueError("Null band")

rasterio/_shim20.pyx

@@ -12,7 +12,7 @@ include "shim_rasterioex.pxi"
 # Declarations and implementations specific for GDAL = 2.0
 cdef extern from "gdal.h" nogil:
 
-    GDALDatasetH GDALOpenEx(const char *filename, int flags, const char **allowed_drivers, const char **options, const char **siblings) # except -1
+    GDALDatasetH GDALOpenEx(const char *filename, int flags, const char **allowed_drivers, const char **options, const char **siblings)
 
 from rasterio._err cimport exc_wrap_pointer
 

rasterio/rio/warp.py

 """$ rio warp"""
 
-import logging
-from math import ceil, log
-import warnings
+from math import ceil
 
 import click
 from cligj import format_opt
 
 import rasterio
 from rasterio.crs import CRS
-from rasterio.env import setenv, GDALVersion
+from rasterio.env import setenv
 from rasterio.errors import CRSError
 from rasterio.rio import options
 from rasterio.rio.helpers import resolve_inout
@@ -150,8 +148,8 @@ def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds,
 
         with rasterio.open(files[0]) as src:
             l, b, r, t = src.bounds
-            out_kwargs = src.profile.copy()
-            out_kwargs['driver'] = driver
+            out_kwargs = src.profile
+            out_kwargs.update(driver=driver)
 
             # Sort out the bounds options.
             if src_bounds and dst_bounds:
@@ -271,9 +269,7 @@ def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds,
             # value to src_nodata (will be overridden by dst_nodata if it is not None
             if src_nodata is not None:
                 # Update the dst nodata value
-                out_kwargs.update({
-                    'nodata': src_nodata
-                })
+                out_kwargs.update(nodata=src_nodata)
 
             # Validate a manually set destination NODATA value
             # against the input datatype.
@@ -283,7 +279,7 @@ def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds,
                         "--src-nodata must be provided because dst-nodata is not None")
                 else:
                     # Update the dst nodata value
-                    out_kwargs.update({'nodata': dst_nodata})
+                    out_kwargs.update(nodata=dst_nodata)
 
             # When the bounds option is misused, extreme values of
             # destination width and height may result.
@@ -294,12 +290,12 @@ def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds,
                     "Invalid output dimensions: {0}.".format(
                         (dst_width, dst_height)))
 
-            out_kwargs.update({
-                'crs': dst_crs,
-                'transform': dst_transform,
-                'width': dst_width,
-                'height': dst_height
-            })
+            out_kwargs.update(
+                crs=dst_crs,
+                transform=dst_transform,
+                width=dst_width,
+                height=dst_height
+            )
 
             # Adjust block size if necessary.
             if ('blockxsize' in out_kwargs and

rasterio/shim_rasterioex.pxi

@@ -31,7 +31,7 @@ cdef extern from "gdal.h" nogil:
 
 
 cdef int io_band(GDALRasterBandH band, int mode, float x0, float y0,
-                 float width, float height, object data, int resampling=0):
+                 float width, float height, object data, int resampling=0) except -1:
     """Read or write a region of data for the band.
 
     Implicit are
@@ -78,7 +78,7 @@ cdef int io_band(GDALRasterBandH band, int mode, float x0, float y0,
 
 cdef int io_multi_band(GDALDatasetH hds, int mode, float x0, float y0,
                        float width, float height, object data,
-                       long[:] indexes, int resampling=0):
+                       Py_ssize_t[:] indexes, int resampling=0) except -1:
     """Read or write a region of data for multiple bands.
 
     Implicit are
@@ -120,7 +120,7 @@ cdef int io_multi_band(GDALDatasetH hds, int mode, float x0, float y0,
     with nogil:
         bandmap = <int *>CPLMalloc(count*sizeof(int))
         for i in range(count):
-            bandmap[i] = indexes[i]
+            bandmap[i] = <int>indexes[i]
         retval = GDALDatasetRasterIOEx(
             hds, <GDALRWFlag>mode, xoff, yoff, xsize, ysize, buf,
             bufxsize, bufysize, buftype, count, bandmap,
@@ -132,7 +132,7 @@ cdef int io_multi_band(GDALDatasetH hds, int mode, float x0, float y0,
 
 cdef int io_multi_mask(GDALDatasetH hds, int mode, float x0, float y0,
                        float width, float height, object data,
-                       long[:] indexes, int resampling=0):
+                       Py_ssize_t[:] indexes, int resampling=0) except -1:
     """Read or write a region of data for multiple band masks.
 
     Implicit are
@@ -173,7 +173,7 @@ cdef int io_multi_mask(GDALDatasetH hds, int mode, float x0, float y0,
     extras.pProgressData = NULL
 
     for i in range(count):
-        j = indexes[i]
+        j = <int>indexes[i]
         band = GDALGetRasterBand(hds, j)
         if band == NULL:
             raise ValueError("Null band")

rasterio/shutil.pyx

-"""Raster management."""
-
+"""Raster file management."""
 
 include "gdal.pxi"
 
 import logging
 
+try:
+    from pathlib import Path
+except ImportError:  # pragma: no cover
+    class Path:
+        pass
+
 from rasterio._io cimport DatasetReaderBase
 from rasterio._err cimport exc_wrap_int, exc_wrap_pointer
 from rasterio.env import ensure_env_with_credentials
@@ -55,9 +60,9 @@ def copy(src, dst, driver='GTiff', strict=True, **creation_options):
 
     Parameters
     ----------
-    src : str or dataset object opened in 'r' mode
+    src : str or pathlib.Path or dataset object opened in 'r' mode
         Source dataset
-    dst : str
+    dst : str or pathlib.Path
         Output dataset path
     driver : str, optional
         Output driver name
@@ -66,6 +71,11 @@ def copy(src, dst, driver='GTiff', strict=True, **creation_options):
         driver may adapt as necessary
     creation_options : **kwargs, optional
         Creation options for output dataset
+
+    Returns
+    -------
+    None
+
     """
 
     cdef bint c_strictness
@@ -87,29 +97,47 @@ def copy(src, dst, driver='GTiff', strict=True, **creation_options):
     c_strictness = strict
     driverb = driver.encode('utf-8')
     drv = GDALGetDriverByName(driverb)
+
     if drv == NULL:
         raise DriverRegistrationError("Unrecognized driver: {}".format(driver))
 
-    # Input is a path or GDAL connection string
+    # Convert src and dst Paths to strings.
+    if isinstance(src, Path):
+        src = str(src)
+    if isinstance(dst, Path):
+        dst = str(dst)
+
+    # Open a new GDAL dataset if src is a string.
     if isinstance(src, str):
+
+        if vsi_path(parse_path(src)) == vsi_path(parse_path(dst)):
+            raise RasterioIOError("{} and {} identify the same dataset.".format(src, dst))
+
         src = src.encode('utf-8')
         c_src_path = src
         with nogil:
             src_dataset = GDALOpenShared(c_src_path, <GDALAccess>0)
         src_dataset = exc_wrap_pointer(src_dataset)
         close_src = True
-    # Input is something like 'rasterio.open()'
+
+    # Try to use the existing GDAL dataset handle otherwise.
     else:
+
+        if src.name == vsi_path(parse_path(dst)):
+            raise RasterioIOError("{} and {} identify the same dataset.".format(src.name, dst))
+
         src_dataset = (<DatasetReaderBase?>src).handle()
         close_src = False
 
     dst = dst.encode('utf-8')
     c_dst_path = dst
+
     try:
         with nogil:
             dst_dataset = GDALCreateCopy(
                 drv, c_dst_path, src_dataset, c_strictness, options, NULL, NULL)
         dst_dataset = exc_wrap_pointer(dst_dataset)
+
     finally:
         CSLDestroy(options)
         with nogil:
@@ -127,19 +155,35 @@ def copyfiles(src, dst):
 
     Parameters
     ----------
-    src : str
+    src : str or pathlib.Path
         Source dataset
-    dst : str
+    dst : str or pathlib.Path
         Target dataset
+
+    Returns
+    -------
+    None
+
     """
 
     cdef GDALDatasetH h_dataset = NULL
     cdef GDALDriverH h_driver = NULL
 
+    # Convert src and dst Paths to strings.
+    if isinstance(src, Path):
+        src = str(src)
+    if isinstance(dst, Path):
+        dst = str(dst)
+
+    src_path = parse_path(src)
+    dst_path = parse_path(dst)
+    if vsi_path(src_path) == vsi_path(dst_path):
+        raise RasterioIOError("{} and {} identify the same dataset.".format(src, dst))
+
     # VFS paths probabaly don't work, but its hard to be completely certain
     # so just attempt to use them.
-    gdal_src_path = vsi_path(parse_path(src))
-    gdal_dst_path = vsi_path(parse_path(dst))
+    gdal_src_path = vsi_path(src_path)
+    gdal_dst_path = vsi_path(dst_path)
     b_gdal_src_path = gdal_src_path.encode('utf-8')
     b_gdal_dst_path = gdal_dst_path.encode('utf-8')
     cdef char* c_gdal_src_path = b_gdal_src_path

rasterio/vrt.py

@@ -71,21 +71,22 @@ class WarpedVRT(WarpedVRTReaderBase, WindowMethodsMixin,
 
 def _boundless_vrt_doc(
         src_dataset, nodata=None, background=None, hidenodata=False,
-        width=None, height=None, transform=None):
+        width=None, height=None, transform=None, masked=False):
     """Make a VRT XML document.
 
     Parameters
     ----------
     src_dataset : Dataset
         The dataset to wrap.
-    background : Dataset, optional
-        A dataset that provides the optional VRT background. NB: this dataset
-        must have the same number of bands as the src_dataset.
+    background : int or float, optional
+        The background fill value for the boundless VRT.
+    masked : book
+        If True, the src_dataset is replaced by its valid data mask.
 
     Returns
     -------
-    bytes
-        An ascii-encoded string (an ElementTree detail)
+    str
+        An XML text string.
     """
 
     nodata = nodata or src_dataset.nodata
@@ -118,56 +119,69 @@ def _boundless_vrt_doc(
         colorinterp.text = ci.name.capitalize()
 
         if background is not None:
-            simplesource = ET.SubElement(vrtrasterband, 'SimpleSource')
-            sourcefilename = ET.SubElement(simplesource, 'SourceFilename')
-            sourcefilename.attrib['relativeToVRT'] = "0"
-            sourcefilename.text = vsi_path(parse_path(background.name))
-            sourceband = ET.SubElement(simplesource, 'SourceBand')
+            complexsource = ET.SubElement(vrtrasterband, 'ComplexSource')
+            sourcefilename = ET.SubElement(complexsource, 'SourceFilename')
+            sourcefilename.attrib['relativeToVRT'] = '1'
+            sourcefilename.text = 'dummy.tif'  # vsi_path(parse_path(background.name))
+            sourceband = ET.SubElement(complexsource, 'SourceBand')
             sourceband.text = str(bidx)
-            sourceproperties = ET.SubElement(simplesource, 'SourceProperties')
+            sourceproperties = ET.SubElement(complexsource, 'SourceProperties')
             sourceproperties.attrib['RasterXSize'] = str(width)
             sourceproperties.attrib['RasterYSize'] = str(height)
             sourceproperties.attrib['dataType'] = _gdal_typename(dtype)
             sourceproperties.attrib['BlockYSize'] = str(block_shape[0])
             sourceproperties.attrib['BlockXSize'] = str(block_shape[1])
-            srcrect = ET.SubElement(simplesource, 'SrcRect')
+            srcrect = ET.SubElement(complexsource, 'SrcRect')
             srcrect.attrib['xOff'] = '0'
             srcrect.attrib['yOff'] = '0'
-            srcrect.attrib['xSize'] = str(background.width)
-            srcrect.attrib['ySize'] = str(background.height)
-            dstrect = ET.SubElement(simplesource, 'DstRect')
+            srcrect.attrib['xSize'] = '1'  # str(background.width)
+            srcrect.attrib['ySize'] = '1'  # str(background.height)
+            dstrect = ET.SubElement(complexsource, 'DstRect')
             dstrect.attrib['xOff'] = '0'
             dstrect.attrib['yOff'] = '0'
-            dstrect.attrib['xSize'] = str(width)
-            dstrect.attrib['ySize'] = str(height)
-
-        simplesource = ET.SubElement(vrtrasterband, 'SimpleSource')
-        sourcefilename = ET.SubElement(simplesource, 'SourceFilename')
+            dstrect.attrib['xSize'] = '1'  # str(width)
+            dstrect.attrib['ySize'] = '1'  # str(height)
+            scaleratio = ET.SubElement(complexsource, 'ScaleRatio')
+            scaleratio.text = '0'
+            scaleoffset = ET.SubElement(complexsource, 'ScaleOffset')
+            scaleoffset.text = str(background)
+
+        complexsource = ET.SubElement(vrtrasterband, 'ComplexSource')
+        sourcefilename = ET.SubElement(complexsource, 'SourceFilename')
         sourcefilename.attrib['relativeToVRT'] = "0"
         sourcefilename.text = vsi_path(parse_path(src_dataset.name))
-        sourceband = ET.SubElement(simplesource, 'SourceBand')
+        sourceband = ET.SubElement(complexsource, 'SourceBand')
         sourceband.text = str(bidx)
-        sourceproperties = ET.SubElement(simplesource, 'SourceProperties')
+        sourceproperties = ET.SubElement(complexsource, 'SourceProperties')
         sourceproperties.attrib['RasterXSize'] = str(width)
         sourceproperties.attrib['RasterYSize'] = str(height)
         sourceproperties.attrib['dataType'] = _gdal_typename(dtype)
         sourceproperties.attrib['BlockYSize'] = str(block_shape[0])
         sourceproperties.attrib['BlockXSize'] = str(block_shape[1])
-        srcrect = ET.SubElement(simplesource, 'SrcRect')
+        srcrect = ET.SubElement(complexsource, 'SrcRect')
         srcrect.attrib['xOff'] = '0'
         srcrect.attrib['yOff'] = '0'
         srcrect.attrib['xSize'] = str(src_dataset.width)
         srcrect.attrib['ySize'] = str(src_dataset.height)
-        dstrect = ET.SubElement(simplesource, 'DstRect')
+        dstrect = ET.SubElement(complexsource, 'DstRect')
         dstrect.attrib['xOff'] = str((src_dataset.transform.xoff - transform.xoff) / transform.a)
         dstrect.attrib['yOff'] = str((src_dataset.transform.yoff - transform.yoff) / transform.e)
         dstrect.attrib['xSize'] = str(src_dataset.width * src_dataset.transform.a / transform.a)
         dstrect.attrib['ySize'] = str(src_dataset.height * src_dataset.transform.e / transform.e)
 
         if src_dataset.nodata is not None:
-            nodata_elem = ET.SubElement(simplesource, 'NODATA')
+            nodata_elem = ET.SubElement(complexsource, 'NODATA')
             nodata_elem.text = str(src_dataset.nodata)
 
+        # Effectively replaces all values of the source dataset with
+        # 255.  Due to GDAL optimizations, the source dataset will not
+        # be read, so we get a performance improvement.
+        if masked:
+            scaleratio = ET.SubElement(complexsource, 'ScaleRatio')
+            scaleratio.text = '0'
+            scaleoffset = ET.SubElement(complexsource, 'ScaleOffset')
+            scaleoffset.text = '255'
+
     if all(MaskFlags.per_dataset in flags for flags in src_dataset.mask_flag_enums):
         maskband = ET.SubElement(vrtdataset, 'MaskBand')
         vrtrasterband = ET.SubElement(maskband, 'VRTRasterBand')
@@ -197,4 +211,4 @@ def _boundless_vrt_doc(
         dstrect.attrib['xSize'] = str(src_dataset.width)
         dstrect.attrib['ySize'] = str(src_dataset.height)
 
-    return ET.tostring(vrtdataset)
+    return ET.tostring(vrtdataset).decode('ascii')

rasterio/warp.py

@@ -213,7 +213,7 @@ def reproject(source, destination, src_transform=None, gcps=None,
         Example: CRS({'init': 'EPSG:4326'})
     src_nodata: int or float, optional
         The source nodata value. Pixels with this value will not be
-        used for interpolation. If not set, it will be default to the
+        used for interpolation. If not set, it will default to the
         nodata value of the source image if a masked ndarray or
         rasterio band, if available.
     dst_transform: affine.Affine(), optional

tests/test_overviews.py

@@ -2,6 +2,8 @@
 
 import pytest
 
+from .conftest import requires_gdal2
+
 import rasterio
 from rasterio.enums import Resampling
 from rasterio.env import GDALVersion
@@ -90,6 +92,7 @@ def test_issue1333(data):
                 overview_factors, resampling=Resampling.average)
 
 
+@requires_gdal2
 def test_build_overviews_new_file(tmpdir, path_rgb_byte_tif):
     """Confirm fix of #1497"""
     dst_file = str(tmpdir.join('test.tif'))

tests/test_shutil.py

@@ -52,6 +52,14 @@ def test_exists(path_rgb_byte_tif):
     assert not rasterio.shutil.exists('trash')
 
 
+def test_copy_fail_same_dataset(data):
+    """A dataset can't be copied to itself."""
+    path = str(data.join('RGB.byte.tif'))
+    with pytest.raises(RasterioIOError):
+        with rasterio.open(path) as src:
+            rasterio.shutil.copy(src, path, **src.profile)
+
+
 @pytest.mark.parametrize("pass_handle", (True, False))
 def test_copy(tmpdir, path_rgb_byte_tif, pass_handle):
 
@@ -143,6 +151,20 @@ def test_copyfiles(data, tmpdir):
         assert all(map(os.path.exists, actual.files))
 
 
+def test_copyfiles_same_dataset(data):
+    """A dataset can't be copied to itself."""
+    path = str(data.join('RGB.byte.tif'))
+    with pytest.raises(RasterioIOError):
+        rasterio.shutil.copyfiles(path, path)
+
+
+def test_copyfiles_same_dataset_another_name(data):
+    """A dataset can't be copied to itself by another name."""
+    path = str(data.join('RGB.byte.tif'))
+    with pytest.raises(RasterioIOError):
+        rasterio.shutil.copyfiles(path, 'file://' + path)
+
+
 def test_copyfiles_fail():
     with pytest.raises(RasterioIOError):
         rasterio.shutil.copyfiles('trash', 'whatever')

tests/test_vrt.py

@@ -7,8 +7,8 @@ import rasterio.vrt
 def test_boundless_vrt(path_rgb_byte_tif):
     with rasterio.open(path_rgb_byte_tif) as rgb:
         doc = rasterio.vrt._boundless_vrt_doc(rgb)
-        assert doc.startswith(b'<VRTDataset')
-        with rasterio.open(doc.decode('ascii')) as vrt:
+        assert doc.startswith('<VRTDataset')
+        with rasterio.open(doc) as vrt:
             assert rgb.count == vrt.count
             assert rgb.dtypes == vrt.dtypes
             assert rgb.mask_flag_enums == vrt.mask_flag_enums
@@ -17,8 +17,8 @@ def test_boundless_vrt(path_rgb_byte_tif):
 def test_boundless_msk_vrt(path_rgb_msk_byte_tif):
     with rasterio.open(path_rgb_msk_byte_tif) as msk:
         doc = rasterio.vrt._boundless_vrt_doc(msk)
-        assert doc.startswith(b'<VRTDataset')
-        with rasterio.open(doc.decode('ascii')) as vrt:
+        assert doc.startswith('<VRTDataset')
+        with rasterio.open(doc) as vrt:
             assert msk.count == vrt.count
             assert msk.dtypes == vrt.dtypes
             assert msk.mask_flag_enums == vrt.mask_flag_enums

tests/test_warp.py

@@ -1472,3 +1472,23 @@ def test_issue_1446_b():
     }
     # Before the fix, this geometry was thrown eastward of 0.0. It should be between -350 and -250.
     assert all([-350 < x < -150 for x, y in transformed_geoms[183519]["coordinates"]])
+
+
+def test_issue_1076():
+    """Confirm fix of #1076"""
+    arr = (np.random.random((20, 30)) * 100).astype('int32')
+    fill_value = 42
+    newarr = np.full((200, 300), fill_value=fill_value, dtype='int32')
+
+    src_crs = {'init': 'EPSG:32632'}
+    src_transform = Affine(600.0, 0.0, 399960.0, 0.0, -600.0, 6100020.0)
+    dst_transform = Affine(60.0, 0.0, 399960.0, 0.0, -60.0, 6100020.0)
+
+    reproject(arr, newarr,
+        src_transform=src_transform,
+        dst_transform=dst_transform,
+        src_crs=src_crs,
+        dst_crs=src_crs,
+        resample=Resampling.nearest)
+
+    assert not (newarr == fill_value).all()