Added rio-rasterize

docs/cli.rst

@@ -143,6 +143,53 @@ datasets.
 
     $ rio merge rasterio/tests/data/R*.tif merged.tif
 
+rasterize
+---------
+
+New in 0.18.
+
+The rasterize command rasterizes GeoJSON features into a new or existing
+raster.
+
+.. code-block:: console
+
+    $ rio rasterize test.tif --res 0.0167 < input.geojson
+
+The resulting file will have an upper left coordinate determined by the bounds
+of the GeoJSON (in EPSG:4326, which is the default), with a
+pixel size of approximately 30 arc seconds.  Pixels whose center is within the
+polygon or that are selected by brezenhams line algorithm will be burned in
+with a default value of 1.
+
+It is possible to rasterize into an existing raster and use an alternative
+default value:
+
+.. code-block:: console
+
+    $ rio rasterize existing.tif --default_value 10 < input.geojson
+
+It is also possible to rasterize using a template raster, which will be used
+to determine the transform, dimensions, and coordinate reference system of the
+output raster:
+
+.. code-block:: console
+
+    $ rio rasterize test.tif --like tests/data/shade.tif < input.geojson
+
+GeoJSON features may be provided using stdin or specified directly as first
+argument, and dimensions may be provided in place of pixel resolution:
+
+.. code-block:: console
+
+    $ rio rasterize input.geojson test.tif --dimensions 1024 1024
+
+Other options are available, see:
+
+.. code-block:: console
+
+    $ rio rasterize --help
+
+
 shapes
 ------
 

rasterio/_features.pyx

@@ -291,6 +291,45 @@ def _rasterize(shapes, image, transform, all_touched):
             _gdal.CSLDestroy(options)
 
 
