basic rio-warp in place

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):

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 '

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([

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)

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

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,

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

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)