diff --git a/.travis.yml b/.travis.yml index 4eecae4e..649f6a74 100644 --- a/.travis.yml +++ b/.travis.yml @@ -48,11 +48,11 @@ jobs: - python: "3.8" env: GDALVERSION="master" - PROJVERSION="7.2.1" + PROJVERSION="8.0.1" allow_failures: - env: GDALVERSION="master" - PROJVERSION="7.2.1" + PROJVERSION="8.0.1" addons: apt: diff --git a/CHANGES.txt b/CHANGES.txt index 530bc7a7..51a42105 100644 --- a/CHANGES.txt +++ b/CHANGES.txt @@ -10,6 +10,28 @@ Changes environment variable. - Add version to CRS.to_wkt() for WKT2 support (#2122). +1.2.4 (2021-05-31) +------------------ + +- Eliminate unneeded marker for CLI nodata options to be ignored. We will stick + with None (#2191). +- Guard against use of gauss resampling when creating a WarpedVRT (#2190). +- Prevent segfaults when buffer and band indexes are mismatched in + io_multi_band and io_multi_mask (#2189). +- Change comparisons of nodata to IgnoreOption to accomodate changes in click + 8.0. +- Use Window constructor instead of from_slices in windows.union to allow a + proper union to be formed from windows extending outside a dataset (#2186). +- Read GDAL CInt16 data as np.complex64 and allow saving complex data to CInt16 + (#2185). +- Skip merge sources which do not overlap the destination. +- Allow unsafe casting in the merge tool, restoring behavior of version 1.2.0 + (#2179). +- Remove a workaround for an old Python 3.4 bug from the BoundingBox + implementation (#2172). +- Add setuptools as an explicit installation requirement. +>>>>>>> 1.2.4 + 1.2.3 (2021-04-26) ------------------ diff --git a/rasterio/_base.pyx b/rasterio/_base.pyx index b36c53ee..c11e3842 100644 --- a/rasterio/_base.pyx +++ b/rasterio/_base.pyx @@ -171,6 +171,7 @@ cdef _band_dtype(GDALRasterBandH band): else: return 'uint8' + return dtypes.dtype_fwd[gdal_dtype] @@ -549,7 +550,9 @@ cdef class DatasetBase: # to check that the return value is within the range of the # data type. If so, the band has a nodata value. If not, # there's no nodata value. - if (success == 0 or + if dtype not in dtypes.dtype_ranges: + pass + elif (success == 0 or val < dtypes.dtype_ranges[dtype][0] or val > dtypes.dtype_ranges[dtype][1]): val = None diff --git a/rasterio/_features.pyx b/rasterio/_features.pyx index 8cd5013a..8f6f50e7 100644 --- a/rasterio/_features.pyx +++ b/rasterio/_features.pyx @@ -9,6 +9,7 @@ import logging import numpy as np from rasterio import dtypes +from rasterio.dtypes import _getnpdtype from rasterio.enums import MergeAlg from rasterio._err cimport exc_wrap_int, exc_wrap_pointer @@ -62,12 +63,12 @@ def _shapes(image, mask, connectivity, transform): cdef ShapeIterator shape_iter = None cdef int fieldtp - is_float = np.dtype(image.dtype).kind == "f" + is_float = _getnpdtype(image.dtype).kind == "f" fieldtp = 2 if is_float else 0 valid_dtypes = ('int16', 'int32', 'uint8', 'uint16', 'float32') - if np.dtype(image.dtype).name not in valid_dtypes: + if _getnpdtype(image.dtype).name not in valid_dtypes: raise ValueError("image dtype must be one of: {0}".format( ', '.join(valid_dtypes))) @@ -89,7 +90,7 @@ def _shapes(image, mask, connectivity, transform): if mask.shape != image.shape: raise ValueError("Mask must have same shape as image") - if np.dtype(mask.dtype).name not in ('bool', 'uint8'): + if _getnpdtype(mask.dtype).name not in ('bool', 'uint8'): raise ValueError("Mask must be dtype rasterio.bool_ or " "rasterio.uint8") @@ -184,7 +185,7 @@ def _sieve(image, size, out, mask, connectivity): valid_dtypes = ('int16', 'int32', 'uint8', 'uint16') - if np.dtype(image.dtype).name not in valid_dtypes: + if _getnpdtype(image.dtype).name not in valid_dtypes: valid_types_str = ', '.join(('rasterio.{0}'.format(t) for t in valid_dtypes)) raise ValueError( @@ -203,7 +204,7 @@ def _sieve(image, size, out, mask, connectivity): if out.shape != image.shape: raise ValueError('out raster shape must be same as image shape') - if np.dtype(image.dtype).name != np.dtype(out.dtype).name: + if _getnpdtype(image.dtype).name != _getnpdtype(out.dtype).name: raise ValueError('out raster must match dtype of image') try: @@ -231,7 +232,7 @@ def _sieve(image, size, out, mask, connectivity): if mask.shape != image.shape: raise ValueError("Mask must have same shape as image") - if np.dtype(mask.dtype) not in ('bool', 'uint8'): + if _getnpdtype(mask.dtype) not in ('bool', 'uint8'): raise ValueError("Mask must be dtype rasterio.bool_ or " "rasterio.uint8") diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx index fdb08477..8581505f 100644 --- a/rasterio/_io.pyx +++ b/rasterio/_io.pyx @@ -26,7 +26,7 @@ from rasterio.errors import ( NotGeoreferencedWarning, NodataShadowWarning, WindowError, UnsupportedOperation, OverviewCreationError, RasterBlockError, InvalidArrayError ) -from rasterio.dtypes import is_ndarray +from rasterio.dtypes import is_ndarray, _is_complex_int, _getnpdtype from rasterio.sample import sample_gen from rasterio.transform import Affine from rasterio.path import parse_path, UnparsedPath @@ -300,7 +300,7 @@ cdef bint in_dtype_range(value, dtype): 105: np.iinfo, 117: np.iinfo} - key = np.dtype(dtype).kind + key = _getnpdtype(dtype).kind if np.isnan(value): return key in ('c', 'f', 99, 102) @@ -437,12 +437,14 @@ cdef class DatasetReaderBase(DatasetBase): check_dtypes = set() nodatavals = [] + # Check each index before processing 3D array for bidx in indexes: + if bidx not in self.indexes: raise IndexError("band index {} out of range (not in {})".format(bidx, self.indexes)) - idx = self.indexes.index(bidx) + idx = self.indexes.index(bidx) dtype = self.dtypes[idx] check_dtypes.add(dtype) @@ -450,16 +452,18 @@ cdef class DatasetReaderBase(DatasetBase): log.debug("Output nodata value read from file: %r", ndv) - if ndv is not None: - kind = np.dtype(dtype).kind - if chr(kind) in "iu": + if ndv is not None and not _is_complex_int(dtype): + kind = _getnpdtype(dtype).kind + + if kind in "iu": info = np.iinfo(dtype) dt_min, dt_max = info.min, info.max - elif chr(kind) in "cf": + elif kind in "cf": info = np.finfo(dtype) dt_min, dt_max = info.min, info.max else: dt_min, dt_max = False, True + if ndv < dt_min: ndv = dt_min elif ndv > dt_max: @@ -480,6 +484,9 @@ cdef class DatasetReaderBase(DatasetBase): if out_dtype is not None: dtype = out_dtype + # Ensure we have a numpy dtype. + dtype = _getnpdtype(dtype) + # Get the natural shape of the read window, boundless or not. # The window can have float values. In this case, we round up # when computing the shape. @@ -553,8 +560,7 @@ cdef class DatasetReaderBase(DatasetBase): log.debug("Jump straight to _read()") log.debug("Window: %r", window) - out = self._read(indexes, out, window, dtype, - resampling=resampling) + out = self._read(indexes, out, window, dtype, resampling=resampling) if masked or fill_value is not None: if all_valid: @@ -751,7 +757,7 @@ cdef class DatasetReaderBase(DatasetBase): out = np.zeros(out_shape, 'uint8') if out is not None: - if out.dtype != np.dtype(dtype): + if out.dtype != _getnpdtype(dtype): raise ValueError( "the out array's dtype '%s' does not match '%s'" % (out.dtype, dtype)) @@ -820,8 +826,7 @@ cdef class DatasetReaderBase(DatasetBase): return out - def _read(self, indexes, out, window, dtype, masks=False, - resampling=Resampling.nearest): + def _read(self, indexes, out, window, dtype, masks=False, resampling=Resampling.nearest): """Read raster bands as a multidimensional array If `indexes` is a list, the result is a 3D array, but @@ -1252,19 +1257,21 @@ cdef class DatasetWriterBase(DatasetReaderBase): try: width = int(width) height = int(height) - except: + except TypeError: raise TypeError("Integer width and height are required.") try: count = int(count) - except: + except TypeError: raise TypeError("Integer band count is required.") - try: - assert dtype is not None - _ = np.dtype(dtype) - except: - raise TypeError("A valid dtype is required.") - self._init_dtype = np.dtype(dtype).name + if _is_complex_int(dtype): + self._init_dtype = dtype + else: + try: + assert dtype is not None + self._init_dtype = _getnpdtype(dtype).name + except Exception: + raise TypeError("A valid dtype is required.") # Make and store a GDAL dataset handle. filename = path.name @@ -1343,7 +1350,9 @@ cdef class DatasetWriterBase(DatasetReaderBase): if nodata is not None: - if not in_dtype_range(nodata, dtype): + if _is_complex_int(dtype): + pass + elif not in_dtype_range(nodata, dtype): raise ValueError( "Given nodata value, %s, is beyond the valid " "range of its data type, %s." % ( @@ -1561,10 +1570,7 @@ cdef class DatasetWriterBase(DatasetReaderBase): else: # unique dtype; normal case dtype = check_dtypes.pop() - if arr is not None and arr.dtype != dtype: - raise ValueError( - "the array's dtype '%s' does not match " - "the file's dtype '%s'" % (arr.dtype, dtype)) + dtype = _getnpdtype(dtype) # Require C-continguous arrays (see #108). arr = np.require(arr, dtype=dtype, requirements='C') @@ -1769,7 +1775,7 @@ cdef class DatasetWriterBase(DatasetReaderBase): GDALFillRaster(mask, 255, 0) elif mask_array is False: GDALFillRaster(mask, 0, 0) - elif mask_array.dtype == np.bool: + elif mask_array.dtype == bool: array = 255 * mask_array.astype(np.uint8) io_band(mask, 1, xoff, yoff, width, height, array) else: @@ -2163,13 +2169,14 @@ cdef class BufferedDatasetWriterBase(DatasetWriterBase): count = int(count) except: raise TypeError("Integer band count is required.") + try: assert dtype is not None - _ = np.dtype(dtype) - except: + _ = _getnpdtype(dtype) + except Exception: raise TypeError("A valid dtype is required.") - self._init_dtype = np.dtype(dtype).name + self._init_dtype = _getnpdtype(dtype).name self.name = path.name self.mode = mode diff --git a/rasterio/_shim1.pyx b/rasterio/_shim1.pyx index 012c1c14..aee0cdcd 100644 --- a/rasterio/_shim1.pyx +++ b/rasterio/_shim1.pyx @@ -16,7 +16,7 @@ from rasterio.enums import Resampling cimport numpy as np from rasterio._err cimport exc_wrap_int, exc_wrap_pointer -from rasterio.errors import GDALOptionNotImplementedError +from rasterio.errors import GDALOptionNotImplementedError, DatasetIOShapeError cdef GDALDatasetH open_dataset( @@ -114,6 +114,9 @@ cdef int io_multi_band( cdef int xsize = width cdef int ysize = height + if len(indexes) != data.shape[0]: + raise DatasetIOShapeError("Dataset indexes and destination buffer are mismatched") + bandmap = CPLMalloc(count*sizeof(int)) for i in range(count): bandmap[i] = indexes[i] @@ -162,6 +165,9 @@ cdef int io_multi_mask( cdef int xsize = width cdef int ysize = height + if len(indexes) != data.shape[0]: + raise DatasetIOShapeError("Dataset indexes and destination buffer are mismatched") + for i in range(count): j = indexes[i] band = GDALGetRasterBand(hds, j) diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx index 1b7d2dc7..ffb309b5 100644 --- a/rasterio/_warp.pyx +++ b/rasterio/_warp.pyx @@ -43,6 +43,17 @@ from rasterio._features cimport GeomBuilder, OGRGeomBuilder log = logging.getLogger(__name__) +# Gauss (7) is not supported for warp +SUPPORTED_RESAMPLING = [r for r in Resampling if r.value < 7] +GDAL2_RESAMPLING = [r for r in Resampling if r.value > 7 and r.value <= 12] +if GDALVersion.runtime().at_least('2.0'): + SUPPORTED_RESAMPLING.extend(GDAL2_RESAMPLING) +# sum supported since GDAL 3.1 +if GDALVersion.runtime().at_least('3.1'): + SUPPORTED_RESAMPLING.append(Resampling.sum) +# rms supported since GDAL 3.3 +if GDALVersion.runtime().at_least('3.3'): + SUPPORTED_RESAMPLING.append(Resampling.rms) def recursive_round(val, precision): """Recursively round coordinates.""" @@ -774,6 +785,16 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase): if src_dataset.mode != "r": warnings.warn("Source dataset should be opened in read-only mode. Use of datasets opened in modes other than 'r' will be disallowed in a future version.", RasterioDeprecationWarning, stacklevel=2) + # Guard against invalid or unsupported resampling algorithms. + try: + if resampling == 7: + raise ValueError("Gauss resampling is not supported") + Resampling(resampling) + + except ValueError: + raise ValueError( + "resampling must be one of: {0}".format(", ".join(['Resampling.{0}'.format(r.name) for r in SUPPORTED_RESAMPLING]))) + self.mode = 'r' self.options = {} self._count = 0 diff --git a/rasterio/coords.py b/rasterio/coords.py index 1931cccc..cef25906 100644 --- a/rasterio/coords.py +++ b/rasterio/coords.py @@ -1,11 +1,9 @@ """Bounding box tuple, and disjoint operator.""" -from collections import namedtuple, OrderedDict +from collections import namedtuple -_BoundingBox = namedtuple('BoundingBox', ('left', 'bottom', 'right', 'top')) - - -class BoundingBox(_BoundingBox): +BoundingBox = namedtuple('BoundingBox', ('left', 'bottom', 'right', 'top')) +BoundingBox.__doc__ = \ """Bounding box named tuple, defining extent in cartesian coordinates. .. code:: @@ -24,10 +22,6 @@ class BoundingBox(_BoundingBox): Top coordinate """ - def _asdict(self): - return OrderedDict(zip(self._fields, self)) - - def disjoint_bounds(bounds1, bounds2): """Compare two bounds and determine if they are disjoint. @@ -43,6 +37,7 @@ def disjoint_bounds(bounds1, bounds2): boolean ``True`` if bounds are disjoint, ``False`` if bounds overlap + """ bounds1_north_up = bounds1[3] > bounds1[1] bounds2_north_up = bounds2[3] > bounds2[1] diff --git a/rasterio/dtypes.py b/rasterio/dtypes.py index 27615566..c73736f3 100644 --- a/rasterio/dtypes.py +++ b/rasterio/dtypes.py @@ -4,10 +4,6 @@ Since 0.13 we are not importing numpy here and data types are strings. Happily strings can be used throughout Numpy and so existing code will not break. -Within Rasterio, to test data types, we use Numpy's dtype() factory to -do something like this: - - if np.dtype(destination.dtype) == np.dtype(rasterio.uint8): ... """ bool_ = 'bool' @@ -23,26 +19,29 @@ complex_ = 'complex' complex64 = 'complex64' complex128 = 'complex128' -# Not supported: -# GDT_CInt16 = 8, GDT_CInt32 = 9, GDT_CFloat32 = 10, GDT_CFloat64 = 11 +complex_int16 = "complex_int16" dtype_fwd = { - 0: None, # GDT_Unknown - 1: ubyte, # GDT_Byte - 2: uint16, # GDT_UInt16 - 3: int16, # GDT_Int16 - 4: uint32, # GDT_UInt32 - 5: int32, # GDT_Int32 - 6: float32, # GDT_Float32 - 7: float64, # GDT_Float64 - 8: complex_, # GDT_CInt16 - 9: complex_, # GDT_CInt32 - 10: complex64, # GDT_CFloat32 - 11: complex128} # GDT_CFloat64 + 0: None, # GDT_Unknown + 1: ubyte, # GDT_Byte + 2: uint16, # GDT_UInt16 + 3: int16, # GDT_Int16 + 4: uint32, # GDT_UInt32 + 5: int32, # GDT_Int32 + 6: float32, # GDT_Float32 + 7: float64, # GDT_Float64 + 8: complex_int16, # GDT_CInt16 + 9: complex64, # GDT_CInt32 + 10: complex64, # GDT_CFloat32 + 11: complex128, # GDT_CFloat64 +} dtype_rev = dict((v, k) for k, v in dtype_fwd.items()) -dtype_rev['uint8'] = 1 -dtype_rev['int8'] = 1 + +dtype_rev["uint8"] = 1 +dtype_rev["int8"] = 1 +dtype_rev["complex"] = 11 +dtype_rev["complex_int16"] = 8 typename_fwd = { 0: 'Unknown', @@ -154,7 +153,7 @@ def can_cast_dtype(values, dtype): if not is_ndarray(values): values = np.array(values) - if values.dtype.name == np.dtype(dtype).name: + if values.dtype.name == _getnpdtype(dtype).name: return True elif values.dtype.kind == 'f': @@ -185,3 +184,15 @@ def validate_dtype(values, valid_dtypes): return (values.dtype.name in valid_dtypes or get_minimum_dtype(values) in valid_dtypes) + + +def _is_complex_int(dtype): + return isinstance(dtype, str) and dtype.startswith("complex_int") + + +def _getnpdtype(dtype): + import numpy as np + if _is_complex_int(dtype): + return np.dtype("complex64") + else: + return np.dtype(dtype) diff --git a/rasterio/errors.py b/rasterio/errors.py index 4f780edd..5452e17f 100644 --- a/rasterio/errors.py +++ b/rasterio/errors.py @@ -87,9 +87,9 @@ class GDALOptionNotImplementedError(RasterioError): by GDAL 1.x. """ + class GDALVersionError(RasterioError): - """Raised if the runtime version of GDAL does not meet the required - version of GDAL.""" + """Raised if the runtime version of GDAL does not meet the required version of GDAL.""" class WindowEvaluationError(ValueError): @@ -142,3 +142,7 @@ class TransformError(RasterioError): class WarpedVRTError(RasterioError): """Raised when WarpedVRT can't be initialized""" + + +class DatasetIOShapeError(RasterioError): + """Raised when data buffer shape is a mismatch when reading and writing""" diff --git a/rasterio/features.py b/rasterio/features.py index c896e8a0..b7a8b4f8 100644 --- a/rasterio/features.py +++ b/rasterio/features.py @@ -1,6 +1,5 @@ """Functions for working with features in a raster dataset.""" - import logging import math import os @@ -8,11 +7,11 @@ import warnings import numpy as np - import rasterio._loading + with rasterio._loading.add_gdal_dll_directories(): import rasterio - from rasterio.dtypes import validate_dtype, can_cast_dtype, get_minimum_dtype + from rasterio.dtypes import validate_dtype, can_cast_dtype, get_minimum_dtype, getnpdtype from rasterio.enums import MergeAlg from rasterio.env import ensure_env from rasterio.errors import ShapeSkipWarning @@ -279,7 +278,7 @@ def rasterize( if dtype is not None and not can_cast_dtype(default_value_array, dtype): raise ValueError(format_cast_error('default_vaue', dtype)) - if dtype is not None and np.dtype(dtype).name not in valid_dtypes: + if dtype is not None and getnpdtype(dtype).name not in valid_dtypes: raise ValueError(format_invalid_dtype('dtype')) valid_shapes = [] @@ -335,7 +334,7 @@ def rasterize( raise ValueError(format_cast_error('shape values', dtype)) if out is not None: - if np.dtype(out.dtype).name not in valid_dtypes: + if getnpdtype(out.dtype).name not in valid_dtypes: raise ValueError(format_invalid_dtype('out')) if not can_cast_dtype(shape_values, out.dtype): diff --git a/rasterio/merge.py b/rasterio/merge.py index 43503d0b..4d0c0811 100644 --- a/rasterio/merge.py +++ b/rasterio/merge.py @@ -9,8 +9,10 @@ import warnings import numpy as np import rasterio._loading + with rasterio._loading.add_gdal_dll_directories(): import rasterio + from rasterio.coords import disjoint_bounds from rasterio.enums import Resampling from rasterio import windows from rasterio.transform import Affine @@ -23,13 +25,13 @@ def copy_first(merged_data, new_data, merged_mask, new_mask, **kwargs): mask = np.empty_like(merged_mask, dtype="bool") np.logical_not(new_mask, out=mask) np.logical_and(merged_mask, mask, out=mask) - np.copyto(merged_data, new_data, where=mask) + np.copyto(merged_data, new_data, where=mask, casting="unsafe") def copy_last(merged_data, new_data, merged_mask, new_mask, **kwargs): mask = np.empty_like(merged_mask, dtype="bool") np.logical_not(new_mask, out=mask) - np.copyto(merged_data, new_data, where=mask) + np.copyto(merged_data, new_data, where=mask, casting="unsafe") def copy_min(merged_data, new_data, merged_mask, new_mask, **kwargs): @@ -39,7 +41,7 @@ def copy_min(merged_data, new_data, merged_mask, new_mask, **kwargs): np.minimum(merged_data, new_data, out=merged_data, where=mask) np.logical_not(new_mask, out=mask) np.logical_and(merged_mask, mask, out=mask) - np.copyto(merged_data, new_data, where=mask) + np.copyto(merged_data, new_data, where=mask, casting="unsafe") def copy_max(merged_data, new_data, merged_mask, new_mask, **kwargs): @@ -49,7 +51,7 @@ def copy_max(merged_data, new_data, merged_mask, new_mask, **kwargs): np.maximum(merged_data, new_data, out=merged_data, where=mask) np.logical_not(new_mask, out=mask) np.logical_and(merged_mask, mask, out=mask) - np.copyto(merged_data, new_data, where=mask) + np.copyto(merged_data, new_data, where=mask, casting="unsafe") MERGE_METHODS = { @@ -292,6 +294,10 @@ def merge( # This approach uses the maximum amount of memory to solve the # problem. Making it more efficient is a TODO. + if disjoint_bounds((dst_w, dst_s, dst_e, dst_n), src.bounds): + logger.debug("Skipping source: src=%r, window=%r", src) + continue + # 1. Compute spatial intersection of destination and source src_w, src_s, src_e, src_n = src.bounds @@ -304,7 +310,6 @@ def merge( src_window = windows.from_bounds( int_w, int_s, int_e, int_n, src.transform, precision=precision ) - logger.debug("Src %s window: %r", src.name, src_window) # 3. Compute the destination window dst_window = windows.from_bounds( diff --git a/rasterio/rio/edit_info.py b/rasterio/rio/edit_info.py index 64d8242e..a9d4a34c 100644 --- a/rasterio/rio/edit_info.py +++ b/rasterio/rio/edit_info.py @@ -190,13 +190,13 @@ def edit(ctx, input, bidx, nodata, unset_nodata, crs, unset_crs, transform, tags = allmd['tags'] colorinterp = allmd['colorinterp'] - if unset_nodata and nodata is not options.IgnoreOption: + if unset_nodata and nodata is not None: raise click.BadParameter( - "--unset-nodata and --nodata cannot be used together.") + "--unset-nodata and --nodata cannot be used together." + ) if unset_crs and crs: - raise click.BadParameter( - "--unset-crs and --crs cannot be used together.") + raise click.BadParameter("--unset-crs and --crs cannot be used together.") if unset_nodata: # Setting nodata to None will raise NotImplementedError @@ -207,17 +207,18 @@ def edit(ctx, input, bidx, nodata, unset_nodata, crs, unset_crs, transform, except NotImplementedError as exc: # pragma: no cover raise click.ClickException(str(exc)) - elif nodata is not options.IgnoreOption: + elif nodata is not None: dtype = dst.dtypes[0] if nodata is not None and not in_dtype_range(nodata, dtype): raise click.BadParameter( - "outside the range of the file's " - "data type (%s)." % dtype, - param=nodata, param_hint='nodata') + "outside the range of the file's data type (%s)." % dtype, + param=nodata, + param_hint="nodata", + ) dst.nodata = nodata if unset_crs: - dst.crs = None # CRS() + dst.crs = None elif crs: dst.crs = crs diff --git a/rasterio/rio/options.py b/rasterio/rio/options.py index c98c1a0d..aba4779c 100644 --- a/rasterio/rio/options.py +++ b/rasterio/rio/options.py @@ -51,26 +51,12 @@ import click import rasterio import rasterio.shutil -from rasterio.path import parse_path, ParsedPath, UnparsedPath +from rasterio.path import parse_path, UnparsedPath logger = logging.getLogger(__name__) -class IgnoreOptionMarker: - """A marker for an option that is to be ignored. - - For use in the case where `None` is a meaningful option value, - such as for nodata where `None` means "unset the nodata value." - """ - - def __repr__(self): - return 'IgnoreOption()' - - -IgnoreOption = IgnoreOptionMarker() - - def _cb_key_val(ctx, param, value): """ @@ -180,29 +166,36 @@ def like_handler(ctx, param, value): def nodata_handler(ctx, param, value): """Return a float or None""" - if value is None or value is IgnoreOption: - return value - elif value.lower() in ['null', 'nil', 'none', 'nada']: + if value is None or value.lower() in ["null", "nil", "none", "nada"]: return None else: try: return float(value) except (TypeError, ValueError): raise click.BadParameter( - "{!r} is not a number".format(value), - param=param, param_hint='nodata') + "{!r} is not a number".format(value), param=param, param_hint="nodata" + ) def edit_nodata_handler(ctx, param, value): """Get nodata value from a template file or command line. - Expected values are 'like', 'null', a numeric value, 'nan', or - IgnoreOption. Anything else should raise BadParameter. + Expected values are 'like', 'null', a numeric value, 'nan', or None. + + Returns + ------- + float or None + + Raises + ------ + click.BadParameter + """ - if value == 'like' or value is IgnoreOption: + if value == "like" or value is None: retval = from_like_context(ctx, param, value) if retval is not None: return retval + return nodata_handler(ctx, param, value) @@ -335,12 +328,20 @@ overwrite_opt = click.option( help="Always overwrite an existing output file.") nodata_opt = click.option( - '--nodata', callback=nodata_handler, default=None, - metavar='NUMBER|nan', help="Set a Nodata value.") + "--nodata", + callback=nodata_handler, + default=None, + metavar="NUMBER|nan", + help="Set a Nodata value.", +) edit_nodata_opt = click.option( - '--nodata', callback=edit_nodata_handler, default=IgnoreOption, - metavar='NUMBER|nan|null', help="Modify the Nodata value.") + "--nodata", + callback=edit_nodata_handler, + default=None, + metavar="NUMBER|nan|null", + help="Modify the Nodata value.", +) like_opt = click.option( '--like', diff --git a/rasterio/windows.py b/rasterio/windows.py index 4fe27c09..3cb47335 100644 --- a/rasterio/windows.py +++ b/rasterio/windows.py @@ -189,9 +189,14 @@ def union(*windows): Window """ stacked = np.dstack([toranges(w) for w in windows]) - return Window.from_slices( - (stacked[0, 0].min(), stacked[0, 1].max()), - (stacked[1, 0].min(), stacked[1, 1]. max())) + row_start, row_stop = stacked[0, 0].min(), stacked[0, 1].max() + col_start, col_stop = stacked[1, 0].min(), stacked[1, 1].max() + return Window( + col_off=col_start, + row_off=row_start, + width=col_stop - col_start, + height=row_stop - row_start, + ) @iter_args @@ -213,9 +218,14 @@ def intersection(*windows): raise WindowError("windows do not intersect") stacked = np.dstack([toranges(w) for w in windows]) - return Window.from_slices( - (stacked[0, 0].max(), stacked[0, 1].min()), - (stacked[1, 0].max(), stacked[1, 1]. min())) + row_start, row_stop = stacked[0, 0].max(), stacked[0, 1].min() + col_start, col_stop = stacked[1, 0].max(), stacked[1, 1].min() + return Window( + col_off=col_start, + row_off=row_start, + width=col_stop - col_start, + height=row_stop - row_start, + ) @iter_args diff --git a/setup.py b/setup.py index 0caa8669..70e964e7 100755 --- a/setup.py +++ b/setup.py @@ -255,7 +255,7 @@ inst_reqs = [ "affine", "attrs", "certifi", - "click>=4.0,<8", + "click>=4.0", "cligj>=0.5", "numpy", "snuggs>=1.4.1", diff --git a/tests/test_complex_dtypes.py b/tests/test_complex_dtypes.py index fc06fe9e..96a9837b 100644 --- a/tests/test_complex_dtypes.py +++ b/tests/test_complex_dtypes.py @@ -1,4 +1,5 @@ import logging +import subprocess import sys import uuid @@ -8,9 +9,6 @@ import pytest import rasterio -logging.basicConfig(stream=sys.stderr, level=logging.DEBUG) - - @pytest.fixture(scope='function') def tempfile(): """A temporary filename in the GDAL '/vsimem' filesystem""" @@ -23,10 +21,8 @@ def complex_image(height, width, dtype): [complex(x, x) for x in range(height * width)], dtype=dtype).reshape(height, width) -dtypes = ['complex', 'complex64', 'complex128'] - -@pytest.mark.parametrize("dtype", dtypes) +@pytest.mark.parametrize("dtype", ["complex", "complex64", "complex128"]) @pytest.mark.parametrize("height,width", [(20, 20)]) def test_read_array(tempfile, dtype, height, width): """_io functions read and write arrays correctly""" @@ -58,3 +54,45 @@ def test_complex_nodata(tmpdir): with rasterio.open(tempfile) as dst: assert dst.nodata == 0 + + +@pytest.mark.gdalbin +def test_complex_int16(tmpdir): + """A cint16 dataset can be created""" + import numpy as np + import rasterio + from rasterio.transform import Affine + + x = np.linspace(-4.0, 4.0, 240) + y = np.linspace(-3.0, 3.0, 180) + X, Y = np.meshgrid(x, y) + Z1 = np.ones_like(X) + 1j + + res = (x[-1] - x[0]) / 240.0 + transform1 = Affine.translation(x[0] - res / 2, y[-1] - res / 2) * Affine.scale( + res, -res + ) + + tempfile = str(tmpdir.join("test.tif")) + + with rasterio.open( + tempfile, + "w", + driver="GTiff", + height=Z1.shape[0], + width=Z1.shape[1], + nodata=0, + count=1, + dtype="complex_int16", + crs="+proj=latlong", + transform=transform1, + ) as dst: + dst.write(Z1, 1) + + assert "Type=CInt16" in subprocess.check_output(["gdalinfo", tempfile]).decode( + "utf-8" + ) + + with rasterio.open(tempfile) as dst: + data = dst.read() + assert data.dtype == np.complex64 diff --git a/tests/test_dtypes.py b/tests/test_dtypes.py index e954f536..e7e24f4a 100644 --- a/tests/test_dtypes.py +++ b/tests/test_dtypes.py @@ -3,10 +3,26 @@ import pytest import rasterio from rasterio import ( - ubyte, uint8, uint16, uint32, int16, int32, float32, float64, complex_) + ubyte, + uint8, + uint16, + uint32, + int16, + int32, + float32, + float64, + complex_, + complex_int16, +) from rasterio.dtypes import ( - _gdal_typename, is_ndarray, check_dtype, get_minimum_dtype, can_cast_dtype, - validate_dtype + _gdal_typename, + is_ndarray, + check_dtype, + get_minimum_dtype, + can_cast_dtype, + validate_dtype, + _is_complex_int, + _getnpdtype, ) @@ -28,10 +44,19 @@ def test_check_dtype_invalid(): assert not check_dtype('foo') -def test_gdal_name(): - assert _gdal_typename(ubyte) == 'Byte' - assert _gdal_typename(np.uint8) == 'Byte' - assert _gdal_typename(np.uint16) == 'UInt16' +@pytest.mark.parametrize( + ("dtype", "name"), + [ + (ubyte, "Byte"), + (np.uint8, "Byte"), + (np.uint16, "UInt16"), + ("uint8", "Byte"), + ("complex_int16", "CInt16"), + (complex_int16, "CInt16"), + ], +) +def test_gdal_name(dtype, name): + assert _gdal_typename(dtype) == name def test_get_minimum_dtype(): @@ -46,8 +71,8 @@ def test_get_minimum_dtype(): assert get_minimum_dtype(np.array([0, 1], dtype=np.uint)) == uint8 assert get_minimum_dtype(np.array([0, 1000], dtype=np.uint)) == uint16 assert get_minimum_dtype(np.array([0, 100000], dtype=np.uint)) == uint32 - assert get_minimum_dtype(np.array([-1, 0, 1], dtype=np.int)) == int16 - assert get_minimum_dtype(np.array([-1, 0, 100000], dtype=np.int)) == int32 + assert get_minimum_dtype(np.array([-1, 0, 1], dtype=int)) == int16 + assert get_minimum_dtype(np.array([-1, 0, 100000], dtype=int)) == int32 assert get_minimum_dtype(np.array([-1.5, 0, 1.5], dtype=np.float64)) == float32 @@ -89,3 +114,17 @@ def test_complex(tmpdir): arr2 = src.read(1) assert np.array_equal(arr1, arr2) + + +def test_is_complex_int(): + assert _is_complex_int("complex_int16") + + +def test_not_is_complex_int(): + assert not _is_complex_int("complex") + + +def test_get_npdtype(): + npdtype = _getnpdtype("complex_int16") + assert npdtype == np.complex64 + assert npdtype.kind == "c" diff --git a/tests/test_merge.py b/tests/test_merge.py index 7bc79139..e43ef908 100644 --- a/tests/test_merge.py +++ b/tests/test_merge.py @@ -59,3 +59,10 @@ def test_issue2163(): data = src.read() result, transform = merge([src]) assert numpy.allclose(data, result) + + +def test_unsafe_casting(): + """Demonstrate fix for issue 2179""" + with rasterio.open("tests/data/float_raster_with_nodata.tif") as src: + result, transform = merge([src], dtype="uint8", nodata=0.0) + assert not result.any() # this is why it's called "unsafe". diff --git a/tests/test_read.py b/tests/test_read.py index 10d87180..f1b2a1c8 100644 --- a/tests/test_read.py +++ b/tests/test_read.py @@ -1,23 +1,18 @@ from hashlib import md5 -import logging -import sys import unittest import numpy as np import pytest import rasterio - - -logging.basicConfig(stream=sys.stderr, level=logging.DEBUG) - +from rasterio.errors import DatasetIOShapeError # Find out if we've got HDF support (needed below). try: with rasterio.open('tests/data/no_band.h5') as s: pass has_hdf = True -except: +except Exception: has_hdf = False @@ -93,13 +88,7 @@ class ReaderContextTest(unittest.TestCase): def test_read_out_dtype_fail(self): with rasterio.open('tests/data/RGB.byte.tif') as s: a = np.zeros((718, 791), dtype=rasterio.float32) - try: - s.read(1, a) - except ValueError as e: - assert ("the array's dtype 'float32' does not match the " - "file's dtype") in str(e) - except: - assert "failed to catch exception" is False + s.read(1, a) def test_read_basic(self): with rasterio.open('tests/data/shade.tif') as s: @@ -315,7 +304,7 @@ def test_out_shape_exceptions(path_rgb_byte_tif): out_shape = (src.count, src.height, src.width) reader(out=out, out_shape=out_shape) - with pytest.raises(ValueError): + with pytest.raises(Exception): out_shape = (5, src.height, src.width) reader(1, out_shape=out_shape) @@ -325,3 +314,10 @@ def test_out_shape_implicit(path_rgb_byte_tif): with rasterio.open(path_rgb_byte_tif) as src: out = src.read(indexes=(1, 2), out_shape=src.shape) assert out.shape == (2,) + src.shape + + +def test_out_shape_no_segfault(path_rgb_byte_tif): + """Prevent segfault as reported in 2189""" + with rasterio.open(path_rgb_byte_tif) as src: + with pytest.raises(DatasetIOShapeError): + src.read(out_shape=(2, src.height, src.width)) diff --git a/tests/test_rio_edit_info.py b/tests/test_rio_edit_info.py index 09d54bd1..94b9f586 100644 --- a/tests/test_rio_edit_info.py +++ b/tests/test_rio_edit_info.py @@ -65,7 +65,7 @@ def test_delete_nodata(data, runner): """Delete a dataset's nodata value""" inputfile = str(data.join('RGB.byte.tif')) result = runner.invoke( - main_group, ['edit-info', inputfile, '--unset-nodata']) + main_group, ['edit-info', inputfile, '--unset-nodata'], catch_exceptions=False) assert result.exit_code == 0 diff --git a/tests/test_rio_options.py b/tests/test_rio_options.py index 29e4476b..a63a4b85 100644 --- a/tests/test_rio_options.py +++ b/tests/test_rio_options.py @@ -9,8 +9,13 @@ import pytest from rasterio.enums import ColorInterp from rasterio.rio.options import ( - IgnoreOption, bounds_handler, file_in_handler, like_handler, - edit_nodata_handler, nodata_handler, _cb_key_val) + bounds_handler, + file_in_handler, + like_handler, + edit_nodata_handler, + nodata_handler, + _cb_key_val, +) class MockContext: @@ -186,14 +191,14 @@ def test_edit_nodata_callback_like(data): def test_edit_nodata_callback_all_like(data): ctx = MockContext() - ctx.obj['like'] = {'nodata': 0.0} - ctx.obj['all_like'] = True - assert edit_nodata_handler(ctx, MockOption('nodata'), IgnoreOption) == 0.0 + ctx.obj["like"] = {"nodata": 0.0} + ctx.obj["all_like"] = True + assert edit_nodata_handler(ctx, MockOption("nodata"), None) == 0.0 def test_edit_nodata_callback_ignore(data): ctx = MockContext() - assert edit_nodata_handler(ctx, MockOption('nodata'), IgnoreOption) is IgnoreOption + assert edit_nodata_handler(ctx, MockOption("nodata"), None) is None def test_edit_nodata_callback_none(data): diff --git a/tests/test_warp.py b/tests/test_warp.py index d2735544..f31f1499 100644 --- a/tests/test_warp.py +++ b/tests/test_warp.py @@ -1,7 +1,6 @@ """rasterio.warp module tests""" import json -import sys import logging from affine import Affine @@ -15,7 +14,6 @@ from rasterio.crs import CRS from rasterio.enums import Resampling from rasterio.env import GDALVersion from rasterio.errors import ( - GDALBehaviorChangeException, CRSError, GDALVersionError, TransformError, @@ -1136,6 +1134,7 @@ def test_reproject_resampling(path_rgb_byte_tif, method): assert np.count_nonzero(out) in expected[method] + @pytest.mark.parametrize("test3d,count_nonzero", [(True, 1309625), (False, 437686)]) def test_reproject_array_interface(test3d, count_nonzero, path_rgb_byte_tif): class DataArray: @@ -1333,11 +1332,12 @@ def test_resample_default_invert_proj(method): assert out.mean() > 0 +@pytest.mark.xfail(reason="Projection extents have changed with PROJ 8") def test_target_aligned_pixels(): """Issue 853 has been resolved""" with rasterio.open("tests/data/world.rgb.tif") as src: source = src.read(1) - profile = src.profile.copy() + profile = src.profile dst_crs = "EPSG:3857" @@ -1348,7 +1348,7 @@ def test_target_aligned_pixels(): ) dst_affine, dst_width, dst_height = aligned_target( - dst_affine, dst_width, dst_height, 100000.0 + dst_affine, dst_width, dst_height, 10000.0 ) profile["height"] = dst_height @@ -1366,7 +1366,7 @@ def test_target_aligned_pixels(): resampling=Resampling.nearest, ) - # Check that there is no black borders + # Check that there are no black borders assert out[:, 0].all() assert out[:, -1].all() assert out[0, :].all() @@ -1694,26 +1694,6 @@ def test_issue_1446_b(): 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 = CRS.from_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() - - def test_reproject_init_dest_nodata(): """No pixels should transfer over""" crs = CRS.from_epsg(4326) @@ -1785,16 +1765,16 @@ def test_reproject_rpcs_with_transformer_options(caplog): transform = Affine(0.001953396267361111, 0.0, -124.00013888888888, 0.0, -0.001953396267361111, 50.000138888888884) with mem.open( driver="GTiff", - width=1024, - height=1024, + width=1024, + height=1024, count=1, transform=transform, - dtype='int16', - crs=crs + dtype="int16", + crs=crs, ) as dem: # we flush dem dataset before letting GDAL read from vsimem dem.write_band(1, 500 * np.ones((1024, 1024), dtype='int16')) - + out = np.zeros( (3, src.profile["width"], src.profile["height"]), dtype=np.uint8 ) @@ -1810,7 +1790,7 @@ def test_reproject_rpcs_with_transformer_options(caplog): resampling=Resampling.nearest, RPC_DEM=dem.name, - ) + ) caplog.set_level(logging.INFO) reproject( rasterio.band(src, src.indexes), @@ -1820,8 +1800,8 @@ def test_reproject_rpcs_with_transformer_options(caplog): dst_crs="EPSG:3857", resampling=Resampling.nearest, - ) - + ) + assert not out.all() assert not out2.all() assert "RPC_DEM" in caplog.text @@ -1853,10 +1833,10 @@ def test_warp_gcps_compute_dst_transform_automatically_array(): assert not out[:, -1, -1].any() assert not out[:, -1, 0].any() + def test_warp_gcps_compute_dst_transform_automatically_reader(tmpdir): """Ensure we don't raise an exception when gcps passed without dst_transform, for a source dataset""" tiffname = str(tmpdir.join('test.tif')) - arr = np.ones((3, 800, 800), dtype=np.uint8) * 255 src_gcps = [ GroundControlPoint(row=0, col=0, x=156113, y=2818720, z=0), GroundControlPoint(row=0, col=800, x=338353, y=2785790, z=0), @@ -1867,7 +1847,7 @@ def test_warp_gcps_compute_dst_transform_automatically_reader(tmpdir): with rasterio.open(tiffname, mode='w', height=800, width=800, count=3, dtype=np.uint8) as source: source.gcps = (src_gcps, CRS.from_epsg(32618)) - + with rasterio.open(tiffname) as source: reproject( rasterio.band(source, source.indexes), @@ -1875,13 +1855,14 @@ def test_warp_gcps_compute_dst_transform_automatically_reader(tmpdir): dst_crs="EPSG:32618", resampling=Resampling.nearest ) - + assert not out.all() assert not out[:, 0, 0].any() assert not out[:, 0, -1].any() assert not out[:, -1, -1].any() assert not out[:, -1, 0].any() + def test_reproject_rpcs_exact_transformer(caplog): """Reproject using rational polynomial coefficients and DEM, requiring that we don't try to make an approximate transformer. @@ -1892,16 +1873,16 @@ def test_reproject_rpcs_exact_transformer(caplog): transform = Affine(0.001953396267361111, 0.0, -124.00013888888888, 0.0, -0.001953396267361111, 50.000138888888884) with mem.open( driver="GTiff", - width=1024, - height=1024, + width=1024, + height=1024, count=1, transform=transform, - dtype='int16', - crs=crs + dtype="int16", + crs=crs, ) as dem: # we flush dem dataset before letting GDAL read from vsimem dem.write_band(1, 500 * np.ones((1024, 1024), dtype='int16')) - + out = np.zeros( (3, src.profile["width"], src.profile["height"]), dtype=np.uint8 ) @@ -1915,10 +1896,10 @@ def test_reproject_rpcs_exact_transformer(caplog): dst_crs="EPSG:3857", resampling=Resampling.nearest, RPC_DEM=dem.name, - ) + ) assert "Created exact transformer" in caplog.text - + def test_reproject_rpcs_approx_transformer(caplog): """Reproject using rational polynomial coefficients without a DEM, for which it's @@ -1937,6 +1918,6 @@ def test_reproject_rpcs_approx_transformer(caplog): rpcs=src_rpcs, dst_crs="EPSG:3857", resampling=Resampling.nearest, - ) + ) - assert "Created approximate transformer" in caplog.text \ No newline at end of file + assert "Created approximate transformer" in caplog.text diff --git a/tests/test_warpedvrt.py b/tests/test_warpedvrt.py index bf2f7c5b..799b1eff 100644 --- a/tests/test_warpedvrt.py +++ b/tests/test_warpedvrt.py @@ -608,3 +608,10 @@ def test_issue2086(): with rasterio.open("tests/data/white-gemini-iv.vrt") as src: with WarpedVRT(src, crs=DST_CRS) as vrt: assert vrt.shape == (1031, 1146) + + +def test_gauss_no(path_rgb_byte_tif): + """Guard against the issue reported in #2190""" + with rasterio.open(path_rgb_byte_tif) as src: + with pytest.raises(Exception): + WarpedVRT(src, resampling=Resampling.gauss) diff --git a/tests/test_windows.py b/tests/test_windows.py index cac6f1ce..932216a4 100644 --- a/tests/test_windows.py +++ b/tests/test_windows.py @@ -1,4 +1,3 @@ -from collections import namedtuple import logging import math import sys @@ -10,7 +9,7 @@ from hypothesis import given, assume, settings, HealthCheck from hypothesis.strategies import floats, integers import rasterio -from rasterio.errors import RasterioDeprecationWarning, WindowError +from rasterio.errors import WindowError from rasterio.windows import ( crop, from_bounds, bounds, transform, evaluate, window_index, shape, Window, intersect, intersection, get_data_window, union, @@ -631,3 +630,27 @@ def test_zero_height(sy): """Permit a zero height window""" transform = Affine.translation(0, 45.0) * Affine.scale(1.0, sy) assert from_bounds(0.0, 44.0, 1.0, 44.0, transform).height == 0 + + +def test_union_boundless_left(): + """Windows entirely to the left of a dataset form a proper union""" + uw = union( + Window(col_off=-10, row_off=0, width=2, height=2), + Window(col_off=-8.5, row_off=0, width=2.5, height=2), + ) + assert uw.col_off == -10 + assert uw.width == 4 + assert uw.height == 2 + assert uw.row_off == 0 + + +def test_union_boundless_above(): + """Windows entirely above a dataset form a proper union""" + uw = union( + Window(col_off=0, row_off=-10, width=2, height=2), + Window(col_off=0, row_off=-8.5, width=2, height=2.5), + ) + assert uw.row_off == -10 + assert uw.height == 4 + assert uw.width == 2 + assert uw.col_off == 0 diff --git a/tests/test_write.py b/tests/test_write.py index eb78ab2a..4edcaa20 100644 --- a/tests/test_write.py +++ b/tests/test_write.py @@ -107,8 +107,12 @@ def test_write_sbyte(tmpdir): with rasterio.open( name, 'w', driver='GTiff', width=100, height=100, count=1, - dtype=a.dtype) as s: - s.write(a, indexes=1) + dtype=a.dtype) as dst: + dst.write(a, indexes=1) + + with rasterio.open(name) as dst: + assert (dst.read() == -33).all() + info = subprocess.check_output(["gdalinfo", "-stats", name]).decode('utf-8') assert "Minimum=-33.000, Maximum=-33.000, Mean=-33.000, StdDev=0.000" in info assert 'SIGNEDBYTE' in info