CHANGES.txt

 Changes
 =======
 
+1.0.27 (2019-09-05)
+-------------------
+
+- Resolve #1744 by adding a `dtype` keyword argument to the WarpedVRT
+  constructor. It allows a user to specify the working data type for the warp
+  operation and output.
+- All cases of deprecated affine right multiplication have been changed to be
+  forward compatible with affine 3.0. The rasterio tests now pass without
+  warnings.
+- The coordinate transformer used in _base._transform() is now properly
+  deleted, fixing the memory leak reported in #1713.
+- An unavoidable warning about 4-channel colormap entries in
+  DatasetWriterBase.write_colormap() has been removed.
+- All deprecated imports of abstract base classes for collections have been
+  corrected, eliminating the warnings reported in #1742 and #1764.
+- DatasetWriterBase no longer requires that GeoTIFF block sizes be smaller than
+  the raster size (#1760). Block sizes are however checked to ensure that they
+  are multiples of 16.
+- DatasetBase.is_tiled has been made more reliable, fixing #1376.
+- Tests have been added to demonstrate that image corruption when writing
+  block-wise to an image with extra large block sizes (#520) is no longer an
+  issue.
+
 1.0.26 (2019-08-26)
 -------------------
 

debian/changelog

+rasterio (1.0.27-1) unstable; urgency=medium
+
+  * Team upload.
+  * New upstream release.
+
+ -- Bas Couwenberg <sebastic@debian.org>  Fri, 06 Sep 2019 07:03:26 +0200
+
 rasterio (1.0.26-1) unstable; urgency=medium
 
   * Team upload.

rasterio/__init__.py

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

rasterio/_base.pyx

@@ -879,7 +879,15 @@ cdef class DatasetBase(object):
     def is_tiled(self):
         if len(self.block_shapes) == 0:
             return False
-        return self.block_shapes[0][1] < self.width and self.block_shapes[0][1] <= 1024
+        else:
+            blockysize, blockxsize = self.block_shapes[0]
+            if blockxsize % 16 or blockysize % 16:
+                return False
+            # Perfectly square is a special case/
+            if blockxsize == blockysize == self.height == self.width:
+                return True
+            else:
+                return blockxsize < self.width or blockxsize > self.width
 
     property profile:
         """Basic metadata and creation options of this dataset.
@@ -1269,17 +1277,6 @@ def _transform(src_crs, dst_crs, xs, ys, zs):
         transform = exc_wrap_pointer(transform)
         exc_wrap_int(OCTTransform(transform, n, x, y, z))
 
-    except CPLE_BaseError as exc:
-        log.debug("{}".format(exc))
-
-    except:
-        CPLFree(x)
-        CPLFree(y)
-        CPLFree(z)
-        _safe_osr_release(src)
-        _safe_osr_release(dst)
-
-    try:
         res_xs = [0]*n
         res_ys = [0]*n
         for i in range(n):
@@ -1297,6 +1294,7 @@ def _transform(src_crs, dst_crs, xs, ys, zs):
         CPLFree(x)
         CPLFree(y)
         CPLFree(z)
+        OCTDestroyCoordinateTransformation(transform)
         _safe_osr_release(src)
         _safe_osr_release(dst)
 

rasterio/_features.pyx

@@ -367,7 +367,7 @@ def _bounds(geometry, north_up=True, transform=None):
     else:
         if transform is not None:
             xyz = list(_explode(geometry['coordinates']))
-            xyz_px = [point * transform for point in xyz]
+            xyz_px = [transform * point for point in xyz]
             xyz = tuple(zip(*xyz_px))
             return min(xyz[0]), max(xyz[1]), max(xyz[0]), min(xyz[1])
         else:

rasterio/_io.pyx

@@ -25,7 +25,7 @@ from rasterio.enums import ColorInterp, MaskFlags, Resampling
 from rasterio.errors import (
     CRSError, DriverRegistrationError, RasterioIOError,
     NotGeoreferencedWarning, NodataShadowWarning, WindowError,
-    UnsupportedOperation, OverviewCreationError
+    UnsupportedOperation, OverviewCreationError, RasterBlockError
 )
 from rasterio.sample import sample_gen
 from rasterio.transform import Affine
@@ -1077,7 +1077,14 @@ cdef class DatasetWriterBase(DatasetReaderBase):
         # Process dataset opening options.
         # "tiled" affects the meaning of blocksize, so we need it
         # before iterating.
-        tiled = bool(kwargs.get('tiled', False))
+        tiled = kwargs.pop("tiled", False) or kwargs.pop("TILED", False)
+
+        if tiled:
+            blockxsize = kwargs.get("blockxsize", None)
+            blockysize = kwargs.get("blockysize", None)
+            if (blockxsize and blockxsize % 16) or (blockysize and blockysize % 16):
+                raise RasterBlockError("The height and width of dataset blocks must be multiples of 16")
+            kwargs["tiled"] = "TRUE"
 
         for k, v in kwargs.items():
             # Skip items that are definitely *not* valid driver
@@ -1087,11 +1094,8 @@ cdef class DatasetWriterBase(DatasetReaderBase):
 
             k, v = k.upper(), str(v).upper()
 
-            # Guard against block size that exceed image size.
-            if k == 'BLOCKXSIZE' and tiled and int(v) > width:
-                raise ValueError("blockxsize exceeds raster width.")
-            if k == 'BLOCKYSIZE' and tiled and int(v) > height:
-                raise ValueError("blockysize exceeds raster height.")
+            if k in ['BLOCKXSIZE', 'BLOCKYSIZE'] and not tiled:
+                continue
 
             key_b = k.encode('utf-8')
             val_b = v.encode('utf-8')
@@ -1505,12 +1509,8 @@ cdef class DatasetWriterBase(DatasetReaderBase):
         vals = range(256)
 
         for i, rgba in colormap.items():
-            if len(rgba) == 4 and self.driver in ('GTiff'):
-                warnings.warn(
-                    "This format doesn't support alpha in colormap entries. "
-                    "The value will be ignored.")
 
-            elif len(rgba) == 3:
+            if len(rgba) == 3:
                 rgba = tuple(rgba) + (255,)
 
             if i not in vals:

rasterio/_warp.pyx

@@ -112,7 +112,7 @@ def _transform_geom(
 
 cdef GDALWarpOptions * create_warp_options(
         GDALResampleAlg resampling, object src_nodata, object dst_nodata, int src_count,
-        object dst_alpha, object src_alpha, int warp_mem_limit, const char **options) except NULL:
+        object dst_alpha, object src_alpha, int warp_mem_limit, GDALDataType working_data_type, const char **options) except NULL:
     """Return a pointer to a GDALWarpOptions composed from input params
 
     This is used in _reproject() and the WarpedVRT constructor. It sets
@@ -145,6 +145,7 @@ cdef GDALWarpOptions * create_warp_options(
 
     warp_extras = CSLMerge(warp_extras, <char **>options)
 
+    psWOptions.eWorkingDataType = <GDALDataType>working_data_type
     psWOptions.eResampleAlg = <GDALResampleAlg>resampling
 
     if warp_mem_limit > 0:
@@ -212,6 +213,7 @@ def _reproject(
         init_dest_nodata=True,
         num_threads=1,
         warp_mem_limit=0,
+        working_data_type=0,
         **kwargs):
     """
     Reproject a source raster to a destination raster.
@@ -461,7 +463,7 @@ def _reproject(
 
     psWOptions = create_warp_options(
         <GDALResampleAlg>resampling, src_nodata,
-        dst_nodata, src_count, dst_alpha, src_alpha, warp_mem_limit,
+        dst_nodata, src_count, dst_alpha, src_alpha, warp_mem_limit, <GDALDataType>working_data_type,
         <const char **>warp_extras)
 
     psWOptions.pfnTransformer = pfnTransformer
@@ -607,7 +609,7 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
                  dst_width=None, width=None, dst_height=None, height=None,
                  src_transform=None, dst_transform=None, transform=None,
                  init_dest_nodata=True, src_alpha=0, add_alpha=False,
-                 warp_mem_limit=0, **warp_extras):
+                 warp_mem_limit=0, dtype=None, **warp_extras):
         """Make a virtual warped dataset
 
         Parameters
@@ -655,6 +657,8 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
         warp_mem_limit : int, optional
             The warp operation's memory limit in MB. The default (0)
             means 64 MB with GDAL 2.2.
+        dtype : str, optional
+            The working data type for warp operation and output.
         warp_extras : dict
             GDAL extra warp options. See
             http://www.gdal.org/structGDALWarpOptions.html.
@@ -735,6 +739,7 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
         self.name = "WarpedVRT({})".format(src_dataset.name)
         self.resampling = resampling
         self.tolerance = tolerance
+        self.working_dtype = dtype
 
         self.src_nodata = self.src_dataset.nodata if src_nodata is DEFAULT_NODATA_FLAG else src_nodata
         self.dst_nodata = self.src_nodata if nodata is DEFAULT_NODATA_FLAG else nodata
@@ -839,7 +844,8 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
         psWOptions = create_warp_options(
             <GDALResampleAlg>c_resampling, self.src_nodata,
             self.dst_nodata, src_dataset.count, dst_alpha,
-            src_alpha_band, warp_mem_limit, <const char **>c_warp_extras)
+            src_alpha_band, warp_mem_limit, <GDALDataType>dtypes.dtype_rev[self.working_dtype],
+            <const char **>c_warp_extras)
 
         if psWOptions == NULL:
             raise RuntimeError("Warp options are NULL")

rasterio/compat.py

@@ -12,7 +12,7 @@ if sys.version_info[0] >= 3:   # pragma: no cover
     import configparser
     from urllib.parse import urlparse
     from collections import UserDict
-    from collections.abc import Mapping
+    from collections.abc import Iterable, Mapping
     from inspect import getfullargspec as getargspec
 else:  # pragma: no cover
     string_types = basestring,
@@ -23,4 +23,4 @@ else:  # pragma: no cover
     from urlparse import urlparse
     from UserDict import UserDict
     from inspect import getargspec
-    from collections import Mapping
+    from collections import Iterable, Mapping

rasterio/crs.py

@@ -10,7 +10,6 @@ used.
 
 """
 
-import collections
 import json
 import pickle
 

rasterio/features.py

@@ -419,8 +419,8 @@ def geometry_window(dataset, shapes, pad_x=0, pad_y=0, north_up=True,
         right = min(dataset.shape[1], right)
         bottom = min(dataset.shape[0], bottom)
         # convert the bounds back to the CRS domain
-        left, top = (left, top) * dataset.transform
-        right, bottom = (right, bottom) * dataset.transform
+        left, top = dataset.transform * (left, top)
+        right, bottom = dataset.transform * (right, bottom)
 
     window = dataset.window(left, bottom, right, top)
     window_floored = window.round_offsets(op='floor', pixel_precision=pixel_precision)

rasterio/rio/stack.py

 """$ rio stack"""
 
 
-import collections
 import logging
 
 import click
 from cligj import format_opt
 
 import rasterio
-from rasterio.compat import zip_longest
+from rasterio.compat import Iterable, zip_longest
 from rasterio.rio import options
 from rasterio.rio.helpers import resolve_inout
 
@@ -104,7 +103,7 @@ def stack(ctx, files, output, driver, bidx, photometric, overwrite,
                             data = src.read(index)
                             dst.write(data, dst_idx)
                             dst_idx += 1
-                        elif isinstance(index, collections.Iterable):
+                        elif isinstance(index, Iterable):
                             data = src.read(index)
                             dst.write(data, range(dst_idx, dst_idx + len(index)))
                             dst_idx += len(index)

rasterio/transform.py

@@ -2,11 +2,12 @@
 
 from __future__ import division
 
-import collections
 import math
 
 from affine import Affine
 
+from rasterio.compat import Iterable
+
 
 IDENTITY = Affine.identity()
 GDAL_IDENTITY = IDENTITY.to_gdal()
@@ -153,10 +154,10 @@ def xy(transform, rows, cols, offset='center'):
 
     single_col = False
     single_row = False
-    if not isinstance(cols, collections.Iterable):
+    if not isinstance(cols, Iterable):
         cols = [cols]
         single_col = True
-    if not isinstance(rows, collections.Iterable):
+    if not isinstance(rows, Iterable):
         rows = [rows]
         single_row = True
 
@@ -221,10 +222,10 @@ def rowcol(transform, xs, ys, op=math.floor, precision=None):
 
     single_x = False
     single_y = False
-    if not isinstance(xs, collections.Iterable):
+    if not isinstance(xs, Iterable):
         xs = [xs]
         single_x = True
-    if not isinstance(ys, collections.Iterable):
+    if not isinstance(ys, Iterable):
         ys = [ys]
         single_y = True
 

rasterio/vrt.py

@@ -21,6 +21,57 @@ class WarpedVRT(WarpedVRTReaderBase, WindowMethodsMixin,
 
     This class is backed by an in-memory GDAL VRTWarpedDataset VRT file.
 
+    Parameters
+    ----------
+    src_dataset : dataset object
+        The warp source.
+    src_crs : CRS or str, optional
+        Overrides the coordinate reference system of `src_dataset`.
+    src_transfrom : Affine, optional
+        Overrides the transform of `src_dataset`.
+    src_nodata : float, optional
+        Overrides the nodata value of `src_dataset`, which is the
+        default.
+    crs : CRS or str, optional
+        The coordinate reference system at the end of the warp
+        operation.  Default: the crs of `src_dataset`. dst_crs is
+        a deprecated alias for this parameter.
+    transform : Affine, optional
+        The transform for the virtual dataset. Default: will be
+        computed from the attributes of `src_dataset`. dst_transform
+        is a deprecated alias for this parameter.
+    height, width: int, optional
+        The dimensions of the virtual dataset. Defaults: will be
+        computed from the attributes of `src_dataset`. dst_height
+        and dst_width are deprecated alias for these parameters.
+    nodata : float, optional
+        Nodata value for the virtual dataset. Default: the nodata
+        value of `src_dataset` or 0.0. dst_nodata is a deprecated
+        alias for this parameter.
+    resampling : Resampling, optional
+        Warp resampling algorithm. Default: `Resampling.nearest`.
+    tolerance : float, optional
+        The maximum error tolerance in input pixels when
+        approximating the warp transformation. Default: 0.125,
+        or one-eigth of a pixel.
+    src_alpha : int, optional
+        Index of a source band to use as an alpha band for warping.
+    add_alpha : bool, optional
+        Whether to add an alpha masking band to the virtual dataset.
+        Default: False. This option will cause deletion of the VRT
+        nodata value.
+    init_dest_nodata : bool, optional
+        Whether or not to initialize output to `nodata`. Default:
+        True.
+    warp_mem_limit : int, optional
+        The warp operation's memory limit in MB. The default (0)
+        means 64 MB with GDAL 2.2.
+    dtype : str, optional
+        The working data type for warp operation and output.
+    warp_extras : dict
+        GDAL extra warp options. See
+        http://www.gdal.org/structGDALWarpOptions.html.
+
     Attributes
     ----------
     src_dataset : dataset
@@ -39,6 +90,8 @@ class WarpedVRT(WarpedVRTReaderBase, WindowMethodsMixin,
         The nodata value used to initialize the destination; it will
         remain in all areas not covered by the reprojected source.
         Defaults to the value of src_nodata, or 0 (gdal default).
+    working_dtype : str, optional
+        The working data type for warp operation and output.
     warp_extras : dict
         GDAL extra warp options. See
         http://www.gdal.org/structGDALWarpOptions.html.

rasterio/windows.py

@@ -27,6 +27,7 @@ import attr
 from affine import Affine
 import numpy as np
 
+from rasterio.compat import Iterable
 from rasterio.errors import WindowError
 from rasterio.transform import rowcol, guard_transform
 
@@ -114,7 +115,7 @@ def iter_args(function):
     """
     @functools.wraps(function)
     def wrapper(*args, **kwargs):
-        if len(args) == 1 and isinstance(args[0], collections.Iterable):
+        if len(args) == 1 and isinstance(args[0], Iterable):
             return function(*args[0])
         else:
             return function(*args)

tests/conftest.py

@@ -27,8 +27,8 @@ if sys.version_info > (3,):
     reduce = functools.reduce
 
 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', 'data/RGBA.byte.tif']]
+    'data/RGB.byte.tif', 'data/float.tif', 'data/float32.tif',
+    'data/float_nan.tif', 'data/shade.tif', 'data/RGBA.byte.tif']]
 
 
 def pytest_cmdline_main(config):
