Issue1744bis (#1768)

* No more right multplication with instances of Affine

That usage is deprecated.

* Add a `dtype` keyword argument to the WarpedVRT constructor

It functions like gdalwarp's `-wt` and resolves #1744.

* Add an xfail-ing test

* Add another test with WarpedVRT nodata set

CHANGES.txt

@@ -1,22 +1,28 @@
 Changes
 =======
 
-1.0.27 (TBD)
-------------
+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.
 - Tear down coordinate transformer at the end of _base._transform(), fixing
-  #1713.
-- Remove unavoidable warning about 4-channel colormap entries in
-  DatasetWriterBase.write_colormap().
-- Another category of Python deprecation warnings has been eliminated (#1742,
-  #1764).
+  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.
+  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.
+  issue.
 
 1.0.26 (2019-08-26)
 -------------------

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/_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.
@@ -315,7 +317,7 @@ def _reproject(
         if not in_dtype_range(dst_nodata, destination.dtype):
             raise ValueError("dst_nodata must be in valid range for "
                              "destination dtype")
-        
+
     def format_transform(in_transform):
         if not in_transform:
             return in_transform
@@ -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/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/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.

tests/conftest.py

@@ -27,8 +27,8 @@
     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):

tests/test_warpedvrt.py

@@ -399,3 +399,37 @@ def test_warpedvrt_float32_preserve(data):
     with rasterio.open("tests/data/float32.tif") as src:
         with WarpedVRT(src, src_crs="EPSG:4326") as vrt:
             assert src.dtypes == vrt.dtypes == ("float32",)
+
+
+def test_warpedvrt_float32_override(data):
+    """Override GDAL defaults for working data type"""
+    float32file = str(data.join("float32.tif"))
+    with rasterio.open(float32file, "r+") as dst:
+        dst.nodata = -3.4028230607370965e+38
+
+    with rasterio.open(float32file) as src:
+        with WarpedVRT(src, src_crs="EPSG:4326", dtype="float32") as vrt:
+            assert src.dtypes == vrt.dtypes == ("float32",)
+
+
+def test_warpedvrt_float32_overridei_nodata(data):
+    """Override GDAL defaults for working data type"""
+    float32file = str(data.join("float32.tif"))
+    with rasterio.open(float32file, "r+") as dst:
+        dst.nodata = -3.4028230607370965e+38
+
+    with rasterio.open(float32file) as src:
+        with WarpedVRT(src, src_crs="EPSG:4326", nodata=0.0001, dtype="float32") as vrt:
+            assert src.dtypes == vrt.dtypes == ("float32",)
+
+
+@pytest.mark.xfail(reason="GDAL's output defaults to float64")
+def test_warpedvrt_issue1744(data):
+    """Reproduce the bug reported in 1744"""
+    float32file = str(data.join("float32.tif"))
+    with rasterio.open(float32file, "r+") as dst:
+        dst.nodata = -3.4028230607370965e+38
+
+    with rasterio.open(float32file) as src:
+        with WarpedVRT(src, src_crs="EPSG:4326") as vrt:
+            assert src.dtypes == vrt.dtypes == ("float32",)