diff --git a/AUTHORS.txt b/AUTHORS.txt index f75c6ac6..3e35b11f 100644 --- a/AUTHORS.txt +++ b/AUTHORS.txt @@ -1,32 +1,34 @@ Authors ======= -Sean Gillies -Brendan Ward -Matthew Perry -Kevin Wurster -James McBride -Amit Kapadia -Kelsey Jordahl -Even Rouault -Mike Toews -Maxim Dubinin -Ryan Grout -Asger Petersen -Joshua Arnott +Aldo Culquicondor Alessandro Amici -Johan Van de Wauw +Alexander +Amit Kapadia +AsgerPetersen +Bas Couwenberg +Brendan Ward +Etienne B. Racine +Even Rouault +Jacques Tardie +James McBride +James Seppi Jeffrey Gerard +Johan Van de Wauw +Joshua Arnott +Kelsey Jordahl +Kevin Wurster +Martijn Visser +Matt Savoie +Matthew Perry +Maxim Dubinin +Mike Toews Nat Wilson Patrick Young Robin Wilson -Jacques Tardie -Etienne B. Racine -Bas Couwenberg +Ryan Grout +Sean Gillies Trevor R.H. Clarke cgohlke -Martijn Visser -Aldo Culquicondor -James Seppi See also https://github.com/mapbox/rasterio/graphs/contributors. diff --git a/CHANGES.txt b/CHANGES.txt index 66f00187..74411eef 100644 --- a/CHANGES.txt +++ b/CHANGES.txt @@ -1,6 +1,41 @@ Changes ======= +0.32.0 (2016-03-22) +------------------- +- Bug fix: geometry factories and warp operations are properly deallocated + in normal and error situations (#494, #568). +- Bug fix: a code block in rio-merge's help has been better formatted (#535). +- Bug fix: the rasterio.vfs module is imported in __init__.py to assist + cx_Freeze (#536). +- Bug fix: old usage of `read_band()` has been replaced by `read()` throughout + the docs (#537). +- Bug fix: accidental overwriting of existing files is now prevented by the + `resolve_inout()` function in `rasterio.rio.helpers`. Commands that take + one or more input files plus an output file should use this helper and force + overwrite either by using a `--force-overwrite` option or by using the + `-o output` option, which implicitly forces overwriting (#539, #540). +- Bug fix: missing support for NaN nodata value in rio-warp added (#542, #544). +- Bug fix: missing documentation of `rasterize()`'s `fill` parameter added + (#543). +- Bug fix: raster dataset bounds are densified before transforming so that + the projected output of rio-bounds is correct (#556, #557). +- Bug fix: add 'line' to the `Interleaving` enum (#560). +- Bug fix: convert `matplotlib` import errors to a `RuntimeWarning` (#562). +- Bug fix: deallocate CPL strings in error cases (#573). +- Bug fix: non-invertable affine transforms are prevented using + `__future__.division` *#580). +- Bug fix: rio-warp clips output regions to the limits of the destination + CRS unless disabled with `--no-check-invert-proj` (#597). +- New feature: the functionality previously available only in rio-mask is now + available as `rasterio.tools.mask.mask()` (#552). +- New feature: raster bounds are used to label axes in `rasterio.tool.show()` + (#553). +- New feature: GDAL's suggested warp bounds algorithm is wrapped and exposed + for use in `warp()` and rio-warp (#574). +- Breaking change: align rio-warp's `--bounds` option with rio-merge's: these + are in destination CRS units (#541, #545). + 0.31.0 (2015-12-18) ------------------- - Warn when rasters have no georeferencing and when the default identity diff --git a/docs/cli.rst b/docs/cli.rst index 7913ee8e..34f31ca7 100644 --- a/docs/cli.rst +++ b/docs/cli.rst @@ -329,7 +329,7 @@ option. .. code-block:: console - $ rio info tests/data/RGB.byte.tif --indent 2 + $ rio info tests/data/RGB.byte.tif --indent 2 --verbose { "count": 3, "crs": "EPSG:32618", diff --git a/rasterio/__init__.py b/rasterio/__init__.py index 35727270..b3a4ec63 100644 --- a/rasterio/__init__.py +++ b/rasterio/__init__.py @@ -23,7 +23,7 @@ from rasterio import _err, coords, enums, vfs __all__ = [ 'band', 'open', 'drivers', 'copy', 'pad'] -__version__ = "0.31.0" +__version__ = "0.32.0" log = logging.getLogger('rasterio') class NullHandler(logging.Handler): diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx index d4ada0fb..1b5b2e7d 100644 --- a/rasterio/_io.pyx +++ b/rasterio/_io.pyx @@ -20,7 +20,7 @@ from rasterio._drivers import driver_count, GDALEnv from rasterio._err import cpl_errs, GDALError from rasterio import dtypes from rasterio.coords import BoundingBox -from rasterio.errors import RasterioDriverRegistrationError +from rasterio.errors import DriverRegistrationError from rasterio.five import text_type, string_types from rasterio.transform import Affine from rasterio.enums import ColorInterp, MaskFlags, Resampling @@ -1851,7 +1851,7 @@ cdef class InMemoryRaster: cdef void *memdriver = _gdal.GDALGetDriverByName("MEM") if memdriver == NULL: - raise RasterioDriverRegistrationError( + raise DriverRegistrationError( "MEM driver is not registered.") self.dataset = _gdal.GDALCreate( diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx index 8965dbd3..e2ad9402 100644 --- a/rasterio/_warp.pyx +++ b/rasterio/_warp.pyx @@ -10,7 +10,7 @@ from rasterio cimport _base, _gdal, _ogr, _io, _features from rasterio import dtypes from rasterio._err import cpl_errs, GDALError from rasterio._io cimport InMemoryRaster -from rasterio.errors import RasterioDriverRegistrationError +from rasterio.errors import DriverRegistrationError from rasterio.transform import Affine, from_bounds @@ -269,7 +269,7 @@ def _reproject( hrdriver = _gdal.GDALGetDriverByName("MEM") if hrdriver == NULL: - raise RasterioDriverRegistrationError( + raise DriverRegistrationError( "'MEM' driver not found. Check that this call is contained " "in a `with rasterio.drivers()` or `with rasterio.open()` " "block.") @@ -317,7 +317,7 @@ def _reproject( hrdriver = _gdal.GDALGetDriverByName("MEM") if hrdriver == NULL: - raise RasterioDriverRegistrationError( + raise DriverRegistrationError( "'MEM' driver not found. Check that this call is contained " "in a `with rasterio.drivers()` or `with rasterio.open()` " "block.") diff --git a/rasterio/errors.py b/rasterio/errors.py index 4cfbddc3..3a4f34e7 100644 --- a/rasterio/errors.py +++ b/rasterio/errors.py @@ -7,7 +7,7 @@ class RasterioIOError(IOError): """A failure to open a dataset using the presently registered drivers.""" -class RasterioDriverRegistrationError(ValueError): +class DriverRegistrationError(ValueError): """To be raised when, eg, _gdal.GDALGetDriverByName("MEM") returns NULL.""" diff --git a/rasterio/rio/warp.py b/rasterio/rio/warp.py index 09f01734..b65f8e87 100644 --- a/rasterio/rio/warp.py +++ b/rasterio/rio/warp.py @@ -68,14 +68,16 @@ def x_dst_bounds_handler(ctx, param, value): @click.option('--resampling', type=click.Choice([r.name for r in Resampling]), default='nearest', help="Resampling method.", show_default=True) -@click.option('--threads', type=int, default=2, +@click.option('--threads', type=int, default=1, help='Number of processing threads.') +@click.option('--check-invert-proj', type=bool, default=True, + help='Constrain output to valid coordinate region in dst-crs') @options.force_overwrite_opt @options.creation_options @click.pass_context def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds, - x_dst_bounds, bounds, res, resampling, threads, force_overwrite, - creation_options): + x_dst_bounds, bounds, res, resampling, threads, check_invert_proj, + force_overwrite, creation_options): """ Warp a raster dataset. @@ -132,7 +134,8 @@ def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds, # Expand one value to two if needed res = (res[0], res[0]) if len(res) == 1 else res - with rasterio.drivers(CPL_DEBUG=verbosity > 2): + with rasterio.drivers(CPL_DEBUG=verbosity > 2, + CHECK_WITH_INVERT_PROJ=check_invert_proj): with rasterio.open(files[0]) as src: l, b, r, t = src.bounds out_kwargs = src.meta.copy() diff --git a/rasterio/warp.py b/rasterio/warp.py index c75dc0f0..81bfc9a1 100644 --- a/rasterio/warp.py +++ b/rasterio/warp.py @@ -1,6 +1,7 @@ """Raster warping and reprojection""" from __future__ import absolute_import +from __future__ import division from math import ceil import warnings @@ -8,8 +9,10 @@ import warnings from affine import Affine import numpy as np +import rasterio from rasterio._base import _transform -from rasterio._warp import _transform_geom, _reproject, Resampling +from rasterio._warp import (_transform_geom, _reproject, Resampling, + _calculate_default_transform) from rasterio.transform import guard_transform @@ -251,8 +254,7 @@ def calculate_default_transform( bottom, right, top, - resolution=None, - densify_pts=21): + resolution=None): """ Transforms bounds to destination coordinate system, calculates resolution if not provided, and returns destination transform and dimensions. @@ -260,13 +262,8 @@ def calculate_default_transform( Destination transform is anchored from the left, top coordinate. - Destination width and height are calculated from the number of pixels on - each dimension required to fit the destination bounds. - - If resolution is not provided, it is calculated using a weighted average - of the relative sizes of source width and height compared to the transformed - bounds (pixels are assumed to be square). - + Destination width and height (and resolution if not provided), are + calculated using GDAL's method for suggest warp output. Parameters ---------- @@ -283,37 +280,50 @@ def calculate_default_transform( Bounding coordinates in src_crs, from the bounds property of a raster. resolution: tuple (x resolution, y resolution) or float, optional Target resolution, in units of target coordinate reference system. - densify_pts: uint, optional - Number of points to add to each edge to account for nonlinear - edges produced by the transform process. Large numbers will produce - worse performance. Default: 21 (gdal default). Returns ------- tuple of destination affine transform, width, and height + + Note + ---- + Must be called within a raster.drivers() context + + Some behavior of this function is determined by the + CHECK_WITH_INVERT_PROJ environment variable + YES: constrain output raster to extents that can be inverted + avoids visual artifacts and coordinate discontinuties. + NO: reproject coordinates beyond valid bound limits """ + dst_affine, dst_width, dst_height = _calculate_default_transform( + src_crs, dst_crs, + width, height, + left, bottom, right, top) - xmin, ymin, xmax, ymax = transform_bounds( - src_crs, dst_crs, left, bottom, right, top, densify_pts) + # If resolution is specified, Keep upper-left anchored + # adjust the transform resolutions + # adjust the width/height by the ratio of estimated:specified res (ceil'd) + if resolution: - x_dif = xmax - xmin - y_dif = ymax - ymin - size = float(width + height) + # resolutions argument into tuple + try: + res = (float(resolution), float(resolution)) + except TypeError: + res = (resolution[0], resolution[0]) \ + if len(resolution) == 1 else resolution[0:2] - if resolution is None: - # TODO: compare to gdalwarp default (see - # _calculate_default_transform() in _warp.pyx. - avg_resolution = ( - (x_dif / float(width)) * (float(width) / size) + - (y_dif / float(height)) * (float(height) / size) - ) - resolution = (avg_resolution, avg_resolution) + # Assume yres is provided as positive, + # needs to be negative for north-up affine + xres = res[0] + yres = -res[1] - elif not isinstance(resolution, (tuple, list)): - resolution = (resolution, resolution) + xratio = dst_affine.a / xres + yratio = dst_affine.e / yres - dst_affine = Affine(resolution[0], 0, xmin, 0, -resolution[1], ymax) - dst_width = max(int(ceil(x_dif / resolution[0])), 1) - dst_height = max(int(ceil(y_dif / resolution[1])), 1) + dst_affine = Affine(xres, dst_affine.b, dst_affine.c, + dst_affine.d, yres, dst_affine.f) + + dst_width = ceil(dst_width * xratio) + dst_height = ceil(dst_height * yratio) return dst_affine, dst_width, dst_height diff --git a/tests/data/world.byte.tif b/tests/data/world.byte.tif new file mode 100644 index 00000000..17f420bd Binary files /dev/null and b/tests/data/world.byte.tif differ diff --git a/tests/test_rio_warp.py b/tests/test_rio_warp.py index 62ced43b..8aa779bb 100644 --- a/tests/test_rio_warp.py +++ b/tests/test_rio_warp.py @@ -114,11 +114,12 @@ def test_warp_reproject_dst_crs(runner, tmpdir): with rasterio.open(outputname) as output: assert output.count == src.count assert output.crs == {'init': 'epsg:4326'} - assert output.width == 824 - assert output.height == 686 + assert output.width == 835 + assert output.height == 696 assert numpy.allclose(output.bounds, - [-78.95864996545055, 23.564424693996177, - -76.57259451863895, 25.550873767433984]) + [-78.95864996545055, 23.564787976164418, + -76.5759177302349, 25.550873767433984]) + def test_warp_reproject_dst_crs_error(runner, tmpdir): srcname = 'tests/data/RGB.byte.tif' @@ -301,3 +302,18 @@ def test_warp_reproject_like(runner, tmpdir): assert numpy.allclose([0.001, 0.001], [output.affine.a, -output.affine.e]) assert output.width == 10 assert output.height == 10 + + +def test_warp_reproject_nolostdata(runner, tmpdir): + srcname = 'tests/data/world.byte.tif' + outputname = str(tmpdir.join('test.tif')) + result = runner.invoke(warp.warp, [srcname, outputname, + '--dst-crs', 'EPSG:3857']) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(outputname) as output: + arr = output.read() + # 50 column swath on the right edge should have some ones (gdalwarped has 7223) + assert arr[0, :, -50:].sum() > 7000 + assert output.crs == {'init': 'epsg:3857'} diff --git a/tests/test_warp.py b/tests/test_warp.py index f6fbfaec..3c1a4b80 100644 --- a/tests/test_warp.py +++ b/tests/test_warp.py @@ -131,9 +131,9 @@ def test_transform_bounds_densify_out_of_bounds(): def test_calculate_default_transform(): target_transform = Affine( - 0.0028956983577810586, 0.0, -78.95864996545055, - 0.0, -0.0028956983577810586, 25.550873767433984 - ) + 0.0028535715391804096, 0.0, -78.95864996545055, + 0.0, -0.0028535715391804096, 25.550873767433984) + with rasterio.drivers(): with rasterio.open('tests/data/RGB.byte.tif') as src: wgs84_crs = {'init': 'EPSG:4326'} @@ -141,8 +141,8 @@ def test_calculate_default_transform(): src.crs, wgs84_crs, src.width, src.height, *src.bounds) assert dst_transform.almost_equals(target_transform) - assert width == 824 - assert height == 686 + assert width == 835 + assert height == 696 def test_calculate_default_transform_single_resolution(): @@ -270,9 +270,9 @@ def test_reproject_nodata(): dst_nodata=nodata ) - assert (out == 1).sum() == 4461 + assert (out == 1).sum() == 6215 assert (out == nodata).sum() == (params.dst_width * - params.dst_height - 4461) + params.dst_height - 6215) def test_reproject_nodata_nan(): @@ -296,9 +296,9 @@ def test_reproject_nodata_nan(): dst_nodata=numpy.nan ) - assert (out == 1).sum() == 4461 + assert (out == 1).sum() == 6215 assert numpy.isnan(out).sum() == (params.dst_width * - params.dst_height - 4461) + params.dst_height - 6215) @@ -325,9 +325,9 @@ def test_reproject_dst_nodata_default(): dst_crs=params.dst_crs ) - assert (out == 1).sum() == 4461 + assert (out == 1).sum() == 6215 assert (out == 0).sum() == (params.dst_width * - params.dst_height - 4461) + params.dst_height - 6215) def test_reproject_invalid_dst_nodata(): diff --git a/tests/test_warp_transform.py b/tests/test_warp_transform.py index c1c3396f..8019df64 100644 --- a/tests/test_warp_transform.py +++ b/tests/test_warp_transform.py @@ -3,7 +3,7 @@ from rasterio._warp import _calculate_default_transform from rasterio.transform import Affine, from_bounds -def test_indentity(): +def test_identity(): """Get the same transform and dimensions back for same crs.""" # Tile: [53, 96, 8] # [-11740727.544603072, 4852834.0517692715, -11584184.510675032, 5009377.085697309] @@ -36,3 +36,28 @@ def test_gdal_transform_notnull(): right=-80, top=70) assert True + + +def test_gdal_transform_fail_dst_crs(): + with rasterio.drivers(): + dt, dw, dh = _calculate_default_transform( + {'init': 'EPSG:4326'}, + '+proj=foobar', + width=80, + height=80, + left=-120, + bottom=30, + right=-80, + top=70) + +def test_gdal_transform_fail_src_crs(): + with rasterio.drivers(): + dt, dw, dh = _calculate_default_transform( + '+proj=foobar', + {'init': 'EPSG:32610'}, + width=80, + height=80, + left=-120, + bottom=30, + right=-80, + top=70)