@@ -601,6 +601,10 @@ requires_gdal22 = pytest.mark.skipif(
     not gdal_version.at_least('2.2'),
     reason="Requires GDAL 2.2.x")
 
+requires_gdal23 = pytest.mark.skipif(
+    not gdal_version.at_least('2.3'),
+    reason="Requires GDAL ~= 2.3")
+
 requires_gdal_lt_3 = pytest.mark.skipif(
     gdal_version.__lt__('3.0'),
     reason="Requires GDAL 1.x/2.x")

tests/test_blocks.py

@@ -11,6 +11,7 @@ import pytest
 import rasterio
 from rasterio import windows
 from rasterio.errors import RasterBlockError
+from rasterio.profiles import default_gtiff_profile
 
 from .conftest import requires_gdal2
 
@@ -51,6 +52,7 @@ class WindowTest(unittest.TestCase):
             rasterio.windows.evaluate(((None, -10), (None, -10)), 100, 90),
             windows.Window.from_slices((0, 90), (0, 80)))
 
+
 def test_window_index():
     idx = rasterio.windows.window_index(((0, 4), (1, 12)))
     assert len(idx) == 2
@@ -191,3 +193,29 @@ def test_block_window_tiff(path_rgb_byte_tif):
     with rasterio.open(path_rgb_byte_tif) as src:
         for (i, j), w in src.block_windows():
             assert src.block_window(1, i, j) == w
