From 1dd8c869820c3201d814cb15afc19e9bc39352db Mon Sep 17 00:00:00 2001 From: Sean Gillies Date: Wed, 18 Feb 2015 09:47:17 -0700 Subject: [PATCH] Add sieve and fix a bug in it. Sieving of images read directly from GDAL was broken, but is now fixed, with tests. --- docs/cli.rst | 7 ++++--- rasterio/_features.pyx | 3 ++- rasterio/features.py | 7 +++++-- rasterio/rio/calc.py | 12 ++++++++++++ requirements-dev.txt | 2 +- tests/test_features_sieve.py | 7 +++++++ tests/test_rio_calc.py | 15 +++++++++++++++ 7 files changed, 46 insertions(+), 7 deletions(-) diff --git a/docs/cli.rst b/docs/cli.rst index aa7cbc75..9764cb75 100644 --- a/docs/cli.rst +++ b/docs/cli.rst @@ -91,7 +91,7 @@ calc The calc command reads files as arrays, evaluates lisp-like expressions in their context, and writes the result as a new file. Members of the numpy -module and arithmetic and logical operators are surfaced as builtin functions +module and arithmetic and logical operators are available builtin functions and operators. It is intended for simple calculations; any calculations requiring multiple steps is better done in Python using the Rasterio and Numpy APIs. @@ -123,8 +123,9 @@ collects a sequence of 2-D arrays into a 3-D array for output. $ rio calc "(asarray (take a 1) (* (take a 2) (/ (mean (take a 1)) (mean (take a 2)))) (* (take a 3) (/ (mean (take a 1)) (mean (take a 3)))))" \ > --name a=tests/data/RGB.byte.tif --dtype ubyte /tmp/out.rgb.tif -The command above is also an example of a calculation that is a bit beyond the design of the calc command and something that could be done much more efficiently in -Python. +The command above is also an example of a calculation that is far beyond the +design of the calc command and something that could be done much more +efficiently in Python. info ---- diff --git a/rasterio/_features.pyx b/rasterio/_features.pyx index c7c63c6d..62d71cc3 100644 --- a/rasterio/_features.pyx +++ b/rasterio/_features.pyx @@ -173,11 +173,12 @@ def _sieve(image, size, output, mask, connectivity): in_band = in_mem_ds.band elif isinstance(image, tuple): rdr = image.ds - hband = rdr.band(image.bidx) + in_band = rdr.band(image.bidx) else: raise ValueError("Invalid source image") if isinstance(output, np.ndarray): + log.debug("Output array: %r", output) out_mem_ds = InMemoryRaster(output) out_band = out_mem_ds.band elif isinstance(output, tuple): diff --git a/rasterio/features.py b/rasterio/features.py index 3673c3a8..34115544 100644 --- a/rasterio/features.py +++ b/rasterio/features.py @@ -127,7 +127,7 @@ def sieve(image, size, out=None, output=None, mask=None, connectivity=4): if np.dtype(image.dtype).name not in valid_dtypes: valid_types_str = ', '.join(('rasterio.{0}'.format(t) for t in valid_dtypes)) - raise ValueError('image dtype must be one of: %' % valid_types_str) + raise ValueError('image dtype must be one of: %s' % valid_types_str) if size <= 0: raise ValueError('size must be greater than 0') @@ -155,7 +155,10 @@ def sieve(image, size, out=None, output=None, mask=None, connectivity=4): out = out if out is not None else output if out is None: - out = np.zeros_like(image) + if isinstance(image, tuple): + out = np.zeros(image.shape, image.dtype) + else: + out = np.zeros_like(image) else: if np.dtype(image.dtype).name != np.dtype(out.dtype).name: raise ValueError('output must match dtype of image') diff --git a/rasterio/rio/calc.py b/rasterio/rio/calc.py index f41f7c48..2c14bcf5 100644 --- a/rasterio/rio/calc.py +++ b/rasterio/rio/calc.py @@ -12,6 +12,8 @@ import snuggs import rasterio from rasterio.fill import fillnodata +from rasterio.features import sieve + from rasterio.rio.cli import cli @@ -25,6 +27,14 @@ def get_bands(inputs, d, i=None): return [rasterio.band(src, i) for i in src.indexes] +def read_array(ix, subix=None, dtype=None): + """Change the type of a read array""" + arr = snuggs._ctx.lookup(ix, subix) + if dtype: + arr = arr.astype(dtype) + return arr + + @cli.command(short_help="Raster data calculator.") @click.argument('command') @click.argument( @@ -115,9 +125,11 @@ def calc(ctx, command, files, name, dtype): src.read(), 'float64', copy=False) # Extend snuggs. + snuggs.func_map['read'] = read_array snuggs.func_map['band'] = lambda d, i: get_bands(inputs, d, i) snuggs.func_map['bands'] = lambda d: get_bands(inputs, d) snuggs.func_map['fillnodata'] = lambda *args: fillnodata(*args) + snuggs.func_map['sieve'] = lambda *args: sieve(*args) res = snuggs.eval(command, **ctxkwds) diff --git a/requirements-dev.txt b/requirements-dev.txt index 1e05e609..8fb6f6b2 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -5,7 +5,7 @@ cython==0.21.2 delocate enum34 numpy>=1.8.0 -snuggs>=1.0 +snuggs>=1.1 pytest setuptools>=0.9.8 wheel diff --git a/tests/test_features_sieve.py b/tests/test_features_sieve.py index baea5804..30dee1a5 100644 --- a/tests/test_features_sieve.py +++ b/tests/test_features_sieve.py @@ -139,3 +139,10 @@ def test_dtypes(): image = numpy.zeros((rows, cols), dtype=dtype) image[2:5, 2:5] = test_value sieved_image = ftrz.sieve(image, 2) + + +def test_sieve_shade(): + with rasterio.drivers(): + with rasterio.open('tests/data/shade.tif') as src: + sieved_image = ftrz.sieve(rasterio.band(src, 1), 42) + assert sieved_image.shape == (1024, 1024) diff --git a/tests/test_rio_calc.py b/tests/test_rio_calc.py index 5a0d5ac7..f6f7f5d6 100644 --- a/tests/test_rio_calc.py +++ b/tests/test_rio_calc.py @@ -154,3 +154,18 @@ def test_fillnodata_map(tmpdir): assert src.meta['dtype'] == 'uint8' data = src.read() assert round(data.mean(), 1) == 58.6 + + +def test_sieve_band(tmpdir): + outfile = str(tmpdir.join('out.tif')) + runner = CliRunner() + result = runner.invoke(calc, [ + '(sieve (band 1 1) 42)', + '--dtype', 'uint8', + 'tests/data/shade.tif', + outfile], + catch_exceptions=False) + assert result.exit_code == 0 + with rasterio.open(outfile) as src: + assert src.count == 1 + assert src.meta['dtype'] == 'uint8'