[rasterio] 01/13: Imported Upstream version 0.27.0
Johan Van de Wauw
johanvdw-guest at moszumanska.debian.org
Tue Sep 29 18:25:15 UTC 2015
This is an automated email from the git hooks/post-receive script.
johanvdw-guest pushed a commit to branch master
in repository rasterio.
commit a423ade4b57ab4e66353d69215da323d871a0582
Author: Johan Van de Wauw <johan at vandewauw.be>
Date: Tue Sep 29 08:47:53 2015 +0200
Imported Upstream version 0.27.0
---
.travis.yml | 29 +++++--
CHANGES.txt | 16 ++++
docs/cli.rst | 26 +++++++
rasterio/__init__.py | 2 +-
rasterio/_base.pyx | 110 +++++++++++++++++++++++---
rasterio/_io.pyx | 16 ++--
rasterio/coords.py | 12 +++
rasterio/crs.py | 2 +-
rasterio/dtypes.py | 10 +++
rasterio/enums.py | 17 ++++
rasterio/rio/convert.py | 92 +++++++++++++++++++++-
rasterio/rio/features.py | 28 ++++---
rasterio/rio/info.py | 2 +-
rasterio/rio/merge.py | 175 ++++++++----------------------------------
rasterio/tool.py | 71 +++++++++++++++++
rasterio/tools/__init__.py | 0
rasterio/tools/merge.py | 164 +++++++++++++++++++++++++++++++++++++++
rasterio/transform.py | 9 +++
requirements-dev.txt | 2 +-
setup.py | 1 +
tests/data/rgb1.tif | Bin 0 -> 481148 bytes
tests/data/rgb2.tif | Bin 0 -> 471548 bytes
tests/data/rgb3.tif | Bin 0 -> 383844 bytes
tests/data/rgb4.tif | Bin 0 -> 376188 bytes
tests/data/rgb_deflate.tif | Bin 0 -> 912290 bytes
tests/test_crs.py | 9 +++
tests/test_image_structure.py | 74 ++++++++++++++++++
tests/test_indexing.py | 24 +++++-
tests/test_read_boundless.py | 15 ++++
tests/test_read_resample.py | 17 ++--
tests/test_rio_convert.py | 65 +++++++++++++++-
tests/test_rio_features.py | 6 +-
tests/test_rio_merge.py | 92 +++++++++++++++++-----
tests/test_tool.py | 23 +++++-
34 files changed, 885 insertions(+), 224 deletions(-)
diff --git a/.travis.yml b/.travis.yml
index bbf5cc5..955a27d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,16 +1,33 @@
language: python
+sudo: false
+cache:
+ directories:
+ - ~/.cache/pip
+env:
+ global:
+ - PIP_WHEEL_DIR=$HOME/.cache/pip/wheels
+ - PIP_FIND_LINKS=file://$HOME/.cache/pip/wheels
+addons:
+ apt:
+ packages:
+ - libgdal1h
+ - gdal-bin
+ - libgdal-dev
+ - libatlas-dev
+ - libatlas-base-dev
+ - gfortran
python:
- "2.7"
- "3.3"
- "3.4"
before_install:
- - sudo add-apt-repository -y ppa:ubuntugis/ppa
- - sudo apt-get update -qq
- - sudo apt-get install -y libgdal1h gdal-bin libgdal-dev
+ - pip install -U pip
+ - pip install wheel
install:
- - pip install --install-option="--no-cython-compile" cython==0.22
- - pip install -r requirements-dev.txt
- - pip install -e .
+ - "pip wheel -r requirements-dev.txt"
+ - "pip install -r requirements-dev.txt"
+ - "pip install coveralls"
+ - "pip install -e ."
script:
- py.test --cov rasterio --cov-report term-missing
after_success:
diff --git a/CHANGES.txt b/CHANGES.txt
index 96ee33a..01ea66e 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,6 +1,22 @@
Changes
=======
+0.27.0 (2015-09-25)
+-------------------
+- Ensure local uniqueness of the rio-shapes feature ids (#479).
+- Surface compression and interleaving as dataset properties and in rio-info
+ (#481). In the module, these are enums (`enums.Compression` and
+ `enums.Interleaving`); the values of the enums correspond to GDAL terms
+ (i.e, "DEFLATE") and the names are what surface in the CLI ("deflate").
+- Change get_window() and DatasetReader.window() to return a window guaranteed
+ to cover the input bounding box (#464, #485).
+- Bug fix for improperly computed transforms of output file in tools.merge and
+ rio-merge (#485).
+- More robust determination of dataset nodata values. In particular, the
+ absence of a nodata value is much more clear: dataset.nodata should never
+ return an out of range value when there is no nodata value, it should always
+ return `None` (#485).
+
0.26.0 (2015-08-11)
-------------------
- Add dependency on click-plugins, a new project that takes over the plugin
diff --git a/docs/cli.rst b/docs/cli.rst
index 56f560c..86e591e 100644
--- a/docs/cli.rst
+++ b/docs/cli.rst
@@ -19,6 +19,7 @@ Rasterio's command line interface is a program named "rio".
Commands:
bounds Write bounding boxes to stdout as GeoJSON.
calc Raster data calculator.
+ clip Clip a raster to given bounds.
convert Copy and convert raster dataset.
edit-info Edit dataset metadata.
env Print information about the rio environment.
@@ -157,6 +158,31 @@ efficiently in Python.
Please see `calc.rst <calc.rst>`__ for more details.
+clip
+----
+
+Added in 0.27
+
+The ``clip`` command clips a raster using bounds input directly or from a
+template raster.
+
+.. code-block:: console
+
+ $ rio clip input.tif output.tif --bounds xmin ymin xmax ymax
+ $ rio clip input.tif output.tif --like template.tif
+
+If using ``--bounds``, values must be in coordinate reference system of input.
+If using ``--like``, bounds will automatically be transformed to match the
+coordinate reference system of the input.
+
+It can also be combined to read bounds of a feature dataset using Fiona:
+
+.. code-block:: console
+
+ $ rio clip input.tif output.tif --bounds $(fio info features.shp --bounds)
+
+
+
convert
-------
diff --git a/rasterio/__init__.py b/rasterio/__init__.py
index 7cf55cf..11934fd 100644
--- a/rasterio/__init__.py
+++ b/rasterio/__init__.py
@@ -23,7 +23,7 @@ from rasterio import _err, coords, enums
__all__ = [
'band', 'open', 'drivers', 'copy', 'pad']
-__version__ = "0.26.0"
+__version__ = "0.27.0"
log = logging.getLogger('rasterio')
class NullHandler(logging.Handler):
diff --git a/rasterio/_base.pyx b/rasterio/_base.pyx
index caf5e85..422db77 100644
--- a/rasterio/_base.pyx
+++ b/rasterio/_base.pyx
@@ -15,7 +15,7 @@ from rasterio._err import cpl_errs
from rasterio import dtypes
from rasterio.coords import BoundingBox
from rasterio.transform import Affine
-from rasterio.enums import ColorInterp
+from rasterio.enums import ColorInterp, Compression, Interleaving
log = logging.getLogger('rasterio')
@@ -277,7 +277,7 @@ cdef class DatasetReader(object):
def get_nodatavals(self):
cdef void *hband = NULL
cdef double nodataval
- cdef int success
+ cdef int success = 0
if not self._nodatavals:
if self._hds == NULL:
@@ -286,11 +286,23 @@ cdef class DatasetReader(object):
hband = _gdal.GDALGetRasterBand(self._hds, i+1)
if hband == NULL:
raise ValueError("Null band")
+ dtype = dtypes.dtype_fwd[_gdal.GDALGetRasterDataType(hband)]
nodataval = _gdal.GDALGetRasterNoDataValue(hband, &success)
val = nodataval
- if not success:
+ # GDALGetRasterNoDataValue() has two ways of telling you that
+ # there's no nodata value. The success flag might come back
+ # 0 (FALSE). Even if it comes back 1 (TRUE), you still need
+ # to check that the return value is within the range of the
+ # data type. If so, the band has a nodata value. If not,
+ # there's no nodata value.
+ if (success == 0 or
+ val < dtypes.dtype_ranges[dtype][0] or
+ val > dtypes.dtype_ranges[dtype][1]):
val = None
+ log.debug("Nodata success: %d", success)
+ log.debug("Nodata value: %f", nodataval)
self._nodatavals.append(val)
+
return self._nodatavals
property nodatavals:
@@ -377,16 +389,15 @@ cdef class DatasetReader(object):
row += self.height
return c+a*col, f+e*row
- def index(self, x, y):
+ def index(self, x, y, op=math.floor):
"""Returns the (row, col) index of the pixel containing (x, y)."""
- a, b, c, d, e, f, _, _, _ = self.affine
- return int(math.floor((y-f)/e)), int(math.floor((x-c)/a))
+ return get_index(x, y, self.affine)
def window(self, left, bottom, right, top, boundless=False):
"""Returns the window corresponding to the world bounding box.
If boundless is False, window is limited to extent of this dataset."""
- EPS = 1.0e-8
- window = tuple(zip(self.index(left + EPS, top - EPS), self.index(right + EPS, bottom - EPS)))
+
+ window = get_window(left, bottom, right, top, self.affine)
if boundless:
return window
else:
@@ -421,6 +432,21 @@ cdef class DatasetReader(object):
self._read = True
return m
+ @property
+ def compression(self):
+ val = self.tags(ns='IMAGE_STRUCTURE').get('COMPRESSION')
+ if val:
+ return Compression(val)
+ else:
+ return None
+
+ @property
+ def interleaving(self):
+ val = self.tags(ns='IMAGE_STRUCTURE').get('INTERLEAVE')
+ if val:
+ return Interleaving(val)
+ else:
+ return None
property profile:
"""Basic metadata and creation options of this dataset.
@@ -435,9 +461,12 @@ cdef class DatasetReader(object):
blockxsize=self.block_shapes[0][1],
blockysize=self.block_shapes[0][0],
tiled=self.block_shapes[0][1] != self.width)
+ if self.compression:
+ m['compress'] = self.compression.name
+ if self.interleaving:
+ m['interleave'] = self.interleaving.name
return m
-
def lnglat(self):
w, s, e, n = self.bounds
cx = (w + e)/2.0
@@ -699,6 +728,69 @@ cpdef eval_window(object window, int height, int width):
return (r_start, r_stop), (c_start, c_stop)
+def get_index(x, y, affine, op=math.floor, precision=6):
+ """
+ Returns the (row, col) index of the pixel containing (x, y) given a
+ coordinate reference system.
+
+ Parameters
+ ----------
+ x : float
+ x value in coordinate reference system
+ y : float
+ y value in coordinate reference system
+ affine : tuple
+ Coefficients mapping pixel coordinates to coordinate reference system.
+ op : function
+ Function to convert fractional pixels to whole numbers (floor, ceiling,
+ round)
+ precision : int
+ Decimal places of precision in indexing, as in `round()`.
+
+ Returns
+ -------
+ row : int
+ row index
+ col : int
+ col index
+ """
+ # Use an epsilon, magnitude determined by the precision parameter
+ # and sign determined by the op function: positive for floor, negative
+ # for ceil.
+ eps = 10.0**-precision * (1.0 - 2.0*op(0.1))
+ row = int(op((y - eps - affine[5]) / affine[4]))
+ col = int(op((x + eps - affine[2]) / affine[0]))
+ return row, col
+
+
+def get_window(left, bottom, right, top, affine, precision=6):
+ """
+ Returns a window tuple given coordinate bounds and the coordinate reference
+ system.
+
+ Parameters
+ ----------
+ left : float
+ Left edge of window
+ bottom : float
+ Bottom edge of window
+ right : float
+ Right edge of window
+ top : float
+ top edge of window
+ affine : tuple
+ Coefficients mapping pixel coordinates to coordinate reference system.
+ precision : int
+ Decimal places of precision in indexing, as in `round()`.
+ """
+ window_start = get_index(
+ left, top, affine, op=math.floor, precision=precision)
+ window_stop = get_index(
+ right, bottom, affine, op=math.ceil, precision=precision)
+ window = tuple(zip(window_start, window_stop))
+ return window
+
+
def window_shape(window, height=-1, width=-1):
"""Returns shape of a window.
diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx
index c317678..1d42712 100644
--- a/rasterio/_io.pyx
+++ b/rasterio/_io.pyx
@@ -816,10 +816,10 @@ cdef class RasterReader(_base.DatasetReader):
else:
# Compute the overlap between the dataset and the boundless window.
overlap = ((
- max(min(window[0][0] or 0, self.height), 0),
- max(min(window[0][1] or self.height, self.height), 0)), (
- max(min(window[1][0] or 0, self.width), 0),
- max(min(window[1][1] or self.width, self.width), 0)))
+ max(min(window[0][0], self.height), 0),
+ max(min(window[0][1], self.height), 0)), (
+ max(min(window[1][0], self.width), 0),
+ max(min(window[1][1], self.width), 0)))
if overlap != ((0, 0), (0, 0)):
# Prepare a buffer.
@@ -978,10 +978,10 @@ cdef class RasterReader(_base.DatasetReader):
else:
# Compute the overlap between the dataset and the boundless window.
overlap = ((
- max(min(window[0][0] or 0, self.height), 0),
- max(min(window[0][1] or self.height, self.height), 0)), (
- max(min(window[1][0] or 0, self.width), 0),
- max(min(window[1][1] or self.width, self.width), 0)))
+ max(min(window[0][0], self.height), 0),
+ max(min(window[0][1], self.height), 0)), (
+ max(min(window[1][0], self.width), 0),
+ max(min(window[1][1], self.width), 0)))
if overlap != ((0, 0), (0, 0)):
# Prepare a buffer.
diff --git a/rasterio/coords.py b/rasterio/coords.py
index 1d89156..50eb59d 100644
--- a/rasterio/coords.py
+++ b/rasterio/coords.py
@@ -3,3 +3,15 @@ from collections import namedtuple
BoundingBox = namedtuple('BoundingBox', ('left', 'bottom', 'right', 'top'))
+
+def disjoint_bounds(bounds1, bounds2):
+ """Returns True if bounds do not overlap
+
+ Parameters
+ ----------
+ bounds1: rasterio bounds tuple (xmin, ymin, xmax, ymax)
+ bounds2: rasterio bounds tuple
+ """
+
+ return (bounds1[0] > bounds2[2] or bounds1[2] < bounds2[0] or
+ bounds1[1] > bounds2[3] or bounds1[3] < bounds2[1])
\ No newline at end of file
diff --git a/rasterio/crs.py b/rasterio/crs.py
index 775a54f..2ee4d45 100644
--- a/rasterio/crs.py
+++ b/rasterio/crs.py
@@ -58,7 +58,7 @@ def from_string(prjs):
except ValueError:
raise ValueError('crs appears to be JSON but is not valid')
- if 'EPSG:' in prjs.upper():
+ if prjs.strip().upper().startswith('EPSG:'):
return from_epsg(prjs.split(':')[1])
parts = [o.lstrip('+') for o in prjs.strip().split()]
diff --git a/rasterio/dtypes.py b/rasterio/dtypes.py
index 449ec5f..90097d8 100644
--- a/rasterio/dtypes.py
+++ b/rasterio/dtypes.py
@@ -58,6 +58,16 @@ typename_fwd = {
typename_rev = dict((v, k) for k, v in typename_fwd.items())
+dtype_ranges = {
+ 'uint8': (0, 255),
+ 'uint16': (0, 65535),
+ 'int16': (-32768, 32767),
+ 'uint32': (0, 4294967295),
+ 'int32': (-2147483648, 2147483647),
+ 'float32': (-3.4028235e+38, 3.4028235e+38),
+ 'float64': (-1.7976931348623157e+308, 1.7976931348623157e+308)}
+
+
def _gdal_typename(dt):
try:
return typename_fwd[dtype_rev[dt]]
diff --git a/rasterio/enums.py b/rasterio/enums.py
index d33f71d..3e347a5 100644
--- a/rasterio/enums.py
+++ b/rasterio/enums.py
@@ -28,3 +28,20 @@ class Resampling(Enum):
mode='MODE'
average_magphase='AVERAGE_MAGPHASE'
none='NONE'
+
+
+class Compression(Enum):
+ jpeg='JPEG'
+ lzw='LZW'
+ packbits='PACKBITS'
+ deflate='DEFLATE'
+ ccittrle='CCITTRLE'
+ ccittfax3='CCITTFAX3'
+ ccittfax4='CCITTFAX4'
+ lzma='LZMA'
+ none='NONE'
+
+
+class Interleaving(Enum):
+ pixel='PIXEL'
+ band='BAND'
diff --git a/rasterio/rio/convert.py b/rasterio/rio/convert.py
index 3443b51..11781f7 100644
--- a/rasterio/rio/convert.py
+++ b/rasterio/rio/convert.py
@@ -4,17 +4,107 @@ import logging
import warnings
import click
-from cligj import files_inout_arg, format_opt
+from cligj import format_opt
import numpy as np
from .helpers import resolve_inout
from . import options
import rasterio
+from rasterio.coords import disjoint_bounds
warnings.simplefilter('default')
+# Clip command
+ at click.command(short_help='Clip a raster to given bounds.')
+ at click.argument(
+ 'files',
+ nargs=-1,
+ type=click.Path(resolve_path=True),
+ required=True,
+ metavar="INPUT OUTPUT")
+ at options.output_opt
+ at options.bounds_opt
+ at click.option(
+ '--like',
+ type=click.Path(exists=True),
+ help='Raster dataset to use as a template for bounds')
+ at format_opt
+ at options.creation_options
+ at click.pass_context
+def clip(
+ ctx,
+ files,
+ output,
+ bounds,
+ like,
+ driver,
+ creation_options):
+ """Clips a raster using bounds input directly or from a template raster.
+
+ \b
+ $ rio clip input.tif output.tif --bounds xmin ymin xmax ymax
+ $ rio clip input.tif output.tif --like template.tif
+
+ If using --bounds, values must be in coordinate reference system of input.
+ If using --like, bounds will automatically be transformed to match the
+ coordinate reference system of the input.
+
+ It can also be combined to read bounds of a feature dataset using Fiona:
+
+ \b
+ $ rio clip input.tif output.tif --bounds $(fio info features.shp --bounds)
+
+ """
+
+ from rasterio.warp import transform_bounds
+
+ verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+
+ with rasterio.drivers(CPL_DEBUG=verbosity > 2):
+
+ output, files = resolve_inout(files=files, output=output)
+ input = files[0]
+
+ with rasterio.open(input) as src:
+ if bounds:
+ if disjoint_bounds(bounds, src.bounds):
+ raise click.BadParameter('must overlap the extent of '
+ 'the input raster',
+ param='--bounds',
+ param_hint='--bounds')
+ elif like:
+ with rasterio.open(like) as template_ds:
+ bounds = template_ds.bounds
+ if template_ds.crs != src.crs:
+ bounds = transform_bounds(template_ds.crs, src.crs,
+ *bounds)
+
+ if disjoint_bounds(bounds, src.bounds):
+ raise click.BadParameter('must overlap the extent of '
+ 'the input raster',
+ param='--like',
+ param_hint='--like')
+
+ else:
+ raise click.UsageError('--bounds or --like required')
+
+ window = src.window(*bounds)
+
+ out_kwargs = src.meta.copy()
+ out_kwargs.update({
+ 'driver': driver,
+ 'height': window[0][1] - window[0][0],
+ 'width': window[1][1] - window[1][0],
+ 'transform': src.window_transform(window)
+ })
+ out_kwargs.update(**creation_options)
+
+ with rasterio.open(output, 'w', **out_kwargs) as out:
+ out.write(src.read(window=window))
+
+
@click.command(short_help="Copy and convert raster dataset.")
@click.argument(
'files',
diff --git a/rasterio/rio/features.py b/rasterio/rio/features.py
index 31d84c2..0345bbb 100644
--- a/rasterio/rio/features.py
+++ b/rasterio/rio/features.py
@@ -3,6 +3,7 @@ import logging
from math import ceil
import os
import shutil
+import re
import click
import cligj
@@ -16,7 +17,7 @@ from .helpers import coords, resolve_inout, write_features, to_lower
from . import options
import rasterio
from rasterio.transform import Affine
-
+from rasterio.coords import disjoint_bounds
logger = logging.getLogger('rio')
@@ -121,10 +122,10 @@ def mask(
bounds = geojson.get('bbox', calculate_bounds(geojson))
with rasterio.open(input) as src:
- disjoint_bounds = _disjoint_bounds(bounds, src.bounds)
+ has_disjoint_bounds = disjoint_bounds(bounds, src.bounds)
if crop:
- if disjoint_bounds:
+ if has_disjoint_bounds:
raise click.BadParameter('not allowed for GeoJSON outside '
'the extent of the input raster',
param=crop, param_hint='--crop')
@@ -134,7 +135,7 @@ def mask(
(r1, r2), (c1, c2) = window
mask_shape = (r2 - r1, c2 - c1)
else:
- if disjoint_bounds:
+ if has_disjoint_bounds:
click.echo('GeoJSON outside bounds of existing output '
'raster. Are they in different coordinate '
'reference systems?',
@@ -329,8 +330,11 @@ def shapes(
if not with_nodata:
kwargs['mask'] = msk
+ src_basename = os.path.basename(src.name)
+
# Yield GeoJSON features.
- for g, i in rasterio.features.shapes(img, **kwargs):
+ for i, (g, val) in enumerate(
+ rasterio.features.shapes(img, **kwargs)):
if projection == 'geographic':
g = rasterio.warp.transform_geom(
src.crs, 'EPSG:4326', g,
@@ -338,9 +342,9 @@ def shapes(
xs, ys = zip(*coords(g))
yield {
'type': 'Feature',
- 'id': str(i),
+ 'id': "{0}:{1}".format(src_basename, i),
'properties': {
- 'val': i, 'filename': os.path.basename(src.name)
+ 'val': val, 'filename': src_basename
},
'bbox': [min(xs), min(ys), max(xs), max(ys)],
'geometry': g
@@ -487,7 +491,7 @@ def rasterize(
'existing output raster',
param='input', param_hint='input')
- if _disjoint_bounds(geojson_bounds, out.bounds):
+ if disjoint_bounds(geojson_bounds, out.bounds):
click.echo("GeoJSON outside bounds of existing output "
"raster. Are they in different coordinate "
"reference systems?",
@@ -521,7 +525,7 @@ def rasterize(
'--like raster',
param='input', param_hint='input')
- if _disjoint_bounds(geojson_bounds, template_ds.bounds):
+ if disjoint_bounds(geojson_bounds, template_ds.bounds):
click.echo("GeoJSON outside bounds of --like raster. "
"Are they in different coordinate reference "
"systems?",
@@ -702,9 +706,3 @@ def bounds(ctx, input, precision, indent, compact, projection, dst_crs,
except Exception:
logger.exception("Exception caught during processing")
raise click.Abort()
-
-
-def _disjoint_bounds(bounds1, bounds2):
- return (bounds1[0] > bounds2[2] or bounds1[2] < bounds2[0] or
- bounds1[1] > bounds2[3] or bounds1[3] < bounds2[1])
-
diff --git a/rasterio/rio/info.py b/rasterio/rio/info.py
index d86d014..71f5a37 100644
--- a/rasterio/rio/info.py
+++ b/rasterio/rio/info.py
@@ -278,7 +278,7 @@ def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx,
try:
with rasterio.drivers(CPL_DEBUG=(verbosity > 2)):
with rasterio.open(input, mode) as src:
- info = src.meta
+ info = src.profile
info['transform'] = info['affine'][:6]
del info['affine']
info['shape'] = info['height'], info['width']
diff --git a/rasterio/rio/merge.py b/rasterio/rio/merge.py
index 0d15e97..d816ac7 100644
--- a/rasterio/rio/merge.py
+++ b/rasterio/rio/merge.py
@@ -1,6 +1,5 @@
# Merge command.
-
import logging
import math
import os.path
@@ -23,9 +22,17 @@ from rasterio.transform import Affine
@options.resolution_opt
@click.option('--nodata', type=float, default=None,
help="Override nodata values defined in input datasets")
+ at click.option('--force-overwrite', '-f', 'force_overwrite', is_flag=True,
+ type=bool, default=False,
+ help="Do not prompt for confirmation before overwriting output "
+ "file")
+ at click.option('--precision', type=int, default=7,
+ help="Number of decimal places of precision in alignment of "
+ "pixels")
@options.creation_options
@click.pass_context
-def merge(ctx, files, output, driver, bounds, res, nodata, creation_options):
+def merge(ctx, files, output, driver, bounds, res, nodata, force_overwrite,
+ precision, creation_options):
"""Copy valid pixels from input files to an output file.
All files must have the same number of bands, data type, and
@@ -43,148 +50,30 @@ def merge(ctx, files, output, driver, bounds, res, nodata, creation_options):
--res 0.1 0.1 => --res 0.1 (square)
--res 0.1 0.2 => --res 0.1 --res 0.2 (rectangular)
"""
- import numpy as np
+
+ from rasterio.tools.merge import merge as merge_tool
verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
logger = logging.getLogger('rio')
- if len(res) == 1:
- res = (res[0], res[0])
-
- try:
- with rasterio.drivers(CPL_DEBUG=verbosity>2):
- output, files = resolve_inout(files=files, output=output)
-
- with rasterio.open(files[0]) as first:
- first_res = first.res
- kwargs = first.meta
- kwargs.update(**creation_options)
- kwargs.pop('affine')
- nodataval = first.nodatavals[0]
- dtype = first.dtypes[0]
-
- if os.path.exists(output):
- # TODO: prompt user to update existing file (-i option) like:
- # overwrite b.tif? (y/n [n]) n
- # not overwritten
- dst = rasterio.open(output, 'r+')
- nodataval = dst.nodatavals[0]
- dtype = dst.dtypes[0]
- dest = np.zeros((dst.count,) + dst.shape, dtype=dtype)
- else:
- # Create new output file.
- # Extent from option or extent of all inputs.
- if not bounds:
- # scan input files.
- xs = []
- ys = []
- for f in files:
- with rasterio.open(f) as src:
- left, bottom, right, top = src.bounds
- xs.extend([left, right])
- ys.extend([bottom, top])
- bounds = min(xs), min(ys), max(xs), max(ys)
- output_transform = Affine.translation(bounds[0], bounds[3])
-
- # Resolution/pixel size.
- if not res:
- res = first_res
- output_transform *= Affine.scale(res[0], -res[1])
-
- # Dataset shape.
- output_width = int(math.ceil((bounds[2]-bounds[0])/res[0]))
- output_height = int(math.ceil((bounds[3]-bounds[1])/res[1]))
-
- kwargs['driver'] = driver
- kwargs['transform'] = output_transform
- kwargs['width'] = output_width
- kwargs['height'] = output_height
-
- logger.debug("Kwargs: %r", kwargs)
- logger.debug("bounds: %r", bounds)
- logger.debug("Res: %r", res)
-
- dst = rasterio.open(output, 'w', **kwargs)
- dest = np.zeros((first.count, output_height, output_width),
- dtype=dtype)
-
- logger.debug("In merge, dest shape: %r", dest.shape)
-
- if nodata is not None:
- nodataval = nodata
-
- if nodataval is not None:
- # Only fill if the nodataval is within dtype's range.
- inrange = False
- if np.dtype(dtype).kind in ('i', 'u'):
- info = np.iinfo(dtype)
- inrange = (info.min <= nodataval <= info.max)
- elif np.dtype(dtype).kind == 'f':
- info = np.finfo(dtype)
- inrange = (info.min <= nodataval <= info.max)
- if inrange:
- dest.fill(nodataval)
- else:
- warnings.warn(
- "Input file's nodata value, %s, is beyond the valid "
- "range of its data type, %s. Consider overriding it "
- "using the --nodata option for better results." % (
- nodataval, dtype))
- else:
- nodataval = 0
-
- dst_w, dst_s, dst_e, dst_n = dst.bounds
-
- for fname in reversed(files):
- with rasterio.open(fname) as src:
- # Real World (tm) use of boundless reads.
- # This approach uses the maximum amount of memory to solve
- # the problem. Making it more efficient is a TODO.
-
- # 1. Compute spatial intersection of destination
- # and source.
- src_w, src_s, src_e, src_n = src.bounds
-
- int_w = src_w if src_w > dst_w else dst_w
- int_s = src_s if src_s > dst_s else dst_s
- int_e = src_e if src_e < dst_e else dst_e
- int_n = src_n if src_n < dst_n else dst_n
-
- # 2. Compute the source window.
- src_window = src.window(int_w, int_s, int_e, int_n)
-
- # 3. Compute the destination window.
- dst_window = dst.window(int_w, int_s, int_e, int_n)
-
- # 4. Initialize temp array.
- temp = np.zeros(
- (first.count,) + tuple(b - a for a, b in dst_window),
- dtype=dtype)
-
- temp = src.read(
- out=temp,
- window=src_window,
- boundless=False,
- masked=True)
-
- # 5. Copy elements of temp into dest.
- roff, coff = dst.index(int_w, int_n)
- h, w = temp.shape[-2:]
-
- region = dest[:,roff:roff+h,coff:coff+w]
- np.copyto(region, temp,
- where=np.logical_and(
- region==nodataval, temp.mask==False))
-
- if dst.mode == 'r+':
- temp = dst.read(masked=True)
- np.copyto(dest, temp,
- where=np.logical_and(
- dest==nodataval, temp.mask==False))
-
- dst.write(dest)
- dst.close()
-
- except Exception:
- logger.exception("Exception caught during processing")
- raise click.Abort()
+ output, files = resolve_inout(files=files, output=output)
+
+ if os.path.exists(output) and not force_overwrite:
+ raise click.ClickException(
+ "Output exists and won't be overwritten without the "
+ "`-f` option")
+
+ sources = [rasterio.open(f) for f in files]
+ dest, output_transform = merge_tool(sources, bounds=bounds, res=res,
+ nodata=nodata, precision=precision)
+
+ profile = sources[0].profile
+ profile.pop('affine')
+ profile['transform'] = output_transform
+ profile['height'] = dest.shape[1]
+ profile['width'] = dest.shape[2]
+ profile['driver'] = driver
+ profile.update(**creation_options)
+
+ with rasterio.open(output, 'w', **profile) as dst:
+ dst.write(dest)
diff --git a/rasterio/tool.py b/rasterio/tool.py
index cb6364f..7804871 100644
--- a/rasterio/tool.py
+++ b/rasterio/tool.py
@@ -11,6 +11,7 @@ except ImportError:
import numpy
import rasterio
+from rasterio.five import zip_longest
logger = logging.getLogger('rasterio')
@@ -48,6 +49,76 @@ def stats(source):
return Stats(numpy.min(arr), numpy.max(arr), numpy.mean(arr))
+def show_hist(source, bins=10, masked=True, title='Histogram'):
+
+ """
+ Easily display a histogram with matplotlib.
+
+ Parameters
+ ----------
+ bins : int, optional
+ Compute histogram across N bins.
+ data : np.array or rasterio.Band or tuple(dataset, bidx)
+ Input data to display. The first three arrays in multi-dimensional
+ arrays are plotted as red, green, and blue.
+ masked : bool, optional
+ When working with a `rasterio.Band()` object, specifies if the data
+ should be masked on read.
+ title : str, optional
+ Title for the figure.
+ """
+
+ if plt is None:
+ raise ImportError("Could not import matplotlib")
+
+ if isinstance(source, (tuple, rasterio.Band)):
+ arr = source[0].read(source[1], masked=masked)
+ else:
+ arr = source
+
+ # The histogram is computed individually for each 'band' in the array
+ # so we need the overall min/max to constrain the plot
+ rng = arr.min(), arr.max()
+
+ if len(arr.shape) is 2:
+ arr = [arr]
+ colors = ['gold']
+ else:
+ colors = ('red', 'green', 'blue', 'violet', 'gold', 'saddlebrown')
+
+ # If a rasterio.Band() is given make sure the proper index is displayed
+ # in the legend.
+ if isinstance(source, (tuple, rasterio.Band)):
+ labels = [str(source[1])]
+ else:
+ labels = (str(i + 1) for i in range(len(arr)))
+
+ # This loop should add a single plot each band in the input array,
+ # regardless of if the number of bands exceeds the number of colors.
+ # The colors slicing ensures that the number of iterations always
+ # matches the number of bands.
+ # The goal is to provide a curated set of colors for working with
+ # smaller datasets and let matplotlib define additional colors when
+ # working with larger datasets.
+ for bnd, color, label in zip_longest(arr, colors[:len(arr)], labels):
+
+ plt.hist(
+ bnd.flatten(),
+ bins=bins,
+ alpha=0.5,
+ color=color,
+ label=label,
+ range=rng
+ )
+
+ plt.legend(loc="upper right")
+ plt.title(title, fontweight='bold')
+ plt.grid(True)
+ plt.xlabel('DN')
+ plt.ylabel('Frequency')
+ plt.show()
+
+
def main(banner, dataset, alt_interpreter=None):
""" Main entry point for use with python interpreter """
local = dict(funcs, src=dataset, np=numpy, rio=rasterio, plt=plt)
diff --git a/rasterio/tools/__init__.py b/rasterio/tools/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/rasterio/tools/merge.py b/rasterio/tools/merge.py
new file mode 100644
index 0000000..0fea083
--- /dev/null
+++ b/rasterio/tools/merge.py
@@ -0,0 +1,164 @@
+import logging
+import math
+import warnings
+
+import numpy as np
+
+import rasterio
+from rasterio._base import get_index, get_window
+from rasterio.transform import Affine
+
+
+logger = logging.getLogger('rasterio')
+
+
+def merge(sources, bounds=None, res=None, nodata=None, precision=7):
+ """Copy valid pixels from input files to an output file.
+
+ All files must have the same number of bands, data type, and
+ coordinate reference system.
+
+ Input files are merged in their listed order using the reverse
+ painter's algorithm. If the output file exists, its values will be
+ overwritten by input values.
+
+ Geospatial bounds and resolution of a new output file in the
+ units of the input file coordinate reference system may be provided
+ and are otherwise taken from the first input file.
+
+ Parameters
+ ----------
+ sources: list of source datasets
+ Open rasterio RasterReader objects to be merged.
+ bounds: tuple, optional
+ Bounds of the output image (left, bottom, right, top).
+ If not set, bounds are determined from bounds of input rasters.
+ res: tuple, optional
+ Output resolution in units of coordinate reference system. If not set,
+ the resolution of the first raster is used. If a single value is passed,
+ output pixels will be square.
+ nodata: float, optional
+ nodata value to use in output file. If not set, uses the nodata value
+ in the first input raster.
+
+ Returns
+ -------
+ dest: numpy ndarray
+ Contents of all input rasters in single array.
+ out_transform: affine object
+ Information for mapping pixel coordinates in `dest` to another
+ coordinate system
+ """
+
+ first = sources[0]
+ first_res = first.res
+ nodataval = first.nodatavals[0]
+ dtype = first.dtypes[0]
+
+ # Extent from option or extent of all inputs.
+ if bounds:
+ dst_w, dst_s, dst_e, dst_n = bounds
+ else:
+ # scan input files.
+ xs = []
+ ys = []
+ for src in sources:
+ left, bottom, right, top = src.bounds
+ xs.extend([left, right])
+ ys.extend([bottom, top])
+ dst_w, dst_s, dst_e, dst_n = min(xs), min(ys), max(xs), max(ys)
+
+ logger.debug("Output bounds: %r", (dst_w, dst_s, dst_e, dst_n))
+ output_transform = Affine.translation(dst_w, dst_n)
+ logger.debug("Output transform, before scaling: %r", output_transform)
+
+ # Resolution/pixel size.
+ if not res:
+ res = first_res
+ elif not np.iterable(res):
+ res = (res, res)
+ elif len(res) == 1:
+ res = (res[0], res[0])
+ output_transform *= Affine.scale(res[0], -res[1])
+ logger.debug("Output transform, after scaling: %r", output_transform)
+
+ # Compute output array shape. We guarantee it will cover the output
+ # bounds completely.
+ output_width = int(math.ceil((dst_e - dst_w) / res[0]))
+ output_height = int(math.ceil((dst_n - dst_s) / res[1]))
+
+ # Adjust bounds to fit.
+ dst_e, dst_s = output_transform * (output_width, output_height)
+ logger.debug("Output width: %d, height: %d", output_width, output_height)
+ logger.debug("Adjusted bounds: %r", (dst_w, dst_s, dst_e, dst_n))
+
+ # create destination array
+ dest = np.zeros((first.count, output_height, output_width), dtype=dtype)
+
+ if nodata is not None:
+ nodataval = nodata
+ logger.debug("Set nodataval: %r", nodataval)
+
+ if nodataval is not None:
+ # Only fill if the nodataval is within dtype's range.
+ inrange = False
+ if np.dtype(dtype).kind in ('i', 'u'):
+ info = np.iinfo(dtype)
+ inrange = (info.min <= nodataval <= info.max)
+ elif np.dtype(dtype).kind == 'f':
+ info = np.finfo(dtype)
+ inrange = (info.min <= nodataval <= info.max)
+ if inrange:
+ dest.fill(nodataval)
+ else:
+ warnings.warn(
+ "Input file's nodata value, %s, is beyond the valid "
+ "range of its data type, %s. Consider overriding it "
+ "using the --nodata option for better results." % (
+ nodataval, dtype))
+ else:
+ nodataval = 0
+
+ for src in sources:
+ # Real World (tm) use of boundless reads.
+ # This approach uses the maximum amount of memory to solve the problem.
+ # Making it more efficient is a TODO.
+
+ # 1. Compute spatial intersection of destination and source.
+ src_w, src_s, src_e, src_n = src.bounds
+
+ int_w = src_w if src_w > dst_w else dst_w
+ int_s = src_s if src_s > dst_s else dst_s
+ int_e = src_e if src_e < dst_e else dst_e
+ int_n = src_n if src_n < dst_n else dst_n
+
+ # 2. Compute the source window.
+ src_window = get_window(
+ int_w, int_s, int_e, int_n, src.affine, precision=precision)
+ logger.debug("Src %s window: %r", src.name, src_window)
+
+ # 3. Compute the destination window.
+ dst_window = get_window(
+ int_w, int_s, int_e, int_n, output_transform, precision=precision)
+ logger.debug("Dst window: %r", dst_window)
+
+ # 4. Initialize temp array.
+ tcount = first.count
+ trows, tcols = tuple(b - a for a, b in dst_window)
+
+ temp_shape = (tcount, trows, tcols)
+ logger.debug("Temp shape: %r", temp_shape)
+
+ temp = np.zeros(temp_shape, dtype=dtype)
+ temp = src.read(out=temp, window=src_window, boundless=False,
+ masked=True)
+
+ # 5. Copy elements of temp into dest.
+ roff, coff = dst_window[0][0], dst_window[1][0]
+
+ region = dest[:, roff:roff + trows, coff:coff + tcols]
+ np.copyto(
+ region, temp,
+ where=np.logical_and(region==nodataval, temp.mask==False))
+
+ return dest, output_transform
diff --git a/rasterio/transform.py b/rasterio/transform.py
index c02bc1c..b832711 100644
--- a/rasterio/transform.py
+++ b/rasterio/transform.py
@@ -2,6 +2,7 @@ import warnings
from affine import Affine
+
IDENTITY = Affine.identity()
@@ -38,3 +39,11 @@ def from_bounds(west, south, east, north, width, height):
`height` in number of pixels."""
return Affine.translation(west, north) * Affine.scale(
(east - west)/width, (south - north)/height)
+
+
+def array_bounds(height, width, transform):
+ """Return the `west, south, east, north` bounds of an array given
+ its height, width, and an affine transform."""
+ w, n = transform.xoff, transform.yoff
+ e, s = transform * (width, height)
+ return w, s, e, n
diff --git a/requirements-dev.txt b/requirements-dev.txt
index a3e8d64..563cdc1 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,7 +1,7 @@
affine
cligj
coveralls>=0.4
-cython==0.22
+cython>=0.23.1
delocate
enum34
numpy>=1.8.0
diff --git a/setup.py b/setup.py
index c5b2c23..fe4cf4d 100755
--- a/setup.py
+++ b/setup.py
@@ -235,6 +235,7 @@ setup_args = dict(
[rasterio.rio_commands]
bounds=rasterio.rio.features:bounds
calc=rasterio.rio.calc:calc
+ clip=rasterio.rio.convert:clip
convert=rasterio.rio.convert:convert
edit-info=rasterio.rio.info:edit
env=rasterio.rio.info:env
diff --git a/tests/data/rgb1.tif b/tests/data/rgb1.tif
new file mode 100644
index 0000000..148b681
Binary files /dev/null and b/tests/data/rgb1.tif differ
diff --git a/tests/data/rgb2.tif b/tests/data/rgb2.tif
new file mode 100644
index 0000000..bc24905
Binary files /dev/null and b/tests/data/rgb2.tif differ
diff --git a/tests/data/rgb3.tif b/tests/data/rgb3.tif
new file mode 100644
index 0000000..ae89ad0
Binary files /dev/null and b/tests/data/rgb3.tif differ
diff --git a/tests/data/rgb4.tif b/tests/data/rgb4.tif
new file mode 100644
index 0000000..30b2c65
Binary files /dev/null and b/tests/data/rgb4.tif differ
diff --git a/tests/data/rgb_deflate.tif b/tests/data/rgb_deflate.tif
new file mode 100644
index 0000000..33eb8e5
Binary files /dev/null and b/tests/data/rgb_deflate.tif differ
diff --git a/tests/test_crs.py b/tests/test_crs.py
index 300b5cd..07b0e0c 100644
--- a/tests/test_crs.py
+++ b/tests/test_crs.py
@@ -89,6 +89,15 @@ def test_from_epsg_string():
assert crs.from_string('epsg:xyz')
+def test_from_string():
+ wgs84_crs = crs.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
+ assert wgs84_crs == {'no_defs': True, 'ellps': 'WGS84', 'datum': 'WGS84', 'proj': 'longlat'}
+
+ # Make sure this doesn't get handled using the from_epsg() even though 'epsg' is in the string
+ epsg_init_crs = crs.from_string('+units=m +init=epsg:26911 +no_defs=True')
+ assert epsg_init_crs == {'units': 'm', 'init': 'epsg:26911', 'no_defs': True}
+
+
def test_bare_parameters():
""" Make sure that bare parameters (e.g., no_defs) are handled properly,
even if they come in with key=True. This covers interaction with pyproj,
diff --git a/tests/test_image_structure.py b/tests/test_image_structure.py
new file mode 100644
index 0000000..ed43973
--- /dev/null
+++ b/tests/test_image_structure.py
@@ -0,0 +1,74 @@
+import rasterio
+
+from rasterio.enums import Compression, Interleaving
+
+
+def test_enum_compression_JPEG():
+ assert Compression('JPEG').name == 'jpeg'
+
+
+def test_enum_compression_LZW():
+ assert Compression('LZW').name == 'lzw'
+
+
+def test_enum_compression_PACKBITS():
+ assert Compression('PACKBITS').name == 'packbits'
+
+
+def test_enum_compression_DEFLATE():
+ assert Compression('DEFLATE').name == 'deflate'
+
+
+def test_enum_compression_CCITTRLE():
+ assert Compression('CCITTRLE').name == 'ccittrle'
+
+
+def test_enum_compression_CCITTFAX3():
+ assert Compression('CCITTFAX3').name == 'ccittfax3'
+
+
+def test_enum_compression_CCITTFAX4():
+ assert Compression('CCITTFAX4').name == 'ccittfax4'
+
+
+def test_enum_compression_LZMA():
+ assert Compression('LZMA').name == 'lzma'
+
+
+def test_enum_compression_NONE():
+ assert Compression('NONE').name == 'none'
+
+
+def test_compression_none():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ assert src.compression is None
+ assert 'compress' not in src.profile
+
+
+def test_compression_deflate():
+ with rasterio.open('tests/data/rgb_deflate.tif') as src:
+ assert src.compression.name == 'deflate'
+ assert src.compression.value == 'DEFLATE'
+ assert src.profile['compress'] == 'deflate'
+
+
+def test_enum_interleaving_BAND():
+ assert Interleaving('BAND').name == 'band'
+
+
+def test_enum_interleaving_PIXEL():
+ assert Interleaving('PIXEL').name == 'pixel'
+
+
+def test_interleaving_pixel():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ assert src.interleaving.name == 'pixel'
+ assert src.interleaving.value == 'PIXEL'
+ assert src.profile['interleave'] == 'pixel'
+
+
+def test_interleaving_pixel():
+ with rasterio.open('tests/data/rgb_deflate.tif') as src:
+ assert src.interleaving.name == 'band'
+ assert src.interleaving.value == 'BAND'
+ assert src.profile['interleave'] == 'band'
diff --git a/tests/test_indexing.py b/tests/test_indexing.py
index c1c5f73..257e7dc 100644
--- a/tests/test_indexing.py
+++ b/tests/test_indexing.py
@@ -24,7 +24,7 @@ def test_window_no_exception():
(0, src.height), (-4, src.width))
-def test_index():
+def test_index_values():
with rasterio.open('tests/data/RGB.byte.tif') as src:
assert src.index(101985.0, 2826915.0) == (0, 0)
assert src.index(101985.0+400.0, 2826915.0) == (0, 1)
@@ -41,11 +41,29 @@ def test_window():
(0, src.width))
assert src.index(left+400, top-400) == (1, 1)
assert src.index(left+dx+eps, top-dy-eps) == (1, 1)
- assert src.window(left, top-400, left+400, top) == ((0, 1), (0, 1))
- assert src.window(left, top-2*dy-eps, left+2*dx+eps, top) == ((0, 2), (0, 2))
+ assert src.window(left, top-400, left+400, top) == ((0, 2), (0, 2))
+ assert src.window(left, top-2*dy-eps, left+2*dx-eps, top) == ((0, 2), (0, 2))
def test_window_bounds_roundtrip():
with rasterio.open('tests/data/RGB.byte.tif') as src:
assert ((100, 200), (100, 200)) == src.window(
*src.window_bounds(((100, 200), (100, 200))))
+
+
+def test_window_full_cover():
+
+ def bound_covers(bounds1, bounds2):
+ """Does bounds1 cover bounds2?
+ """
+ return (bounds1[0] <= bounds2[0] and bounds1[1] <= bounds2[1] and
+ bounds1[2] >= bounds2[2] and bounds1[3] >= bounds2[3])
+
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ bounds = list(src.window_bounds(((100, 200), (100, 200))))
+ bounds[1] = bounds[1] - 10.0 # extend south
+ bounds[2] = bounds[2] + 10.0 # extend east
+
+ win = src.window(*bounds)
+ bounds_calc = list(src.window_bounds(win))
+ assert bound_covers(bounds_calc, bounds)
diff --git a/tests/test_read_boundless.py b/tests/test_read_boundless.py
index dc29c36..117479f 100644
--- a/tests/test_read_boundless.py
+++ b/tests/test_read_boundless.py
@@ -69,3 +69,18 @@ def test_read_boundless_masked_overlap():
assert not data.mask.all()
assert data.mask[0,399,399] == False
assert data.mask[0,0,0] == True
+
+
+def test_read_boundless_zero_stop():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ data = src.read(
+ window=((-200, 0), (-200, 0)), boundless=True, masked=True)
+ assert data.shape == (3, 200, 200)
+ assert data.mask.all()
+
+
+def test_read_boundless_masks_zero_stop():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ data = src.read_masks(window=((-200, 0), (-200, 0)), boundless=True)
+ assert data.shape == (3, 200, 200)
+ assert data.min() == data.max() == src.nodata
diff --git a/tests/test_read_resample.py b/tests/test_read_resample.py
index d21c0a4..5593cf7 100644
--- a/tests/test_read_resample.py
+++ b/tests/test_read_resample.py
@@ -10,16 +10,17 @@ import rasterio
def test_read_out_shape_resample_down():
with rasterio.open('tests/data/RGB.byte.tif') as s:
- out = numpy.zeros((7, 8), dtype=rasterio.ubyte)
+ out = numpy.zeros((8, 8), dtype=rasterio.ubyte)
data = s.read(1, out=out)
expected = numpy.array([
- [ 0, 8, 5, 7, 0, 0, 0, 0],
- [ 0, 6, 61, 15, 27, 15, 24, 128],
- [ 0, 20, 152, 23, 15, 19, 28, 0],
- [ 0, 17, 255, 25, 255, 22, 32, 0],
- [ 9, 7, 14, 16, 19, 18, 36, 0],
- [ 6, 27, 43, 207, 38, 31, 73, 0],
- [ 0, 0, 0, 0, 74, 23, 0, 0]], dtype=numpy.uint8)
+ [ 0, 0, 20, 15, 0, 0, 0, 0],
+ [ 0, 6, 193, 9, 255, 127, 23, 39],
+ [ 0, 7, 27, 255, 193, 14, 28, 34],
+ [ 0, 31, 29, 44, 14, 22, 43, 0],
+ [ 0, 9, 69, 49, 17, 22, 255, 0],
+ [ 11, 7, 13, 25, 13, 29, 33, 0],
+ [ 8, 10, 88, 27, 20, 33, 25, 0],
+ [ 0, 0, 0, 0, 98, 23, 0, 0]], dtype=numpy.uint8)
assert (data == expected).all() # all True.
diff --git a/tests/test_rio_convert.py b/tests/test_rio_convert.py
index 91590a2..e905180 100644
--- a/tests/test_rio_convert.py
+++ b/tests/test_rio_convert.py
@@ -1,15 +1,76 @@
import sys
import os
import logging
-import click
+import numpy
from click.testing import CliRunner
import rasterio
-from rasterio.rio.convert import convert
+from rasterio.rio.main import main_group
+from rasterio.rio.convert import convert, clip
logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+TEST_BBOX = [-11850000, 4804000, -11840000, 4808000]
+
+
+def test_clip_bounds(runner, tmpdir):
+ output = str(tmpdir.join('test.tif'))
+ result = runner.invoke(
+ main_group,
+ ['clip', 'tests/data/shade.tif', output, '--bounds'] + TEST_BBOX
+ )
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+
+ with rasterio.open(output) as out:
+ assert out.shape == (420, 173)
+
+
+def test_clip_like(runner, tmpdir):
+ output = str(tmpdir.join('test.tif'))
+ result = runner.invoke(
+ clip,
+ ['tests/data/shade.tif', output, '--like', 'tests/data/shade.tif']
+ )
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+
+ with rasterio.open('tests/data/shade.tif') as template_ds:
+ with rasterio.open(output) as out:
+ assert out.shape == template_ds.shape
+ assert numpy.allclose(out.bounds, template_ds.bounds)
+
+
+def test_clip_missing_params(runner, tmpdir):
+ output = str(tmpdir.join('test.tif'))
+ result = runner.invoke(
+ clip,
+ ['tests/data/shade.tif', output]
+ )
+ assert result.exit_code == 2
+ assert '--bounds or --like required' in result.output
+
+
+def test_clip_bounds_disjunct(runner, tmpdir):
+ output = str(tmpdir.join('test.tif'))
+ result = runner.invoke(
+ clip,
+ ['tests/data/shade.tif', output, '--bounds'] + [0, 0, 10, 10]
+ )
+ assert result.exit_code == 2
+ assert '--bounds' in result.output
+
+
+def test_clip_like_disjunct(runner, tmpdir):
+ output = str(tmpdir.join('test.tif'))
+ result = runner.invoke(
+ clip,
+ ['tests/data/shade.tif', output, '--like', 'tests/data/RGB.byte.tif']
+ )
+ assert result.exit_code == 2
+ assert '--like' in result.output
+
# Tests: format and type conversion, --format and --dtype
diff --git a/tests/test_rio_features.py b/tests/test_rio_features.py
index dd3da8e..a1a4b15 100644
--- a/tests/test_rio_features.py
+++ b/tests/test_rio_features.py
@@ -213,7 +213,7 @@ def test_mask_crop(runner, tmpdir):
with rasterio.open(output) as out:
assert out.shape[1] == src.shape[1]
assert out.shape[0] < src.shape[0]
- assert out.shape[0] == 823
+ assert out.shape[0] == 824
# Adding invert option after crop should be ignored
result = runner.invoke(
@@ -310,10 +310,10 @@ def test_shapes_compact(runner):
def test_shapes_sampling(runner):
result = runner.invoke(
- features.shapes, ['tests/data/shade.tif', '--sampling', '10'])
+ features.shapes, ['tests/data/shade.tif', '--sampling', '11'])
assert result.exit_code == 0
assert result.output.count('"FeatureCollection"') == 1
- assert result.output.count('"Feature"') == 124
+ assert result.output.count('"Feature"') == 117
def test_shapes_precision(runner):
diff --git a/tests/test_rio_merge.py b/tests/test_rio_merge.py
index d492373..db29e97 100644
--- a/tests/test_rio_merge.py
+++ b/tests/test_rio_merge.py
@@ -30,12 +30,12 @@ def test_data_dir_1(tmpdir):
with rasterio.drivers():
- with rasterio.open(str(tmpdir.join('a.tif')), 'w', **kwargs) as dst:
+ with rasterio.open(str(tmpdir.join('b.tif')), 'w', **kwargs) as dst:
data = numpy.ones((10, 10), dtype=rasterio.uint8)
data[0:6, 0:6] = 255
dst.write_band(1, data)
- with rasterio.open(str(tmpdir.join('b.tif')), 'w', **kwargs) as dst:
+ with rasterio.open(str(tmpdir.join('a.tif')), 'w', **kwargs) as dst:
data = numpy.ones((10, 10), dtype=rasterio.uint8)
data[4:8, 4:8] = 254
dst.write_band(1, data)
@@ -58,12 +58,12 @@ def test_data_dir_2(tmpdir):
with rasterio.drivers():
- with rasterio.open(str(tmpdir.join('a.tif')), 'w', **kwargs) as dst:
+ with rasterio.open(str(tmpdir.join('b.tif')), 'w', **kwargs) as dst:
data = numpy.zeros((10, 10), dtype=rasterio.uint8)
data[0:6, 0:6] = 255
dst.write_band(1, data)
- with rasterio.open(str(tmpdir.join('b.tif')), 'w', **kwargs) as dst:
+ with rasterio.open(str(tmpdir.join('a.tif')), 'w', **kwargs) as dst:
data = numpy.zeros((10, 10), dtype=rasterio.uint8)
data[4:8, 4:8] = 254
dst.write_band(1, data)
@@ -95,6 +95,7 @@ def test_merge_warn(test_data_dir_1):
runner = CliRunner()
result = runner.invoke(merge, inputs + [outputname] + ['--nodata', '-1'])
assert result.exit_code == 0
+ assert os.path.exists(outputname)
assert "using the --nodata option for better results" in result.output
@@ -130,12 +131,23 @@ def test_merge_output_exists(tmpdir):
assert out.count == 3
-def test_merge_output_exists_without_nodata(test_data_dir_2):
+def test_merge_output_exists_without_nodata_fails(test_data_dir_2):
+ """Fails without -f or --force-overwrite"""
runner = CliRunner()
result = runner.invoke(
merge,
[str(test_data_dir_2.join('a.tif')),
str(test_data_dir_2.join('b.tif'))])
+ assert result.exit_code == 1
+
+
+def test_merge_output_exists_without_nodata(test_data_dir_2):
+ """Succeeds with -f"""
+ runner = CliRunner()
+ result = runner.invoke(
+ merge,
+ ['-f', str(test_data_dir_2.join('a.tif')),
+ str(test_data_dir_2.join('b.tif'))])
assert result.exit_code == 0
@@ -173,12 +185,12 @@ def test_data_dir_overlapping(tmpdir):
}
with rasterio.drivers():
- with rasterio.open(str(tmpdir.join('nw.tif')), 'w', **kwargs) as dst:
+ with rasterio.open(str(tmpdir.join('se.tif')), 'w', **kwargs) as dst:
data = numpy.ones((10, 10), dtype=rasterio.uint8)
dst.write_band(1, data)
kwargs['transform'] = (-113, 0.2, 0, 45, 0, -0.2)
- with rasterio.open(str(tmpdir.join('se.tif')), 'w', **kwargs) as dst:
+ with rasterio.open(str(tmpdir.join('nw.tif')), 'w', **kwargs) as dst:
data = numpy.ones((10, 10), dtype=rasterio.uint8) * 2
dst.write_band(1, data)
@@ -219,12 +231,12 @@ def test_data_dir_float(tmpdir):
}
with rasterio.drivers():
- with rasterio.open(str(tmpdir.join('one.tif')), 'w', **kwargs) as dst:
+ with rasterio.open(str(tmpdir.join('two.tif')), 'w', **kwargs) as dst:
data = numpy.zeros((10, 10), dtype=rasterio.float64)
data[0:6, 0:6] = 255
dst.write_band(1, data)
- with rasterio.open(str(tmpdir.join('two.tif')), 'w', **kwargs) as dst:
+ with rasterio.open(str(tmpdir.join('one.tif')), 'w', **kwargs) as dst:
data = numpy.zeros((10, 10), dtype=rasterio.float64)
data[4:8, 4:8] = 254
dst.write_band(1, data)
@@ -295,15 +307,15 @@ def test_merge_tiny(tiffs):
# Output should be
#
- # [[ 0 120 120 90]
- # [ 0 120 120 90]
+ # [[ 0 120 90 90]
+ # [ 0 120 90 90]
# [ 0 60 0 0]
# [ 40 0 0 0]]
with rasterio.open(outputname) as src:
data = src.read()
- assert (data[0][0:2,1:3] == 120).all()
- assert (data[0][0:2,3] == 90).all()
+ assert (data[0][0:2,1] == 120).all()
+ assert (data[0][0:2,2:4] == 90).all()
assert data[0][2][1] == 60
assert data[0][3][0] == 40
@@ -318,34 +330,72 @@ def test_merge_tiny_output_opt(tiffs):
# Output should be
#
- # [[ 0 120 120 90]
- # [ 0 120 120 90]
+ # [[ 0 120 90 90]
+ # [ 0 120 90 90]
# [ 0 60 0 0]
# [ 40 0 0 0]]
with rasterio.open(outputname) as src:
data = src.read()
- assert (data[0][0:2,1:3] == 120).all()
- assert (data[0][0:2,3] == 90).all()
+ assert (data[0][0:2,1] == 120).all()
+ assert (data[0][0:2,2:4] == 90).all()
assert data[0][2][1] == 60
assert data[0][3][0] == 40
-def test_merge_tiny_res(tiffs):
+def test_merge_tiny_res_bounds(tiffs):
outputname = str(tiffs.join('merged.tif'))
inputs = [str(x) for x in tiffs.listdir()]
inputs.sort()
runner = CliRunner()
- result = runner.invoke(merge, inputs + [outputname, '--res', 2])
+ result = runner.invoke(merge, inputs + [outputname, '--res', 2, '--bounds', 1, 0, 5, 4])
assert result.exit_code == 0
# Output should be
# [[[120 90]
- # [ 0 0]]]
+ # [ 40 0]]]
with rasterio.open(outputname) as src:
data = src.read()
print(data)
assert data[0, 0, 0] == 120
assert data[0, 0, 1] == 90
- assert (data[0, 1, 0:1] == 0).all()
+ assert data[0, 1, 0] == 40
+ assert data[0, 1, 1] == 0
+
+
+def test_merge_tiny_res_high_precision(tiffs):
+ outputname = str(tiffs.join('merged.tif'))
+ inputs = [str(x) for x in tiffs.listdir()]
+ inputs.sort()
+ runner = CliRunner()
+ result = runner.invoke(merge, inputs + [outputname, '--res', 2, '--precision', 15])
+ assert result.exit_code == 0
+
+ # Output should be
+ # [[[120 90]
+ # [ 40 0]]]
+
+ with rasterio.open(outputname) as src:
+ data = src.read()
+ print(data)
+ assert data[0, 0, 0] == 120
+ assert data[0, 0, 1] == 90
+ assert data[0, 1, 0] == 40
+ assert data[0, 1, 1] == 0
+
+
+def test_merge_rgb(tmpdir):
+ """Get back original image"""
+ outputname = str(tmpdir.join('merged.tif'))
+ inputs = [
+ 'tests/data/rgb1.tif',
+ 'tests/data/rgb2.tif',
+ 'tests/data/rgb3.tif',
+ 'tests/data/rgb4.tif']
+ runner = CliRunner()
+ result = runner.invoke(merge, inputs + [outputname])
+ assert result.exit_code == 0
+
+ with rasterio.open(outputname) as src:
+ assert [src.checksum(i) for i in src.indexes] == [25420, 29131, 37860]
diff --git a/tests/test_tool.py b/tests/test_tool.py
index f579c92..b0e5285 100644
--- a/tests/test_tool.py
+++ b/tests/test_tool.py
@@ -6,7 +6,7 @@ except ImportError:
plt = None
import rasterio
-from rasterio.tool import show, stats
+from rasterio.tool import show, show_hist, stats
def test_stats():
@@ -42,3 +42,24 @@ def test_show():
show(src.read(1))
except ImportError:
pass
+
+
+def test_show_hist():
+ """
+ This test only verifies that code up to the point of plotting with
+ matplotlib works correctly. Tests do not exercise matplotlib.
+ """
+ if plt:
+ # Return because plotting causes the tests to block until the plot
+ # window is closed.
+ return
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ try:
+ show_hist((src, 1), bins=256)
+ except ImportError:
+ pass
+
+ try:
+ show_hist(src.read(), bins=256)
+ except ImportError:
+ pass
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/rasterio.git
More information about the Pkg-grass-devel
mailing list