+
+
+@pytest.mark.parametrize("blocksize", [16, 32, 256, 1024])
+def test_block_windows_bigger_blocksize(tmpdir, blocksize):
+    """Ensure that block sizes greater than raster size are ok"""
+    tempfile = str(tmpdir.join("test.tif"))
+    profile = default_gtiff_profile.copy()
+    profile.update(height=16, width=16, count=1, blockxsize=blocksize, blockysize=blocksize)
+    with rasterio.open(tempfile, "w", **profile) as dst:
+        assert dst.is_tiled
+        for ij, window in dst.block_windows():
+            dst.write(np.ones((1, 1), dtype="uint8"), 1, window=window)
+
+    with rasterio.open(tempfile) as dst:
+        assert list(dst.block_windows()) == [((0, 0), windows.Window(0, 0, 16, 16))]
+        assert (dst.read(1) == 1).all()
+
+
+@pytest.mark.parametrize("blocksizes", [{"blockxsize": 33, "blockysize": 32}, {"blockxsize": 32, "blockysize": 33}])
+def test_odd_blocksize_error(tmpdir, blocksizes):
+    """For a tiled TIFF block sizes must be multiples of 16"""
+    tempfile = str(tmpdir.join("test.tif"))
+    profile = default_gtiff_profile.copy()
+    profile.update(height=64, width=64, count=1, **blocksizes)
+    with pytest.raises(RasterBlockError):
+        rasterio.open(tempfile, "w", **profile)

