CHANGES.txt

 Changes
 =======
 
+1.0.6 (2018-09-24)
+------------------
+
+Bug fixes:
+
+- If the build_overviews method of a dataset is passed a list of factors that
+  specify more than one 1x1 pixel overview (#1333), Rasterio raises an
+  exception.
+- Calling calculate_default_transform for large extents should no longer result
+  in the out of memory error reported in #1131. The rio-warp command should
+  also now run more quickly and with a smaller memory footprint.
+- We have a more general fix for the problem of filling the empty regions of
+  boundless reads (#1471).
+
 1.0.5 (2018-09-19)
 ------------------
 
@@ -8,7 +22,7 @@ Bug fixes:
 
 - The fill value for boundless reads was ignored in Rasterio versions 1-1.0.4
   but now applies (#1471).
-- An invalid shortcut has been eliminated a Rasterio now prroduces a proper
+- An invalid shortcut has been eliminated and Rasterio now produces a proper
   mask in the boundless masked read case (#1449).
 - Loss of a row or column in geometry_window() and mask() has been fixed
   (#1472).

debian/changelog

+rasterio (1.0.6-1) unstable; urgency=medium
+
+  * Team upload.
+  * New upstream release.
+
+ -- Bas Couwenberg <sebastic@debian.org>  Tue, 25 Sep 2018 07:03:48 +0200
+
 rasterio (1.0.5-1) unstable; urgency=medium
 
   * Team upload.

docs/topics/plotting.rst

@@ -12,7 +12,6 @@ two dimensional data can be accomplished directly with ``pyplot``.
     >>> src = rasterio.open("tests/data/RGB.byte.tif")
     >>> pyplot.imshow(src.read(1), cmap='pink')
     <matplotlib.image.AxesImage object at 0x...>
-    >>> pyplot.show = lambda : None  # prevents showing during doctests
     >>> pyplot.show()
 
 

rasterio/__init__.py

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

rasterio/_io.pyx

@@ -6,6 +6,7 @@ from __future__ import absolute_import
 include "directives.pxi"
 include "gdal.pxi"
 
+from collections import Counter
 import logging
 import sys
 import uuid
@@ -24,7 +25,7 @@ from rasterio.enums import ColorInterp, MaskFlags, Resampling
 from rasterio.errors import (
     CRSError, DriverRegistrationError, RasterioIOError,
     NotGeoreferencedWarning, NodataShadowWarning, WindowError,
-    UnsupportedOperation
+    UnsupportedOperation, OverviewCreationError
 )
 from rasterio.sample import sample_gen
 from rasterio.transform import Affine
@@ -235,14 +236,7 @@ cdef class DatasetReaderBase(DatasetBase):
 
             log.debug("Output nodata value read from file: %r", ndv)
 
-            # Change given nodatavals to the closest value that
-            # can be represented by this band's data type to
-            # match GDAL's strategy.
-            if fill_value:
-                ndv = fill_value
-                log.debug("Output nodata value set from fill value")
-
-            elif ndv is not None:
+            if ndv is not None:
                 if np.dtype(dtype).kind in ('i', 'u'):
                     info = np.iinfo(dtype)
                     dt_min, dt_max = info.min, info.max
@@ -321,12 +315,6 @@ cdef class DatasetReaderBase(DatasetBase):
             else:
                 out = np.empty(out_shape, dtype=dtype)
 
-            for i, (ndv, arr) in enumerate(zip(
-                    nodatavals, out if len(out.shape) == 3 else [out])):
-
-                if ndv is not None:
-                    arr.fill(ndv)
-
         # Masking
         # -------
         #
@@ -337,7 +325,6 @@ cdef class DatasetReaderBase(DatasetBase):
         # read_masks(), invert them and use them in constructing masked
         # arrays.
 
-        if masked:
         enums = self.mask_flag_enums
         all_valid = all([MaskFlags.all_valid in flags for flags in enums])
         log.debug("all_valid: %s", all_valid)
@@ -353,7 +340,7 @@ cdef class DatasetReaderBase(DatasetBase):
             out = self._read(indexes, out, window, dtype,
                              resampling=resampling)
 
-            if masked:
+            if masked or fill_value is not None:
                 if all_valid:
                     mask = np.ma.nomask
                 else:
@@ -365,17 +352,23 @@ cdef class DatasetReaderBase(DatasetBase):
                 kwds = {'mask': mask}
                 # Set a fill value only if the read bands share a
                 # single nodata value.
-                if len(set(nodatavals)) == 1:
+                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 not masked:
+                    out = out.filled(fill_value)
+
         # If this is a boundless read we will create an in-memory VRT
         # in order to use GDAL's windowing and compositing logic.
         else:
 
             vrt_doc = _boundless_vrt_doc(
-                self, nodata=ndv, hidenodata=bool(fill_value), width=max(self.width, window.width) + 1,
+                self, nodata=ndv, width=max(self.width, window.width) + 1,
                 height=max(self.height, window.height) + 1,
                 transform=self.window_transform(window)).decode('ascii')
 
@@ -389,7 +382,7 @@ cdef class DatasetReaderBase(DatasetBase):
                     indexes, out, Window(0, 0, window.width, window.height),
                     None, resampling=resampling)
 
-                if masked:
+                if masked or fill_value is not None:
                     mask = np.zeros(out.shape, 'uint8')
                     mask = ~vrt._read(
                         indexes, mask, Window(0, 0, window.width, window.height), None, masks=True).astype('bool')
@@ -397,11 +390,17 @@ cdef class DatasetReaderBase(DatasetBase):
                     kwds = {'mask': mask}
                     # Set a fill value only if the read bands share a
                     # single nodata value.
-                    if len(set(nodatavals)) == 1:
+                    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 not masked:
+                        out = out.filled(fill_value)
+
         if return2d:
             out.shape = out.shape[1:]
 
@@ -1551,6 +1550,11 @@ cdef class DatasetWriterBase(DatasetReaderBase):
                     ['Resampling.{0}'.format(Resampling(k).name) for k in
                      resampling_map.keys()])))
 
+        # Check factors
+        ovr_shapes = Counter([(int((self.height + f - 1) / f), int((self.width + f - 1) / f)) for f in factors])
+        if ovr_shapes[(1, 1)] > 1:
+            raise OverviewCreationError("Too many overviews levels of 1x1 dimension were requested")
+
         # Allocate arrays.
         if factors:
             factors_c = <int *>CPLMalloc(len(factors)*sizeof(int))

rasterio/_warp.pyx

@@ -6,10 +6,12 @@ include "gdal.pxi"
 import logging
 import uuid
 import warnings
+import xml.etree.ElementTree as ET
 
 from affine import identity
 import numpy as np
 
+import rasterio
 from rasterio._base import gdal_version
 from rasterio._err import (
     CPLE_IllegalArgError, CPLE_NotSupportedError,
@@ -20,6 +22,7 @@ from rasterio.enums import Resampling, MaskFlags, ColorInterp
 from rasterio.env import GDALVersion
 from rasterio.crs import CRS
 from rasterio.errors import (
+    GDALOptionNotImplementedError,
     DriverRegistrationError, CRSError, RasterioIOError,
     RasterioDeprecationWarning, WarpOptionsError)
 from rasterio.transform import Affine, from_bounds, guard_transform, tastes_like_gdal
@@ -31,7 +34,7 @@ from rasterio._err cimport exc_wrap_pointer, exc_wrap_int
 from rasterio._io cimport (
     DatasetReaderBase, InMemoryRaster, in_dtype_range, io_auto)
 from rasterio._features cimport GeomBuilder, OGRGeomBuilder
-from rasterio._shim cimport delete_nodata_value
+from rasterio._shim cimport delete_nodata_value, open_dataset
 
 
 log = logging.getLogger(__name__)
@@ -520,7 +523,7 @@ def _calculate_default_transform(src_crs, dst_crs, width, height,
     cdef double geotransform[6]
     cdef OGRSpatialReferenceH osr = NULL
     cdef char *wkt = NULL
-    cdef InMemoryRaster temp = None
+    cdef GDALDatasetH hds = NULL
 
     extent[:] = [0.0, 0.0, 0.0, 0.0]
     geotransform[:] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
@@ -540,15 +543,25 @@ def _calculate_default_transform(src_crs, dst_crs, width, height,
     OSRExportToWkt(osr, &wkt)
     _safe_osr_release(osr)
 
-    with InMemoryRaster(width=width, height=height, transform=transform,
-                        gcps=gcps, crs=src_crs) as temp:
+    if isinstance(src_crs, str):
+        src_crs = CRS.from_string(src_crs)
+    elif isinstance(src_crs, dict):
+        src_crs = CRS(**src_crs)
+
+    vrt_doc = _suggested_proxy_vrt_doc(width, height, transform=transform, crs=src_crs, gcps=gcps).decode('ascii')
+
     try:
+        try:
+            hds = open_dataset(vrt_doc, 0x00 | 0x02 | 0x04, ['VRT'], {}, None)
+        except GDALOptionNotImplementedError:
+            hds = open_dataset(vrt_doc, 0x00 | 0x02 | 0x04, None, None, None)
+
         hTransformArg = exc_wrap_pointer(
             GDALCreateGenImgProjTransformer(
-                    temp._hds, NULL, NULL, wkt, 1, 1000.0,0))
+                hds, NULL, NULL, wkt, 1, 1000.0,0))
         exc_wrap_int(
             GDALSuggestedWarpOutput2(
-                    temp._hds, GDALGenImgProjTransform, hTransformArg,
+                hds, GDALGenImgProjTransform, hTransformArg,
                 geotransform, &npixels, &nlines, extent, 0))
 
         log.debug("Created transformer and warp output.")
@@ -922,3 +935,32 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
             raise ValueError("WarpedVRT does not permit boundless reads")
         else:
             return super(WarpedVRTReaderBase, self).read_masks(indexes=indexes, out=out, window=window, out_shape=out_shape, resampling=resampling)
+
+
+def _suggested_proxy_vrt_doc(width, height, transform=None, crs=None, gcps=None):
+    """Make a VRT XML document to serve _calculate_default_transform."""
+    vrtdataset = ET.Element('VRTDataset')
+    vrtdataset.attrib['rasterYSize'] = str(height)
+    vrtdataset.attrib['rasterXSize'] = str(width)
+    vrtrasterband = ET.SubElement(vrtdataset, 'VRTRasterBand')
+
+    srs = ET.SubElement(vrtdataset, 'SRS')
+    srs.text = crs.wkt if crs else ""
+
+    if gcps:
+        gcplist = ET.SubElement(vrtdataset, 'GCPList')
+        gcplist.attrib['Projection'] = crs.wkt if crs else ""
+        for point in gcps:
+            gcp = ET.SubElement(gcplist, 'GCP')
+            gcp.attrib['Id'] = str(point.id)
+            gcp.attrib['Info'] = str(point.info)
+            gcp.attrib['Pixel'] = str(point.col)
+            gcp.attrib['Line'] = str(point.row)
+            gcp.attrib['X'] = str(point.x)
+            gcp.attrib['Y'] = str(point.y)
+            gcp.attrib['Z'] = str(point.z)
+    else:
+        geotransform = ET.SubElement(vrtdataset, 'GeoTransform')
+        geotransform.text = ','.join([str(v) for v in transform.to_gdal()])
+
+    return ET.tostring(vrtdataset)

rasterio/errors.py

@@ -98,3 +98,7 @@ class WarpOptionsError(RasterioError):
 
 class UnsupportedOperation(RasterioError):
     """Raised when reading from a file opened in 'w' mode"""
+
+
+class OverviewCreationError(RasterioError):
+    """Raised when creation of an overview fails"""

tests/test_boundless_read.py

@@ -80,3 +80,18 @@ def test_boundless_fill_value():
     assert (filled[:, -1] == 5).all()
     assert (filled[0, :] == 5).all()
     assert (filled[-1, :] == 5).all()
+
+
+def test_boundless_fill_value_overview_masks():
+    """Confirm a more general resolution to issue #1471"""
+    with rasterio.open("tests/data/cogeo.tif") as src:
+        data = src.read(1, boundless=True, window=Window(-300, -335, 1000, 1000), fill_value=5, out_shape=(512, 512))
+    assert (data[:, 0] == 5).all()
+
+
+def test_boundless_masked_fill_value_overview_masks():
+    """Confirm a more general resolution to issue #1471"""
+    with rasterio.open("tests/data/cogeo.tif") as src:
+        data = src.read(1, masked=True, boundless=True, window=Window(-300, -335, 1000, 1000), fill_value=5, out_shape=(512, 512))
+    assert data.fill_value == 5
+    assert data.mask[:, 0].all()

tests/test_overviews.py

@@ -9,6 +9,8 @@ import pytest
 import rasterio
 from rasterio.enums import Resampling
 from rasterio.env import GDALVersion
+from rasterio.errors import OverviewCreationError
+
 
 gdal_version = GDALVersion()
 
@@ -79,3 +81,12 @@ def test_test_unsupported_algo(data):
         with rasterio.open(inputfile, 'r+') as src:
             overview_factors = [2, 4]
             src.build_overviews(overview_factors, resampling=Resampling.q1)
+
+
+def test_issue1333(data):
+    """Fail if asked to create more than one 1x1 overview"""
+    inputfile = str(data.join('RGB.byte.tif'))
+    with pytest.raises(OverviewCreationError):
+        with rasterio.open(inputfile, 'r+') as src:
+            overview_factors = [1024, 2048]
+            src.build_overviews(overview_factors, resampling=Resampling.average)

tests/test_warp_transform.py

@@ -4,6 +4,7 @@ import pytest
 
 from rasterio._warp import _calculate_default_transform
 from rasterio.control import GroundControlPoint
+from rasterio.crs import CRS
 from rasterio.errors import CRSError
 from rasterio.transform import from_bounds
 from rasterio.warp import calculate_default_transform, transform_bounds
@@ -60,7 +61,7 @@ def test_identity():
     assert res_width == width
     assert res_height == height
     for res, exp in zip(res_transform, transform):
-        assert round(res, 7) == round(exp, 7)
+        assert round(res, 3) == round(exp, 3)
 
 
 def test_identity_gcps():
@@ -88,7 +89,7 @@ def test_identity_gcps():
     assert res_width == width
     assert res_height == height
     for res, exp in zip(res_transform, transform):
-        assert round(res, 7) == round(exp, 7)
+        assert round(res, 3) == round(exp, 3)
 
 
 def test_transform_bounds():
@@ -173,3 +174,9 @@ def test_transform_bounds_identity():
     """Confirm fix of #1411"""
     bounds = (12978395.906596646, 146759.09430753812, 12983287.876406897, 151651.06411778927)
     assert transform_bounds("+init=epsg:3857", "+init=epsg:3857", *bounds) == bounds
+
+
+def test_issue1131():
+    """Confirm that we don't run out of memory"""
+    transform, w, h = calculate_default_transform(CRS.from_epsg(4326), CRS.from_epsg(3857), 455880, 454450, 13.0460235139, 42.6925552354, 13.2511695428, 42.8970561511)
+    assert (w, h) == (381595, 518398)