From ced7fa3bcd199f46625d4f312c8a19fcf4e73ff4 Mon Sep 17 00:00:00 2001 From: Brendan Ward Date: Wed, 1 Jul 2015 20:12:28 -0700 Subject: [PATCH] basic rio-warp in place --- rasterio/crs.py | 6 + rasterio/rio/features.py | 3 +- rasterio/rio/options.py | 5 + rasterio/rio/warp.py | 176 ++++++++++++++++++++++++++++++ rasterio/warp.py | 4 + setup.py | 1 + tests/test_rio_warp.py | 229 +++++++++++++++++++++++++++++++++++++++ tests/test_warp.py | 8 +- 8 files changed, 425 insertions(+), 7 deletions(-) create mode 100644 rasterio/rio/warp.py create mode 100644 tests/test_rio_warp.py diff --git a/rasterio/crs.py b/rasterio/crs.py index 7d3492ca..02767680 100644 --- a/rasterio/crs.py +++ b/rasterio/crs.py @@ -44,7 +44,13 @@ def from_string(prjs): Bare parameters like "+no_defs" are given a value of ``True``. All keys are checked against the ``all_proj_keys`` list. + + EPSG:XXXX is also allowed. """ + + if 'EPSG:' in prjs.upper(): + return from_epsg(prjs.split(':')[1]) + parts = [o.lstrip('+') for o in prjs.strip().split()] def parse(v): diff --git a/rasterio/rio/features.py b/rasterio/rio/features.py index 789c8181..47da9c32 100644 --- a/rasterio/rio/features.py +++ b/rasterio/rio/features.py @@ -363,8 +363,7 @@ def shapes( @format_opt @options.like_file_opt @options.bounds_opt -@click.option('--dimensions', nargs=2, type=int, default=None, - help='Output dataset width, height in number of pixels.') +@options.dimensions_opt @options.resolution_opt @click.option('--src-crs', '--src_crs', 'src_crs', default=None, help='Source coordinate reference system. Limited to EPSG ' diff --git a/rasterio/rio/options.py b/rasterio/rio/options.py index e53ce74f..76784861 100644 --- a/rasterio/rio/options.py +++ b/rasterio/rio/options.py @@ -105,6 +105,11 @@ bounds_opt = click.option( nargs=4, type=float, default=None, help='Output bounds: left, bottom, right, top.') +dimensions_opt = click.option( + '--dimensions', + nargs=2, type=int, default=None, + help='Output dataset width, height in number of pixels.') + dtype_opt = click.option( '-t', '--dtype', type=click.Choice([ diff --git a/rasterio/rio/warp.py b/rasterio/rio/warp.py new file mode 100644 index 00000000..12447add --- /dev/null +++ b/rasterio/rio/warp.py @@ -0,0 +1,176 @@ +import logging +from math import ceil +import click +from cligj import format_opt + +from . import options +import rasterio +from rasterio import crs +from rasterio.transform import Affine +from rasterio.warp import (reproject, RESAMPLING, calculate_default_transform, + transform_bounds) + + +logger = logging.getLogger('rio') + + +@click.command(short_help='Warp a raster dataset.') +@options.file_in_arg +@options.file_out_arg +@format_opt +@options.like_file_opt +@click.option('--dst-crs', default=None, + help='Target coordinate reference system. Default: EPSG:4326') +@options.dimensions_opt +@options.bounds_opt +#TODO: flag for bounds in target +@options.resolution_opt +@click.option('--resampling', type=click.Choice(['nearest', 'bilinear', 'cubic', + 'cubic_spline','lanczos', 'average', 'mode']), + default='nearest', help='Resampling method') +@click.option('--threads', type=int, default=1, + help='Number of processing threads.') +@click.pass_context +def warp( + ctx, + input, + output, + driver, + like, + dst_crs, + dimensions, + bounds, + res, + resampling, + threads): + """ + Warp a raster dataset. + + The output is always overwritten. + + If a template raster is provided using the --like option, the coordinate + reference system, affine transform, and dimensions of that raster will + be used for the output. In this case --dst-crs, --bounds, --res, and + --dimensions options are ignored. + + If --dimensions are provided, --res and --bounds are ignored. Resolution + is calculated based on the relationship between the source raster bounds + and dimensions, and may produce rectangular rather than square pixels. + + If --bounds are provided, --res is required if --dst-crs is provided + (defaults to source raster resolution otherwise). Bounds are in the source + coordinate reference system. + """ + + verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1 + resampling = getattr(RESAMPLING, resampling) # get integer code for method + + if not len(res): + # Click sets this as an empty tuple if not provided + res = None + else: + # 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.open(input) as src: + l, b, r, t = src.bounds + out_kwargs = src.meta.copy() + out_kwargs['driver'] = driver + + if like: + with rasterio.open(like) as template_ds: + dst_crs = template_ds.crs + dst_transform = template_ds.affine + dst_height = template_ds.height + dst_width = template_ds.width + + elif dst_crs: + if dimensions: + # Calculate resolution appropriate for dimensions in target + dst_width, dst_height = dimensions + xmin, ymin, xmax, ymax = transform_bounds(src.crs, dst_crs, + *src.bounds) + dst_transform = Affine( + (xmax - xmin) / float(dst_width), + 0, xmin, 0, + (ymin - ymax) / float(dst_height), + ymax + ) + + elif bounds: + if not res: + raise click.BadParameter('Required when using --bounds', + param='res', param_hint='res') + + xmin, ymin, xmax, ymax = transform_bounds(src.crs, dst_crs, + *bounds) + dst_transform = Affine(res[0], 0, xmin, 0, -res[1], ymax) + dst_width = max(int(ceil((xmax - xmin) / res[0])), 1) + dst_height = max(int(ceil((ymax - ymin) / res[1])), 1) + + else: + dst_crs = crs.from_string(dst_crs) + dst_transform, dst_width, dst_height = calculate_default_transform( + src.crs, dst_crs, src.width, src.height, *src.bounds, + resolution=res) + + elif dimensions: + # Same projection, different dimensions, calculate resolution + dst_crs = src.crs + dst_width, dst_height = dimensions + dst_transform = Affine( + (r - l) / float(dst_width), + 0, l, 0, + (b - t) / float(dst_height), + t + ) + + elif bounds: + # Same projection, different dimensions and possibly different + # resolution + if not res: + res = (src.affine.a, -src.affine.e) + + dst_crs = src.crs + xmin, ymin, xmax, ymax = bounds + dst_transform = Affine(res[0], 0, xmin, 0, -res[1], ymax) + dst_width = max(int(ceil((xmax - xmin) / res[0])), 1) + dst_height = max(int(ceil((ymax - ymin) / res[1])), 1) + + elif res: + # Same projection, different resolution + dst_crs = src.crs + dst_transform = Affine(res[0], 0, l, 0, -res[1], t) + dst_width = max(int(ceil((r - l) / res[0])), 1) + dst_height = max(int(ceil((t - b) / res[1])), 1) + + else: + dst_crs = src.crs + dst_transform = src.affine + dst_width = src.width + dst_height = src.height + + out_kwargs.update({ + 'crs': dst_crs, + 'transform': dst_transform, + 'affine': dst_transform, + 'width': dst_width, + 'height': dst_height + }) + + with rasterio.open(output, 'w', **out_kwargs) as dst: + for i in range(1, src.count + 1): + click.echo('Warping band {0}...'.format(i)) + + reproject( + source=rasterio.band(src, i), + destination=rasterio.band(dst, i), + src_transform=src.affine, + src_crs=src.crs, + # src_nodata=#TODO + dst_transform=out_kwargs['transform'], + dst_crs=out_kwargs['crs'], + # dst_nodata=#TODO + resampling=resampling, + num_threads=threads) \ No newline at end of file diff --git a/rasterio/warp.py b/rasterio/warp.py index 07e0ec78..34d384fb 100644 --- a/rasterio/warp.py +++ b/rasterio/warp.py @@ -262,6 +262,10 @@ def calculate_default_transform( Example: {'init': 'EPSG:4326'} dst_crs: dict Target coordinate reference system. + width: int + Source raster width. + height: int + Source raster height. left, bottom, right, top: float Bounding coordinates in src_crs, from the bounds property of a raster. resolution: tuple (x resolution, y resolution) or float, optional diff --git a/setup.py b/setup.py index 1b20b1f5..d5809052 100755 --- a/setup.py +++ b/setup.py @@ -244,6 +244,7 @@ setup_args = dict( sample=rasterio.rio.sample:sample shapes=rasterio.rio.features:shapes stack=rasterio.rio.bands:stack + warp=rasterio.rio.warp:warp transform=rasterio.rio.info:transform ''', include_package_data=True, diff --git a/tests/test_rio_warp.py b/tests/test_rio_warp.py new file mode 100644 index 00000000..86bb074f --- /dev/null +++ b/tests/test_rio_warp.py @@ -0,0 +1,229 @@ +import logging +import os +import re +import sys +import numpy + +import rasterio +from rasterio.rio import warp + + +logging.basicConfig(stream=sys.stderr, level=logging.DEBUG) + + +def test_warp_no_reproject(runner, tmpdir): + """ When called without parameters, output should be same as source """ + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + result = runner.invoke(warp.warp, [srcname, outputname]) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(srcname) as src: + with rasterio.open(outputname) as output: + assert output.count == src.count + assert output.crs == src.crs + assert output.nodata == src.nodata + assert numpy.allclose(output.bounds, src.bounds) + assert output.affine.almost_equals(src.affine) + assert numpy.allclose(output.read(1), src.read(1)) + + +def test_warp_no_reproject_dimensions(runner, tmpdir): + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + result = runner.invoke(warp.warp, [srcname, outputname, + '--dimensions', '100', '100']) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(srcname) as src: + with rasterio.open(outputname) as output: + assert output.crs == src.crs + assert output.width == 100 + assert output.height == 100 + assert numpy.allclose([97.839396, 97.839396], + [output.affine.a, -output.affine.e]) + + +def test_warp_no_reproject_res(runner, tmpdir): + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + result = runner.invoke(warp.warp, [srcname, outputname, + '--res', '30']) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(srcname) as src: + with rasterio.open(outputname) as output: + assert output.crs == src.crs + assert numpy.allclose([30, 30], [output.affine.a, -output.affine.e]) + assert output.width == 327 + assert output.height == 327 + + +def test_warp_no_reproject_bounds(runner, tmpdir): + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + out_bounds = [-11850000, 4810000, -11849000, 4812000] + result = runner.invoke(warp.warp,[srcname, outputname, + '--bounds'] + out_bounds) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(srcname) as src: + with rasterio.open(outputname) as output: + assert output.crs == src.crs + assert numpy.allclose(output.bounds, out_bounds) + assert numpy.allclose([src.affine.a, src.affine.e], + [output.affine.a, output.affine.e]) + assert output.width == 105 + assert output.height == 210 + + +def test_warp_no_reproject_bounds_res(runner, tmpdir): + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + out_bounds = [-11850000, 4810000, -11849000, 4812000] + result = runner.invoke(warp.warp,[srcname, outputname, + '--res', 30, '--bounds', ] + out_bounds) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(srcname) as src: + with rasterio.open(outputname) as output: + assert output.crs == src.crs + assert numpy.allclose(output.bounds, out_bounds) + assert numpy.allclose([30, 30], [output.affine.a, -output.affine.e]) + assert output.width == 34 + assert output.height == 67 + + +def test_warp_reproject_dst_crs(runner, tmpdir): + srcname = 'tests/data/RGB.byte.tif' + outputname = str(tmpdir.join('test.tif')) + result = runner.invoke(warp.warp, [srcname, outputname, + '--dst-crs', 'EPSG:4326']) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(srcname) as src: + 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 numpy.allclose(output.bounds, + [-78.95864996545055, 23.564424693996177, + -76.57259451863895, 25.550873767433984]) + + +def test_warp_reproject_dst_crs_proj4(runner, tmpdir): + proj4 = '+proj=longlat +ellps=WGS84 +datum=WGS84' + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + result = runner.invoke(warp.warp, [srcname, outputname, + '--dst-crs', proj4]) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(outputname) as output: + assert output.crs == {'init': 'epsg:4326'} # rasterio converts to EPSG + + +def test_warp_reproject_res(runner, tmpdir): + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + result = runner.invoke(warp.warp, [srcname, outputname, + '--dst-crs', 'EPSG:4326', + '--res', '0.01']) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(outputname) as output: + assert output.crs == {'init': 'epsg:4326'} + assert numpy.allclose([0.01, 0.01], [output.affine.a, -output.affine.e]) + assert output.width == 9 + assert output.height == 7 + + +def test_warp_reproject_dimensions(runner, tmpdir): + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + result = runner.invoke(warp.warp, [srcname, outputname, + '--dst-crs', 'EPSG:4326', + '--dimensions', '100', '100']) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(srcname) as src: + with rasterio.open(outputname) as output: + assert output.crs == {'init': 'epsg:4326'} + assert output.width == 100 + assert output.height == 100 + assert numpy.allclose([0.0008789062498762235, 0.0006771676143921468], + [output.affine.a, -output.affine.e]) + + +def test_warp_reproject_bounds_no_res(runner, tmpdir): + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + out_bounds = [-11850000, 4810000, -11849000, 4812000] + result = runner.invoke(warp.warp, [srcname, outputname, + '--dst-crs', 'EPSG:4326', + '--bounds', ] + out_bounds) + assert result.exit_code == 2 + + +def test_warp_reproject_bounds_res(runner, tmpdir): + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + out_bounds = [-11850000, 4810000, -11849000, 4812000] + result = runner.invoke(warp.warp, [srcname, outputname, + '--dst-crs', 'EPSG:4326', + '--res', 0.001, '--bounds', ] + + out_bounds) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(srcname) as src: + with rasterio.open(outputname) as output: + assert output.crs == {'init': 'epsg:4326'} + assert numpy.allclose(output.bounds[:], + [-106.45036, 39.6138, -106.44136, 39.6278]) + assert numpy.allclose([0.001, 0.001], + [output.affine.a, -output.affine.e]) + assert output.width == 9 + assert output.height == 14 + + +def test_warp_reproject_like(runner, tmpdir): + likename = str(tmpdir.join('like.tif')) + kwargs = { + "crs": {'init': 'epsg:4326'}, + "transform": (-106.523, 0.001, 0, 39.6395, 0, -0.001), + "count": 1, + "dtype": rasterio.uint8, + "driver": "GTiff", + "width": 10, + "height": 10, + "nodata": 0 + } + + with rasterio.drivers(): + with rasterio.open(likename, 'w', **kwargs) as dst: + data = numpy.zeros((10, 10), dtype=rasterio.uint8) + dst.write_band(1, data) + + srcname = 'tests/data/shade.tif' + outputname = str(tmpdir.join('test.tif')) + result = runner.invoke(warp.warp, [srcname, outputname, + '--like', likename]) + assert result.exit_code == 0 + assert os.path.exists(outputname) + + with rasterio.open(outputname) as output: + assert output.crs == {'init': 'epsg:4326'} + assert numpy.allclose([0.001, 0.001], [output.affine.a, -output.affine.e]) + assert output.width == 10 + assert output.height == 10 \ No newline at end of file diff --git a/tests/test_warp.py b/tests/test_warp.py index 9368f4e0..d64dfae7 100644 --- a/tests/test_warp.py +++ b/tests/test_warp.py @@ -136,10 +136,9 @@ def test_calculate_default_transform(): ) with rasterio.drivers(): with rasterio.open('tests/data/RGB.byte.tif') as src: - l, b, r, t = src.bounds wgs84_crs = {'init': 'EPSG:4326'} dst_transform, width, height = calculate_default_transform( - src.crs, wgs84_crs, src.width, src.height, l, b, r, t) + src.crs, wgs84_crs, src.width, src.height, *src.bounds) assert dst_transform.almost_equals(target_transform) assert width == 824 @@ -157,7 +156,7 @@ def test_calculate_default_transform_single_resolution(): ) dst_transform, width, height = calculate_default_transform( src.crs, {'init': 'EPSG:4326'}, src.width, src.height, - l, b, r, t, resolution=target_resolution + *src.bounds, resolution=target_resolution ) assert dst_transform.almost_equals(target_transform) @@ -168,7 +167,6 @@ def test_calculate_default_transform_single_resolution(): def test_calculate_default_transform_multiple_resolutions(): with rasterio.drivers(): with rasterio.open('tests/data/RGB.byte.tif') as src: - l, b, r, t = src.bounds target_resolution = (0.2, 0.1) target_transform = Affine( target_resolution[0], 0.0, -78.95864996545055, @@ -177,7 +175,7 @@ def test_calculate_default_transform_multiple_resolutions(): dst_transform, width, height = calculate_default_transform( src.crs, {'init': 'EPSG:4326'}, src.width, src.height, - l, b, r, t, resolution=target_resolution + *src.bounds, resolution=target_resolution ) assert dst_transform.almost_equals(target_transform)