tests/test_colormap.py

@@ -5,17 +5,12 @@ import sys
 import rasterio
 
 
-logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
-
-
 def test_write_colormap_warn(tmpdir, recwarn):
     with rasterio.open('tests/data/shade.tif') as src:
         profile = src.meta
     tiffname = str(tmpdir.join('foo.tif'))
     with rasterio.open(tiffname, 'w', **profile) as dst:
         dst.write_colormap(1, {0: (255, 0, 0, 255), 255: (0, 0, 0, 0)})
-    w = recwarn.pop(UserWarning)
-    assert "The value will be ignored" in str(w.message)
 
 
 def test_write_colormap(tmpdir):

tests/test_dataset.py

@@ -47,23 +47,16 @@ def test_untiled_dataset_blocksize(tmpdir):
     """Blocksize is not relevant to untiled datasets (see #1689)"""
     tmpfile = str(tmpdir.join("test.tif"))
     with rasterio.open(
-            tmpfile, "w", driver="GTiff", count=1, height=13, width=13, dtype="uint8", crs="epsg:3857",
-            transform=Affine.identity(), blockxsize=256, blockysize=256) as dataset:
+            tmpfile, "w", driver="GTiff", count=1, height=13, width=23, dtype="uint8", crs="epsg:3857",
+            transform=Affine.identity(), blockxsize=64, blockysize=64) as dataset:
         pass
 
     with rasterio.open(tmpfile) as dataset:
         assert not dataset.profile["tiled"]
-        assert dataset.shape == (13, 13)
+        assert dataset.shape == (13, 23)
+        assert dataset.block_shapes == [(13, 23)]
 
 
-def test_tiled_dataset_blocksize_guard(tmpdir):
-    """Tiled datasets with dimensions less than blocksize are not permitted"""
-    tmpfile = str(tmpdir.join("test.tif"))
-    with pytest.raises(ValueError):
-        rasterio.open(
-            tmpfile, "w", driver="GTiff", count=1, height=13, width=13, dtype="uint8", crs="epsg:3857",
-            transform=Affine.identity(), tiled=True, blockxsize=256, blockysize=256)
-
 def test_dataset_readonly_attributes(path_rgb_byte_tif):
     """Attempts to set read-only attributes fail with DatasetAttributeError"""
     with pytest.raises(DatasetAttributeError):

tests/test_plot.py

@@ -58,6 +58,7 @@ def test_show_cmyk_interp(tmpdir):
         meta = src.meta
     meta['photometric'] = 'cmyk'
     meta['count'] = 4
+    del meta["nodata"]
     tiffname = str(tmpdir.join('foo.tif'))
     with rasterio.open(tiffname, 'w', **meta) as dst:
         assert dst.colorinterp == (