+def _explode(coords):
+    """Explode a GeoJSON geometry's coordinates object and yield
+    coordinate tuples. As long as the input is conforming, the type of
+    the geometry doesn't matter.  From Fiona 1.4.8"""
+    for e in coords:
+        if isinstance(e, (float, int)):
+            yield coords
+            break
+        else:
+            for f in _explode(e):
+                yield f
+
+
+def _bounds(geometry):
+    """Bounding box of a GeoJSON geometry.  
+    From Fiona 1.4.8 with updates here to handle feature collections.
+    TODO: add to Fiona.
+    """
+
+    try:  
+        if 'features' in geometry:
+            xmins = []
+            ymins = []
+            xmaxs = []
+            ymaxs = []
+            for feature in geometry['features']:
+                xmin, ymin, xmax, ymax = _bounds(feature['geometry'])
+                xmins.append(xmin)
+                ymins.append(ymin)
+                xmaxs.append(xmax)
+                ymaxs.append(ymax)
+            return min(xmins), min(ymins), max(xmaxs), max(ymaxs)
+        else:
+            xyz = tuple(zip(*list(_explode(geometry['coordinates']))))
+            return min(xyz[0]), min(xyz[1]), max(xyz[0]), max(xyz[1])
+    except (KeyError, TypeError):
+        return None
+
+
 # Mapping of OGR integer geometry types to GeoJSON type names.
 GEOMETRY_TYPES = {
     0: 'Unknown',

rasterio/features.py

@@ -8,7 +8,7 @@ import warnings
 import numpy as np
 
 import rasterio
-from rasterio._features import _shapes, _sieve, _rasterize
+from rasterio._features import _shapes, _sieve, _rasterize, _bounds
 from rasterio.transform import IDENTITY, guard_transform
 from rasterio.dtypes import get_minimum_int_dtype
 
@@ -321,3 +321,24 @@ def rasterize(
         _rasterize(valid_shapes, out, transform.to_gdal(), all_touched)
 
     return out
+
+
+def bounds(geometry):
+    """Returns a (minx, miny, maxx, maxy) bounding box.  From Fiona 1.4.8.
+    Modified to return bbox from geometry if available.
+
+    Parameters
+    ----------
+    geometry: GeoJSON-like feature, feature collection, or geometry.
+
+    Returns
+    -------
+    tuple
+        Bounding box: (minx, miny, maxx, maxy)
+    """
+
+    if 'bbox' in geometry:
+        return tuple(geometry['bbox'])
+
+    geom = geometry.get('geometry') or geometry
+    return _bounds(geom)

rasterio/rio/features.py

@@ -1,6 +1,8 @@
 import functools
 import json
 import logging
+from math import ceil, floor
+import os
 import sys
 import warnings
 
@@ -8,13 +10,15 @@ import click
 from cligj import (
     precision_opt, indent_opt, compact_opt, projection_geographic_opt,
     projection_projected_opt, sequence_opt, use_rs_opt,
-    geojson_type_feature_opt, geojson_type_bbox_opt)
+    geojson_type_feature_opt, geojson_type_bbox_opt, files_inout_arg,
+    format_opt)
 
 import rasterio
 from rasterio.transform import Affine
 from rasterio.rio.cli import cli, coords, write_features
 
 
+logger = logging.getLogger('rio')
 warnings.simplefilter('default')
 
 
@@ -142,3 +146,211 @@ def shapes(
     except Exception:
         logger.exception("Failed. Exception caught")
         sys.exit(1)
+
+
+# Rasterize command.
+@cli.command(short_help='Rasterize features.')
+@files_inout_arg
+@format_opt
+@click.option('--like', type=click.Path(exists=True),
+              help='Raster dataset to use as a template for obtaining affine '
+              'transform (bounds and resolution), crs, data type, and driver '
+              'used to create the output.  Only a single band will be output.')
+@click.option('--bounds', nargs=4, type=float, default=None,
+              help='Output bounds: left, bottom, right, top.')
+@click.option('--dimensions', nargs=2, type=int, default=None,
+              help='Output dataset width, height in number of pixels.')
+@click.option('--res', multiple=True, type=float, default=None,
+              help='Output dataset resolution in units of coordinate '
+              'reference system. Pixels assumed to be square if this option is '
+              'used once, otherwise use: '
+              '--res pixel_width --res pixel_height')
+@click.option('--src_crs', default='EPSG:4326',
+              help='Source coordinate reference system.  Limited to EPSG '
+              'codes for now.  Used as output coordinate system if output does '
+              'not exist or --like option is not used. Default: EPSG:4326')
+@click.option('--all_touched', is_flag=True, default=False)
+@click.option('--default_value', type=float, default=1, help='Default value '
+              'for rasterized pixels')
+@click.option('--fill', type=float, default=0, help='Fill value for all pixels '
+              'not overlapping features.  Will be evaluated as NoData pixels '
+              'for output.  Default: 0')
+@click.option('--property', type=str, default=None, help='Property in '
+              'GeoJSON features to use for rasterized values.  Any features '
+              'that lack this property will be given --default_value instead.')
+@click.pass_context
+def rasterize(
+        ctx,
+        files,
+        driver,
+        like,
+        bounds,
+        dimensions,
+        res,
+        src_crs,
+        all_touched,
+        default_value,
+        fill,
+        property):
+
+    """Rasterize GeoJSON into a new or existing raster.
+
+    If the output raster exists, rio-rasterize will rasterize feature values
+    into all bands of that raster.  The GeoJSON is assumed to be in the same
+    coordinate reference system as the output.
+
+    --default_value or property values when using --property must be using a
+    data type valid for the data type of that raster.
+
+
+    If a template raster is provided using the --like option, the affine
+    transform, coordinate reference system, and data type from that raster will
+    be used to create the output.
+
+    --default_value or property values when using --property must be using a
+    data type valid for the data type of that raster.
+
+    --driver, --bounds, --dimensions, --res, and --src_crs are ignored when
+    output exists or --like raster is provided
+
+
+    If the output does not exist and --like raster is not provided, the input
+    GeoJSON will be used to determine the bounds of the output unless
+    provided using --bounds.
+
+    --dimensions or --res are required in this case.
+
+    If --res is provided, the bottom and right coordinates of bounds are
+    ignored.
+
+
+    Note:
+    The GeoJSON is not projected to match the coordinate reference system
+    of the output or --like rasters at this time.  This functionality may be
+    added in the future.
+    """
+
+    import numpy as np
+    from rasterio.features import rasterize
+    from rasterio.features import bounds as calculate_bounds
+
+    verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+
+    files = list(files)
+    output = files.pop()
+    input = click.open_file(files.pop(0) if files else '-')
+
+    # If values are actually meant to be integers, we need to cast them
+    # as such or rasterize creates floating point outputs
+    if default_value == int(default_value):
+        default_value = int(default_value)
+    if fill == int(fill):
+        fill = int(fill)
+
+    with rasterio.drivers(CPL_DEBUG=verbosity > 2):
+
+        def feature_value(feature):
+            if property and 'properties' in feature:
+                return feature['properties'].get(property, default_value)
+            return default_value
+
+        geojson = json.loads(input.read())
+        if 'features' in geojson:
+            geometries = []
+            for f in geojson['features']:
+                geometries.append((f['geometry'], feature_value(f)))
+        elif 'geometry' in geojson:
+            geometries = ((geojson['geometry'], feature_value(geojson)), )
+        else:
+            raise click.BadParameter('Invalid GeoJSON', param=input,
+                                     param_hint='input')
+
+        if os.path.exists(output):
+            with rasterio.open(output, 'r+') as out:
+                meta = out.meta.copy()
+
+                result = rasterize(
+                    geometries,
+                    out_shape=(meta['height'], meta['width']),
+                    transform=meta.get('affine', meta['transform']),
+                    all_touched=all_touched,
+                    dtype=meta.get('dtype', None),
+                    default_value=default_value,
+                    fill = fill)
+
+                for bidx in range(1, meta['count'] + 1):
+                    data = out.read_band(bidx, masked=True)
+                    # Burn in any non-fill pixels, and update mask accordingly
+                    ne = result != fill
+                    data[ne] = result[ne]
+                    data.mask[ne] = False
+                    out.write_band(bidx, data)
+
+        else:
+            if like is not None:
+                template_ds = rasterio.open(like)
+                kwargs = template_ds.meta.copy()
+                kwargs['count'] = 1
+
+            else:
+                bounds = bounds or geojson.get('bbox', calculate_bounds(geojson))
+
+                if src_crs == 'EPSG:4326':
+                    if (bounds[0] < -180 or bounds[2] > 180 or
+                        bounds[1] < -80 or bounds[3] > 80):
+                        raise click.BadParameter(
+                            'Bounds are beyond the valid extent for EPSG:4326.',
+                            ctx, param=bounds, param_hint='--bounds')
+
+                if dimensions:
+                    width, height = dimensions
+                    res = (
+                        (bounds[2] - bounds[0]) / float(width),
+                        (bounds[3] - bounds[1]) / float(height)
+                    )
+
+                else:
+                    if not res:
+                        raise click.BadParameter(
+                            'pixel dimensions are required',
+                            ctx, param=res, param_hint='--res')
+                    elif len(res) == 1:
+                        res = (res[0], res[0])
+
+                    width = max(int(ceil((bounds[2] - bounds[0]) /
+                                float(res[0]))), 1)
+                    height = max(int(ceil((bounds[3] - bounds[1]) /
+                                 float(res[1]))), 1)
+
+                src_crs = src_crs.upper()
+                if not src_crs.count('EPSG:'):
+                    raise click.BadParameter(
+                        'invalid CRS.  Must be an EPSG code.',
+                        ctx, param=src_crs, param_hint='--src_crs')
+
+                kwargs = {
+                    'count': 1,
+                    'crs': src_crs,
+                    'width': width,
+                    'height': height,
+                    'transform': Affine(res[0], 0, bounds[0], 0, -res[1],
+                                        bounds[3]),
+                    'driver': driver
+                }
+
+            result = rasterize(
+                geometries,
+                out_shape=(kwargs['height'], kwargs['width']),
+                transform=kwargs.get('affine', kwargs['transform']),
+                all_touched=all_touched,
+                dtype=kwargs.get('dtype', None),
+                default_value=default_value,
+                fill = fill)
+
+            if 'dtype' not in kwargs:
+                kwargs['dtype'] = result.dtype
+
+            kwargs['nodata'] = fill
+
+            with rasterio.open(output, 'w', **kwargs) as out:
+                out.write_band(1, result)

rasterio/rio/main.py

@@ -2,7 +2,7 @@
 
 from rasterio.rio.cli import cli
 from rasterio.rio.bands import stack
-from rasterio.rio.features import shapes
+from rasterio.rio.features import shapes, rasterize
 from rasterio.rio.info import env, info
 from rasterio.rio.merge import merge
 from rasterio.rio.rio import bounds, insp, transform

tests/conftest.py

@@ -4,6 +4,8 @@ import os
 import sys
 
 import pytest
+from click.testing import CliRunner
+
 
 if sys.version_info > (3,):
     reduce = functools.reduce
@@ -11,6 +13,7 @@ if sys.version_info > (3,):
 test_files = [os.path.join(os.path.dirname(__file__), p) for p in [
     'data/RGB.byte.tif', 'data/float.tif', 'data/float_nan.tif', 'data/shade.tif']]
 
+
 def pytest_cmdline_main(config):
     # Bail if the test raster data is not present. Test data is not 
     # distributed with sdists since 0.12.
@@ -19,3 +22,8 @@ def pytest_cmdline_main(config):
     else:
         print("Test data not present. See download directions in tests/README.txt")
         sys.exit(1)
+
+
+@pytest.fixture(scope='function')
+def runner():
+    return CliRunner()

tests/test_features_bounds.py

@@ -0,0 +1,65 @@
+from rasterio.features import bounds
+
+
+# Tests copied from Fiona 1.4.1
+
+def test_bounds_point():
+    g = {'type': 'Point', 'coordinates': [10, 10]}
+    assert bounds(g) == (10, 10, 10, 10)
+
+
+def test_bounds_line():
+    g = {'type': 'LineString', 'coordinates': [[0, 0], [10, 10]]}
+    assert bounds(g) == (0, 0, 10, 10)
+
+
+def test_bounds_polygon():
+    g = {'type': 'Polygon', 'coordinates': [[[0, 0], [10, 10], [10, 0]]]}
+    assert bounds(g) == (0, 0, 10, 10)
+
+
+def test_bounds_z():
+    g = {'type': 'Point', 'coordinates': [10, 10, 10]}
+    assert bounds(g) == (10, 10, 10, 10)
+
+
+# TODO: add these to Fiona with update to bounds
+def test_bounds_existing_bbox():
+    """ Test with existing bbox in geojson, similar to that produced by
+    rasterio.  Values specifically modified here for testing, bboxes are not
+    valid as written.
+    """
+
+    fc = {
+        'bbox': [-107, 40, -105, 41],
+        'features': [{
+            'bbox': [-107, 40, -104, 42],
+            'geometry': {
+                'coordinates': [
+                    [[-107, 40], [-106, 40], [-106, 41], [-107, 41], [-107, 40]]
+                ],
+                'type': 'Polygon'
+            },
+            'type': 'Feature'
+        }],
+        'type': 'FeatureCollection'
+    }
+    assert bounds(fc['features'][0]) == (-107, 40, -104, 42)
+    assert bounds(fc) == (-107, 40, -105, 41)
+
+
+def test_feature_collection():
+    fc = {
+        'features': [{
+            'geometry': {
+                'coordinates': [
+                    [[-107, 40], [-106, 40], [-106, 41], [-107, 41], [-107, 40]]
+                ],
+                'type': 'Polygon'
+            },
+            'type': 'Feature'
+        }],
+        'type': 'FeatureCollection'
+    }
+    assert bounds(fc['features'][0]) == (-107, 40, -106, 41)
+    assert bounds(fc) == (-107, 40, -106, 41)

tests/test_rio_features.py

@@ -1,14 +1,81 @@
 import logging
+import os
 import re
 import sys
 
 import click
 from click.testing import CliRunner
 
+import rasterio
 from rasterio.rio import features
 
 
-logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+# logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+TEST_FEATURES = """{
+    "geometry": {
+        "coordinates": [
+            [
+                [-110, 40],
+                [-100, 40],
+                [-100, 45],
+                [-105, 45],
+                [-110, 40]
+            ]
+        ],
+        "type": "Polygon"
+    },
+    "properties": {
+        "val": 15
+    },
+    "type": "Feature"
+}"""
+
+
+# > rio shapes tests/data/shade.tif --mask --sampling 500 --projected --precision 0
+TEST_MERC_FEATURECOLLECTION = """{
+    "bbox": [-11858135.0, 4803914.0, -11848351.0, 4813698.0],
+    "features": [{
+        "bbox": [-11853357.504145855, 4808920.97837715,
+                 -11848580.189878704, 4813698.2926443005],
+        "geometry": {
+            "coordinates": [
+                [
+                    [-11853357.504145855, 4813698.2926443005],
+                    [-11853357.504145855, 4808920.97837715],
+                    [-11848580.189878704, 4808920.97837715],
+                    [-11848580.189878704, 4813698.2926443005],
+                    [-11853357.504145855, 4813698.2926443005]
+                ]
+            ],
+            "type": "Polygon"
+        },
+        "properties": {
+            "val": 2
+        },
+        "type": "Feature"
+    }, {
+        "bbox": [-11858134.818413004, 4804143.66411,
+                 -11853357.504145855, 4808920.97837715],
+        "geometry": {
+            "coordinates": [
+                [
+                    [-11858134.818413004, 4808920.97837715],
+                    [-11858134.818413004, 4804143.66411],
+                    [-11853357.504145855, 4804143.66411],
+                    [-11853357.504145855, 4808920.97837715],
+                    [-11858134.818413004, 4808920.97837715]
+                ]
+            ],
+            "type": "Polygon"
+        },
+        "properties": {
+            "val": 3
+        },
+        "type": "Feature"
+    }],
+    "type": "FeatureCollection"
+}"""
 
 
 def test_err():
@@ -18,24 +85,21 @@ def test_err():
     assert result.exit_code == 1
 
 
-def test_shapes():
-    runner = CliRunner()
+def test_shapes(runner):
     result = runner.invoke(features.shapes, ['tests/data/shade.tif'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
     assert result.output.count('"Feature"') == 232
 
 
-def test_shapes_sequence():
-    runner = CliRunner()
+def test_shapes_sequence(runner):
     result = runner.invoke(features.shapes, ['tests/data/shade.tif', '--sequence'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 0
     assert result.output.count('"Feature"') == 232
 
 
-def test_shapes_sequence_rs():
-    runner = CliRunner()
+def test_shapes_sequence_rs(runner):
     result = runner.invoke(
         features.shapes, [
             'tests/data/shade.tif',
@@ -47,24 +111,21 @@ def test_shapes_sequence_rs():
     assert result.output.count(u'\u001e') == 232
 
 
-def test_shapes_with_nodata():
-    runner = CliRunner()
+def test_shapes_with_nodata(runner):
     result = runner.invoke(features.shapes, ['tests/data/shade.tif', '--with-nodata'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
     assert result.output.count('"Feature"') == 288
 
 
-def test_shapes_indent():
-    runner = CliRunner()
+def test_shapes_indent(runner):
     result = runner.invoke(features.shapes, ['tests/data/shade.tif', '--indent', '2'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
     assert result.output.count('\n') == 70139
 
 
-def test_shapes_compact():
-    runner = CliRunner()
+def test_shapes_compact(runner):
     result = runner.invoke(features.shapes, ['tests/data/shade.tif', '--compact'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
@@ -72,8 +133,7 @@ def test_shapes_compact():
     assert result.output.count(': ') == 0
 
 
-def test_shapes_sampling():
-    runner = CliRunner()
+def test_shapes_sampling(runner):
     result = runner.invoke(
         features.shapes, ['tests/data/shade.tif', '--sampling', '10'])
     assert result.exit_code == 0
@@ -81,8 +141,7 @@ def test_shapes_sampling():
     assert result.output.count('"Feature"') == 124
 
 
-def test_shapes_precision():
-    runner = CliRunner()
+def test_shapes_precision(runner):
     result = runner.invoke(
         features.shapes, ['tests/data/shade.tif', '--precision', '1'])
     assert result.exit_code == 0
@@ -91,9 +150,171 @@ def test_shapes_precision():
     assert re.search(r'\d*\.\d{2,}', result.output) is None
 
 
-def test_shapes_mask():
-    runner = CliRunner()
+def test_shapes_mask(runner):
     result = runner.invoke(features.shapes, ['tests/data/RGB.byte.tif', '--mask'])
     assert result.exit_code == 0
     assert result.output.count('"FeatureCollection"') == 1
     assert result.output.count('"Feature"') == 9
+
+
+def test_rasterize_err(tmpdir, runner):
+    output = str(tmpdir.join('test.tif'))
+    # Test invalid stdin
+    result = runner.invoke(features.rasterize, [output], input='BOGUS')
+    assert result.exit_code == -1
+
+    # Test invalid GeoJSON
+    result = runner.invoke(features.rasterize, [output],
+                           input='{"foo": "bar"}')
+    assert result.exit_code == 2
+
+    # Test invalid res
+    result = runner.invoke(features.rasterize, [output], input=TEST_FEATURES)
+    assert result.exit_code == 2
+
+    # Test invalid CRS for bounds
+    result = runner.invoke(features.rasterize, [output, '--res', 1],
+                           input=TEST_MERC_FEATURECOLLECTION)
+    assert result.exit_code == 2
+
+    # Tets invalid CRS value
+    result = runner.invoke(features.rasterize, [output,
+                                                '--res', 1,
+                                                '--src_crs', 'BOGUS'],
+                           input=TEST_MERC_FEATURECOLLECTION)
+    assert result.exit_code == 2
+
+
+def test_rasterize(tmpdir, runner):
+    # Test dimensions
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, '--dimensions', 20, 10],
+                           input=TEST_FEATURES)
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        assert out.meta['width'] == 20
+        assert out.meta['height'] == 10
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 55
+        assert (data == 1).sum() == 145
+
+    # Test dimensions and bounds
+    output = str(tmpdir.join('test2.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output,
+                            '--dimensions', 40, 20,
+                            '--bounds', -120, 30, -90, 50
+                           ], input=TEST_FEATURES)
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        assert out.meta['width'] == 40
+        assert out.meta['height'] == 20
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 748
+        assert (data == 1).sum() == 52
+
+    # Test resolution
+    output = str(tmpdir.join('test3.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, '--res', 0.5], input=TEST_FEATURES)
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        assert out.meta['width'] == 20
+        assert out.meta['height'] == 10
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 55
+        assert (data == 1).sum() == 145
+
+
+def test_rasterize_existing_output(tmpdir, runner):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, '--res', 0.5], input=TEST_FEATURES)
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+
+    geojson = """{
+        "geometry": {
+            "coordinates": [
+                [
+                    [-102, 40],
+                    [-98, 40],
+                    [-98, 45],
+                    [-100, 45],
+                    [-102, 40]
+                ]
+            ],
+            "type": "Polygon"
+        },
+        "type": "Feature"
+    }"""
+
+    result = runner.invoke(features.rasterize, [output, '--default_value', 2],
+                           input=geojson)
+
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 55
+        assert (data == 1).sum() == 125
+        assert (data == 2).sum() == 20
+
+
+def test_rasterize_like(tmpdir, runner):
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, '--like', 'tests/data/shade.tif'],
+                           input=TEST_MERC_FEATURECOLLECTION)
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 548576
+        assert (data == 1).sum() == 500000
+
+    # Test invalid like raster
+    output = str(tmpdir.join('test2.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, '--like', 'foo.tif'], input=TEST_FEATURES)
+    assert result.exit_code == 2
+
+
+def test_rasterize_property_value(tmpdir, runner):
+    # Test feature collection property values
+    output = str(tmpdir.join('test.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, 
+                            '--res', 1000,
+                            '--property', 'val',
+                            '--src_crs', 'EPSG:3857'
+                           ],
+                           input=TEST_MERC_FEATURECOLLECTION)
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 50
+        assert (data == 2).sum() == 25
+        assert (data == 3).sum() == 25
+
+    # Test feature property values
+    output = str(tmpdir.join('test2.tif'))
+    result = runner.invoke(features.rasterize,
+                           [output, '--res', 0.5, '--property', 'val'],
+                           input=TEST_FEATURES)
+    assert result.exit_code == 0
+    assert os.path.exists(output)
+    with rasterio.open(output) as out:
+        assert out.count == 1
+        data = out.read_band(1, masked=False)
+        assert (data == 0).sum() == 55
+        assert (data == 15).sum() == 145