[rasterio] 01/04: Imported Upstream version 0.24.0
Johan Van de Wauw
johanvdw-guest at moszumanska.debian.org
Tue Jun 2 20:23:23 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 eca81bae25b2e238dc9eb526285425a7b934ed5f
Author: Johan Van de Wauw <johan.vandewauw at gmail.com>
Date: Tue Jun 2 20:27:32 2015 +0200
Imported Upstream version 0.24.0
---
.travis.yml | 4 +-
AUTHORS.txt | 23 +-
CHANGES.txt | 48 ++++
README.rst | 19 ++
benchmarks/calc.sh | 21 ++
benchmarks/ndarray.py | 37 +--
docs/calc.rst | 119 +++++++++
docs/cli.rst | 117 ++++++++-
docs/img/mask_band.png | Bin 0 -> 24324 bytes
docs/img/mask_bands_rgb.png | Bin 0 -> 24606 bytes
docs/img/mask_conj.png | Bin 0 -> 24989 bytes
docs/img/mask_sieved.png | Bin 0 -> 21704 bytes
docs/masks.rst | 253 ++++++++++++-------
examples/reproject.py | 13 +-
rasterio/__init__.py | 12 +-
rasterio/_base.pyx | 79 +++++-
rasterio/_features.pyx | 11 +-
rasterio/_fill.pyx | 27 +-
rasterio/_gdal.pxd | 29 ++-
rasterio/_io.pxd | 2 +
rasterio/_io.pyx | 515 ++++++++++++++++++++++++++++++++-------
rasterio/_warp.pyx | 262 +++++++++++---------
rasterio/crs.py | 1 +
rasterio/features.py | 77 ++++--
rasterio/fill.py | 52 ++--
rasterio/profiles.py | 48 ++++
rasterio/rio/bands.py | 17 +-
rasterio/rio/calc.py | 146 +++++++++++
rasterio/rio/cli.py | 177 +++++++++++++-
rasterio/rio/features.py | 403 +++++++++++++++++++++++-------
rasterio/rio/info.py | 128 ++++++++--
rasterio/rio/main.py | 41 +++-
rasterio/rio/merge.py | 75 ++++--
rasterio/rio/rio.py | 33 +--
rasterio/rio/sample.py | 8 +-
rasterio/tool.py | 4 +-
rasterio/transform.py | 15 ++
rasterio/warp.py | 158 ++++++++++++
requirements-dev.txt | 3 +-
requirements.txt | 1 +
setup.py | 94 ++++---
tests/conftest.py | 15 +-
tests/test_band_masks.py | 100 ++++++++
tests/test_coords.py | 16 --
tests/test_crs.py | 54 ++++
tests/test_features_rasterize.py | 22 +-
tests/test_features_shapes.py | 36 ++-
tests/test_features_sieve.py | 50 +++-
tests/test_indexing.py | 31 ++-
tests/test_mask_creation.py | 45 ++++
tests/test_meta.py | 28 +++
tests/test_nodata.py | 5 +-
tests/test_profile.py | 79 ++++++
tests/test_read.py | 17 +-
tests/test_read_boundless.py | 20 ++
tests/test_read_resample.py | 2 +-
tests/test_rio_calc.py | 180 ++++++++++++++
tests/test_rio_cli.py | 18 ++
tests/test_rio_features.py | 280 ++++++++++++++++++---
tests/test_rio_info.py | 107 +++++++-
tests/test_rio_merge.py | 88 ++++++-
tests/test_rio_rio.py | 10 +-
tests/test_rio_sample.py | 12 +-
tests/test_sampling.py | 4 +-
tests/test_tool.py | 10 +-
tests/test_transform.py | 69 +++++-
tests/test_update.py | 82 +++++--
tests/test_warp.py | 371 ++++++++++++++++++++++++----
tests/test_write.py | 18 +-
69 files changed, 4048 insertions(+), 793 deletions(-)
diff --git a/.travis.yml b/.travis.yml
index 97cb081..440f792 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -10,10 +10,10 @@ before_install:
install:
- pip install --install-option="--no-cython-compile" cython
- "pip install -r requirements-dev.txt"
- - "pip install pytest"
+ - "pip install pytest pytest-cov"
- "pip install coveralls"
- "pip install -e ."
script:
- - coverage run --source=rasterio --omit='*.pxd,*.pyx,*/tests/*,*/docs/*,*/examples/*,*/benchmarks/*,*/rio/main.py,*/rio/__init__.py' -m py.test
+ - py.test --cov rasterio --cov-report term-missing
after_success:
- coveralls
diff --git a/AUTHORS.txt b/AUTHORS.txt
index 3cd8e34..cfbfc24 100644
--- a/AUTHORS.txt
+++ b/AUTHORS.txt
@@ -2,13 +2,20 @@ Authors
=======
Sean Gillies <sean at mapbox.com>
-Brendan Ward https://github.com/brendan-ward
-Asger Skovbo Petersen https://github.com/AsgerPetersen
-James Seppi https://github.com/jseppi
-Chrisophe Gohlke https://github.com/cgohlke
-Robin Wilson https://github.com/robintw
-Mike Toews https://github.com/mwtoews
-Amit Kapadia https://github.com/kapadia
-Alessandro Amici https://github.com/alexamici
+Brendan Ward <bcward at consbio.org>
+Ryan Grout <rgrout at continuum.io>
+Mike Toews <mwtoews at gmail.com>
+AsgerPetersen <asgerpetersen at gmail.com>
+Alessandro Amici <alexamici at gmail.com>
+Joshua Arnott <josh at snorfalorpagus.net>
+Amit Kapadia <amit at planet.com>
+Johan Van de Wauw <johan.vandewauw at gmail.com>
+Robin Wilson <robin at rtwilson.com>
+James Seppi <james.seppi at gmail.com>
+Etienne B. Racine <etiennebr at gmail.com>
+cgohlke <cgohlke at uci.edu>
+Kevin Wurster <wursterk at gmail.com>
+Aldo Culquicondor <alculquicondor at gmail.com>
+Martijn Visser <mgvisser at gmail.com>
See also https://github.com/mapbox/rasterio/graphs/contributors.
diff --git a/CHANGES.txt b/CHANGES.txt
index 9dce1f6..3ecfe74 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,6 +1,54 @@
Changes
=======
+0.24.0 (2015-05-27)
+-------------------
+- New rio-edit-info command (#358).
+- Add option to package GDAL data in distributions (#362).
+- Remove check that the path given to `rasterio.open()` in read mode is an
+ existing file, turning on some non-file formats (#364).
+- Addition of a `window_bounds()` method to dataset objects (#366).
+- Delegation of command exiting to Click (#367).
+
+0.23.0 (2015-05-08)
+-------------------
+- Redesign CLI as dynamically loaded entry points (#346).
+
+0.22.0 (2015-05-01)
+-------------------
+- Return masked arrays in the boundless read case (#338).
+- Add -o/--output option to rio-calc,merge,stack,mask,shapes,rasterize (#333).
+
+0.21.0 (2015-04-22)
+-------------------
+- New rio-mask command (#323).
+- Masking bug fix for rio-shapes (#335).
+- Addition of single valued nodata property to be used instead of nodatavals
+ (#329).
+
+0.20.0 (2015-04-08)
+-------------------
+- Switch read() default to masked=False (#300, #317).
+- Fix documentation of masking throughout module (#305).
+- Remove option for in place nodata filling (#309).
+- Enhancements for valid data footprint extraction in rio-shapes (#316, #318).
+
+0.19.1 (2015-03-30)
+-------------------
+- Add missing blockxsize, blockysize, tiled keywords (#301).
+
+0.19.0 (2015-03-25)
+-------------------
+- New rio-calc command (#175).
+- Added a file band shortcut to fillnodata() (#271).
+- Added fillnodata() to rio-calc functions (#277).
+- New approach to masking arrays on read that conforms more closely to GDAL's
+ RFC 15 (#282, #284, #285).
+- New read_masks() method (#284).
+- Deprecation of read_mask() and read_band (#284).
+- New affine transform factory functions from_origin(), from_bounds() (#287).
+- Improve correctness of indexing and rio-merge logic (#288, #290).
+
0.18.0 (2015-02-10)
-------------------
- New rio-rasterize command (#187).
diff --git a/README.rst b/README.rst
index 14dfc7e..5951445 100644
--- a/README.rst
+++ b/README.rst
@@ -236,6 +236,25 @@ Windows
Windows binary packages created by Christoph Gohlke are available `here
<http://www.lfd.uci.edu/~gohlke/pythonlibs/#rasterio>`_.
+You can download a binary distribution of GDAL from `here
+<http://www.gisinternals.com/release.php>`_. You will also need to download
+the compiled libraries and headers (include files).
+
+When building from source on Windows, it is important to know that setup.py
+cannot rely on gdal-config, which is only present on UNIX systems, to discover
+the locations of header files and libraries that rasterio needs to compile its
+C extensions. On Windows, these paths need to be provided by the user.
+You will need to find the include files and the library files for gdal and
+use setup.py as follows.
+
+.. code-block:: console
+
+ $ python setup.py build_ext -I<path to gdal include files> -lgdal_i -L<path to gdal library>
+ $ python setup.py install
+
+Note: The GDAL dll (gdal111.dll) and gdal-data directory need to be in your
+Windows PATH otherwise rasterio will fail to work.
+
Testing
-------
diff --git a/benchmarks/calc.sh b/benchmarks/calc.sh
new file mode 100755
index 0000000..04d06a0
--- /dev/null
+++ b/benchmarks/calc.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+echo "1. gdal_calc.py mult: 0.95 * a"
+echo "------------------------------"
+time gdal_calc.py --calc "0.95*A" -A tests/data/RGB.byte.tif --allBands A --overwrite --outfile out_gdal.tif
+echo ""
+
+echo "2. rio calc mult: 0.95 * a"
+echo "--------------------------"
+time rio calc "(* (read 1) 0.95)" tests/data/RGB.byte.tif out_rio.tif
+echo ""
+
+echo "3. gdal_calc.py mult add: 0.95 * a + 10"
+echo "---------------------------------------"
+time gdal_calc.py --calc "0.95*A + 10" -A tests/data/RGB.byte.tif --allBands A --overwrite --outfile out_gdal.tif
+echo ""
+
+echo "4. rio calc mult add: 0.95 * a + 10"
+echo "-----------------------------------"
+time rio calc "(+ (* (read 1) 0.95) 10)" tests/data/RGB.byte.tif out_rio.tif
+echo ""
diff --git a/benchmarks/ndarray.py b/benchmarks/ndarray.py
index b898138..da50d93 100644
--- a/benchmarks/ndarray.py
+++ b/benchmarks/ndarray.py
@@ -7,30 +7,30 @@ from osgeo import gdal
# GDAL
s = """
-src = gdal.Open('rasterio/tests/data/RGB.byte.tif')
+src = gdal.Open('tests/data/RGB.byte.tif')
arr = src.GetRasterBand(1).ReadAsArray()
src = None
"""
-n = 100
+n = 1000
t = timeit.timeit(s, setup='from osgeo import gdal', number=n)
print("GDAL:")
-print("%f msec\n" % (1000*t/n))
+print("%f usec\n" % (1000*t/n))
# Rasterio
s = """
-with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
- arr = src.read_band(1)
+with rasterio.open('tests/data/RGB.byte.tif') as src:
+ arr = src.read(1, masked=False)
"""
t = timeit.timeit(s, setup='import rasterio', number=n)
print("Rasterio:")
-print("%f msec\n" % (1000*t/n))
+print("%f usec\n" % (1000*t/n))
# GDAL Extras
s = """
-src = gdal.Open('rasterio/tests/data/RGB.byte.tif')
+src = gdal.Open('tests/data/RGB.byte.tif')
transform = src.GetGeoTransform()
srs = osr.SpatialReference()
srs.ImportFromWkt(src.GetProjectionRef())
@@ -44,30 +44,17 @@ n = 1000
t = timeit.timeit(s, setup='from osgeo import gdal; from osgeo import osr', number=n)
print("GDAL + Extras:\n")
-print("%f usec\n" % (t/n))
+print("%f usec\n" % (1000*t/n))
# Rasterio
s = """
-with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
- transform = src.transform
+with rasterio.open('tests/data/RGB.byte.tif') as src:
+ transform = src.affine
proj = src.crs
wkt = src.crs_wkt
- arr = src.read_band(1)
+ arr = src.read(1, masked=False)
"""
t = timeit.timeit(s, setup='import rasterio', number=n)
print("Rasterio:\n")
-print("%f usec\n" % (t/n))
-
-
-import pstats, cProfile
-
-s = """
-with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
- arr = src.read_band(1, window=(10, 10, 10, 10))
-"""
-
-cProfile.runctx(s, globals(), locals(), "Profile.prof")
-
-s = pstats.Stats("Profile.prof")
-s.strip_dirs().sort_stats("time").print_stats()
+print("%f usec\n" % (1000*t/n))
diff --git a/docs/calc.rst b/docs/calc.rst
new file mode 100644
index 0000000..7ee8c3b
--- /dev/null
+++ b/docs/calc.rst
@@ -0,0 +1,119 @@
+Using rio-calc
+==============
+
+Simple raster data processing on the command line is possible using Rasterio's
+rio-calc command. It uses the `snuggs <https://github.com/mapbox/snuggs>`__
+Numpy S-expression engine. The `snuggs README
+<https://github.com/mapbox/snuggs/blob/master/README.rst>`__ explains how
+expressions are written and evaluated in general. This document explains
+Rasterio-specific details of rio-calc and offers some examples.
+
+Expressions
+-----------
+
+Rio-calc expressions look like
+
+.. code-block::
+
+ (func|operator arg [*args])
+
+where ``func`` may be the name of any function in the module ``numpy`` or
+one of the rio-calc builtins: ``read``, ``fillnodata``, or ``sieve``; and
+``operator`` may be any of the standard Python arithmetic or logical operators.
+The arguments may themselves be expressions.
+
+Copying a file
+--------------
+
+Here's a trivial example of copying a dataset. The expression ``(read 1)``
+evaluates to all bands of the first input dataset, an array with shape
+``(3, 718, 791)`` in this case.
+
+Note: rio-calc's indexes start at ``1``.
+
+.. code-block:: console
+
+ $ rio calc "(read 1)" tests/data/RGB.byte.tif out.tif
+
+Reversing the band order of a file
+----------------------------------
+
+The expression ``(read i j)`` evaluates to the j-th band of the i-th input
+dataset. The ``asarray`` function collects bands read in reverse order into
+an array with shape ``(3, 718, 791)`` for output.
+
+.. code-block:: console
+
+ $ rio calc "(asarray (read 1 3) (read 1 2) (read 1 1))" tests/data/RGB.byte.tif out.tif
+
+Stacking bands of multiple files
+--------------------------------
+
+Bands can be read from multiple input files. This example is another (slower)
+way to copy a file.
+
+.. code-block:: console
+
+ $ rio calc "(asarray (read 1 1) (read 2 2) (read 3 3))" \
+ > tests/data/RGB.byte.tif tests/data/RGB.byte.tif tests/data/RGB.byte.tif \
+ > out.tif
+
+Named datasets
+--------------
+
+Datasets can be referenced in expressions by name and single bands picked out
+using the ``take`` function.
+
+.. code-block:: console
+
+ $ rio calc "(asarray (take a 3) (take a 2) (take a 1))" \
+ > --name "a=tests/data/RGB.byte.tif" out.tif
+
+The third example, re-done using names, is:
+
+.. code-block:: console
+
+ $ rio calc "(asarray (take a 1) (take b 2) (take b 3))" \
+ > --name "a=tests/data/RGB.byte.tif" --name "b=tests/data/RGB.byte.tif" \
+ > --name "c=tests/data/RGB.byte.tif" out.tif
+
+Read and take
+-------------
+
+The functions ``read`` and ``take`` overlap a bit in the previous examples but
+are rather different. The former involves I/O and the latter does not. You may
+also ``take`` from any array, as in this example.
+
+.. code-block:: console
+
+ $ rio calc "(take (read 1) 1)" tests/data/RGB.byte.tif out.tif
+
+Arithmetic operations
+---------------------
+
+Arithmetic operations can be performed as with Numpy. Here is an example of
+scaling all three bands of a dataset by the same factors.
+
+.. code-block:: console
+
+ $ rio calc "(+ 2 (* 0.95 (read 1)))" tests/data/RGB.byte.tif out.tif
+
+
+Here is a more complicated example of scaling bands by different factors.
+
+
+.. code-block:: console
+
+ $ rio calc "(asarray (+ 2 (* 0.95 (read 1 1))) (+ 3 (* 0.9 (read 1 2))) (+ 4 (* 0.85 (read 1 3))))" tests/data/RGB.byte.tif out.tif
+
+Logical operations
+------------------
+
+Logical operations can be used in conjunction with arithemtic operations. In
+this example, the output values are 255 wherever the input values are greater
+than or equal to 40.
+
+.. code-block:: console
+
+ $ rio calc "(* (>= (read 1) 40) 255)" tests/data/RGB.byte.tif out.tif
+
diff --git a/docs/cli.rst b/docs/cli.rst
index 4e6adf2..f32ab95 100644
--- a/docs/cli.rst
+++ b/docs/cli.rst
@@ -21,6 +21,7 @@ Rasterio's new command line interface is a program named "rio".
env Print information about the rio environment.
info Print information about a data file.
insp Open a data file and start an interpreter.
+ mask Mask in raster using features.
merge Merge a stack of raster datasets.
rasterize Rasterize features.
sample Sample a dataset.
@@ -29,6 +30,9 @@ Rasterio's new command line interface is a program named "rio".
It is developed using `Click <http://click.pocoo.org/3/>`__.
+Commands are shown below. See ``--help`` of individual commands for more
+details.
+
bounds
------
@@ -83,6 +87,115 @@ use with, e.g., `geojsonio-cli <https://github.com/mapbox/geojsonio-cli>`__.
Shoot the GeoJSON into a Leaflet map using geojsonio-cli by typing
``rio bounds tests/data/RGB.byte.tif | geojsonio``.
+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 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.
+
+Input files may have different numbers of bands but should have the same
+number of rows and columns. The output file will have the same number of rows
+and columns as the inputs and one band per element of the expression result.
+An expression involving arithmetic operations on N-D arrays will produce a
+N-D array and result in an N-band output file.
+
+The following produces a 3-band GeoTIFF with all values scaled by 0.95 and
+incremented by 2. In the expression, ``(read 1)`` evaluates to the first
+input dataset (3 bands) as a 3-D array.
+
+.. code-block:: console
+
+ $ rio calc "(+ 2 (* 0.95 (read 1)))" tests/data/RGB.byte.tif /tmp/out.tif
+
+The following produces a 3-band GeoTIFF in which the first band is copied from
+the first band of the input and the next two bands are scaled (down) by the
+ratio of the first band's mean to their own means. The ``--name`` option is
+used to bind datasets to a name within the expression. ``(take a 1)`` gets the
+first band of the dataset named ``a`` as a 2-D array and ``(asarray ...)``
+collects a sequence of 2-D arrays into a 3-D array for output.
+
+.. code-block:: console
+
+ $ 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 /tmp/out.rgb.tif
+
+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.
+
+Please see `calc.rst <calc.rst>`__ for more details.
+
+
+edit-info
+---------
+
+New in 0.24
+
+The edit-info command allows you edit a raster dataset's metadata, namely
+
+- coordinate reference system
+- affine transformation matrix
+- nodata value
+- tags
+
+A TIFF created by spatially-unaware image processing software like Photoshop
+or Imagemagick can be turned into a GeoTIFF by editing these metadata items.
+
+You can set or change a dataset's coordinate reference system to, e.g.,
+EPSG:3857 (Web Mercator),
+
+.. code-block:: console
+
+ $ rio edit-info --crs EPSG:3857 example.tif
+
+set its `affine transformation matrix <https://github.com/mapbox/rasterio/blob/master/docs/georeferencing.rst#coordinate-transformation>`__,
+
+.. code-block:: console
+
+ $ rio edit-info --transform "[300.0, 0.0, 101985.0, 0.0, -300.0, 2826915.0]"
+
+or set its nodata value to, e.g., `0`:
+
+.. code-block:: console
+
+ $ rio edit-info --nodata 0 example.tif
+
+
+mask
+----
+
+New in 0.21
+
+The mask command masks in pixels from all bands of a raster using features
+(masking out all areas not covered by features) and optionally crops the output
+raster to the extent of the features. Features are assumed to be in the same
+coordinate reference system as the input raster.
+
+A common use case is masking in raster data by political or other boundaries.
+
+.. code-block:: console
+
+ $ rio mask input.tif output.tif --geojson-mask input.geojson
+
+GeoJSON features may be provided using stdin or specified directly as first
+argument, and output can be cropped to the extent of the features.
+
+.. code-block:: console
+
+ $ rio mask input.tif output.tif --crop --geojson-mask - < input.geojson
+
+The feature mask can be inverted to mask out pixels covered by features and
+keep pixels not covered by features.
+
+.. code-block:: console
+
+ $ rio mask input.tif output.tif --invert --geojson-mask input.geojson
+
+
info
----
@@ -322,7 +435,7 @@ following.
.. code-block:: console
- $ echo "[-78.0, 23.0]" | rio transform - --dst_crs EPSG:32618 --precision 2
+ $ echo "[-78.0, 23.0]" | rio transform - --dst-crs EPSG:32618 --precision 2
[192457.13, 2546667.68]
To transform a longitude, latitude bounding box to the coordinate system of
@@ -330,7 +443,7 @@ a raster dataset, do the following.
.. code-block:: console
- $ echo "[-78.0, 23.0, -76.0, 25.0]" | rio transform - --dst_crs tests/data/RGB.byte.tif --precision 2
+ $ echo "[-78.0, 23.0, -76.0, 25.0]" | rio transform - --dst-crs tests/data/RGB.byte.tif --precision 2
[192457.13, 2546667.68, 399086.97, 2765319.94]
Suggestions for other commands are welcome!
diff --git a/docs/img/mask_band.png b/docs/img/mask_band.png
new file mode 100644
index 0000000..246ef06
Binary files /dev/null and b/docs/img/mask_band.png differ
diff --git a/docs/img/mask_bands_rgb.png b/docs/img/mask_bands_rgb.png
new file mode 100644
index 0000000..35ab7f7
Binary files /dev/null and b/docs/img/mask_bands_rgb.png differ
diff --git a/docs/img/mask_conj.png b/docs/img/mask_conj.png
new file mode 100644
index 0000000..7658edc
Binary files /dev/null and b/docs/img/mask_conj.png differ
diff --git a/docs/img/mask_sieved.png b/docs/img/mask_sieved.png
new file mode 100644
index 0000000..78681b5
Binary files /dev/null and b/docs/img/mask_sieved.png differ
diff --git a/docs/masks.rst b/docs/masks.rst
index 10a6239..fb4da60 100644
--- a/docs/masks.rst
+++ b/docs/masks.rst
@@ -1,122 +1,205 @@
Masks
=====
-Reading masks
--------------
+In using Rasterio, you'll encounter two different kinds of masks. One is the
+the valid data mask from GDAL, an unsigned byte array with the same number of
+rows and columns as the dataset in which non-zero elements indicate that the
+corresponding data elements are valid. Other elements are invalid, or *nodata*
+elements. The other kind of mask is the mask in Numpy's [masked
+arrays](http://docs.scipy.org/doc/numpy/reference/maskedarray.generic.html),
+which have the inverse sense: `True` values in a masked array's mask indicate
+that the corresponding data elements are invalid. With care, you can safely
+navigate this divide.
+
+Consider Rasterio's RGB.byte.tif test dataset. It has 718 rows and 791
+columns of pixels. Each pixel has 3 8-bit (uint8) channels or bands. It has a
+trapezoid of image data within a rectangular background of 0,0,0 value pixels.
+
+.. image:: https://www.dropbox.com/s/sg7qejccih5m4ah/RGB.byte.jpg?dl=1
+
+Metadata in the dataset declares that values of 0 shall be interpreted as
+invalid data or *nodata* pixels. In, e.g., merging the image with adjacent
+scenes, we'd like to ignore the nodata pixels and have only valid image data in
+our final mosaic.
+
+Let's use the rio-insp command to look at the two kinds of masks and their
+inverse relationship in the context of RGB.byte.tif.
+
+.. code-block:: console
+
+ $ rio insp tests/data/RGB.byte.tif
+ Rasterio 0.19.0 Interactive Inspector (Python 2.7.9)
+ Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
+ >>> src.shape
+ (718, 791)
+ >>> src.count
+ 3
+ >>> src.dtypes
+ ['uint8', 'uint8', 'uint8']
+ >>> src.nodatavals
+ [0.0, 0.0, 0.0]
+
+Reading dataset masks
+---------------------
+
+For every band of a dataset there is a mask. These masks can be had as arrays
+using the dataset's `read_masks()`` method. Below, ``msk`` is the valid data
+mask corresponding to the first dataset band.
-There are a few different ways for raster datasets to carry valid data masks.
-Rasterio subscribes to GDAL's abstract mask band interface, so although the
-module's main test dataset has no mask band, GDAL will build one based upon
-its declared nodata value of 0.
+.. code-block:: python
+
+ >>> msk = src.read_masks(1)
+ >>> msk
+ array([[0, 0, 0, ..., 0, 0, 0],
+ [0, 0, 0, ..., 0, 0, 0],
+ [0, 0, 0, ..., 0, 0, 0],
+ ...,
+ [0, 0, 0, ..., 0, 0, 0],
+ [0, 0, 0, ..., 0, 0, 0],
+ [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)
+
+This array is a valid data mask in the sense of `GDAL RFC 15
+<https://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask>`__. The 0 values in its
+corners represent *nodata* regions. Zooming in on the interior of the mask
+array shows the ``255`` values that indicate valid data regions.
.. code-block:: python
-
- import rasterio
- with rasterio.open('tests/data/RGB.byte.tif') as src:
- mask = src.read_mask()
- print mask.any()
- count = mask.shape[0] * mask.shape[1]
- print float((mask > 0).sum())/count
- print float((mask == 0).sum())/count
+ >>> m[200:250,200:250]
+ array([[255, 255, 255, ..., 255, 255, 255],
+ [255, 255, 255, ..., 255, 255, 255],
+ [255, 255, 255, ..., 255, 255, 255],
+ ...,
+ [255, 255, 255, ..., 255, 255, 255],
+ [255, 255, 255, ..., 255, 255, 255],
+ [255, 255, 255, ..., 255, 255, 255]], dtype=uint8)
-Some of the elements of the mask evaluate to ``True``, meaning that there is some
-valid data. Just over 2/3 of the dataset's pixels (use of sum being a neat trick for
-computing the number of pixels in a selection) have valid data.
+Displayed using Matplotlib's `imshow()`, the mask looks like this:
-.. code-block:: console
+.. image:: img/mask_band.png
- True
- 0.673974976142
- 0.326025023858
+Wait, what are these 0 values in the mask interior? This is an example of
+a problem inherent in 8-bit raster data: lack of dynamic range. The dataset
+creator has said that 0 values represent missing data (see the
+``nodatavals`` property in the first code block of this document), but some of
+the valid data have values so low they've been rounded during processing to
+zero. This can very easily happen in scaling 16-bit data to 8 bits. There's
+no magic nodata value bullet for this. Using 16 bits per band helps, but you
+really have to be careful with 8-bit per band datasets and their nodata values.
Writing masks
-------------
-Writing a mask is just as straightforward: pass an ndarray with ``True`` (or values
-that evaluate to ``True`` to indicate valid data and ``False`` to indicate no data
-to ``write_mask()``.
+Writing a mask that applies to all dataset bands is just as straightforward:
+pass an ndarray with ``True`` (or values that evaluate to ``True`` to indicate
+valid data and ``False`` to indicate no data to ``write_mask()``. Consider a
+copy of the test data opened using rio-insp in "r+" (update) mode.
.. code-block:: python
- import os
- import shutil
- import tempfile
+ $ rio insp copy.tif --mode r+
+ Rasterio 0.19.0 Interactive Inspector (Python 2.7.9)
+ Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
+ >>>
- import numpy
- import rasterio
+To mark that all pixels of all bands are valid (i.e., to override nodata
+metadata values that can't be unset), you'd do this.
- tempdir = tempfile.mkdtemp()
+.. code-block::
- with rasterio.open(
- os.path.join(tempdir, 'example.tif'),
- 'w',
- driver='GTiff',
- count=1,
- dtype=rasterio.uint8,
- width=10,
- height=10) as dst:
-
- dst.write_band(1, numpy.ones(dst.shape, dtype=rasterio.uint8))
+ >>> src.write_mask(True)
+ >>> src.read_masks(1).all()
+ True
+
+No data have been altered, nor have the dataset's nodata values been changed.
+A new band has been added to the dataset to store the valid data mask. By
+default it is saved to a "sidecar" GeoTIFF alongside the dataset file. When
+such a .msk GeoTIFF exists, Rasterio will ignore the nodata metadata values and
+return mask arrays based on the .msk file.
- mask = numpy.zeros(dst.shape, rasterio.uint8)
- mask[5:,5:] = 255
- dst.write_mask(mask)
+.. code-block:: console
- print os.listdir(tempdir)
- shutil.rmtree(tempdir)
+ $ ls -l copy.tif*
+ -rw-r--r--@ 1 sean staff 1713704 Mar 24 14:19 copy.tif
+ -rw-r--r-- 1 sean staff 916 Mar 24 14:25 copy.tif.msk
-The code above masks out all of the file except the lower right quadrant and
-writes a file with a sidecar TIFF to hold the mask.
+Can Rasterio help fix buggy nodata masks like the ones in RGB.byte.tif? It
+certainly can. Consider a fresh copy of that file. This time we'll read all
+3 band masks (based on the nodata values, not a .msk GeoTIFF) and show them
+as an RGB image (with the help of `numpy.dstack()`):
-.. code-block:: console
+.. code-block:: python
+
+ $rio insp copy.tif --mode r+
+ Rasterio 0.19.0 Interactive Inspector (Python 2.7.9)
+ Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
+ >>> msk = src.read_masks()
+ >>> show(np.dstack(msk))
- ['example.tif', 'example.tif.msk']
+.. image:: img/mask_bands_rgb.png
-To use an internal TIFF mask, use the ``drivers()`` option shown below:
+Colored regions appear where valid data pixels don't quite coincide. This is,
+again, an artifact of scaling data down to 8 bits per band. We'll begin by
+constructing a new mask array from the logical conjunction of the three band
+masks we've read.
.. code-block:: python
- tempdir = tempfile.mkdtemp()
- tiffname = os.path.join(tempdir, 'example.tif')
+ >>> new_msk = (msk[0] & msk[1] & msk[2])
+ >>> show(new_msk)
- with rasterio.drivers(GDAL_TIFF_INTERNAL_MASK=True):
+.. image:: img/mask_conj.png
- with rasterio.open(
- tiffname,
- 'w',
- driver='GTiff',
- count=1,
- dtype=rasterio.uint8,
- width=10,
- height=10) as dst:
-
- dst.write_band(1, numpy.ones(dst.shape, dtype=rasterio.uint8))
+Now we'll use `sieve()` to shake out the small buggy regions of the mask. I've
+found the right value for the ``size`` argument empirically.
- mask = numpy.zeros(dst.shape, rasterio.uint8)
- mask[5:,5:] = 255
- dst.write_mask(mask)
+.. code-block:: python
- print os.listdir(tempdir)
- print subprocess.check_output(['gdalinfo', tiffname])
+ >>> from rasterio.features import sieve
+ >>> sieved_msk = sieve(new_msk, size=800)
+ >>> show(sieved_msk)
-The output:
+.. image:: img/mask_sieved.png
-.. code-block:: console
+Last thing to do is write that sieved mask back to the dataset.
- ['example.tif']
- Driver: GTiff/GeoTIFF
- Files: /var/folders/jh/w0mgrfqd1t37n0bcqzt16bnc0000gn/T/tmpcnGV_r/example.tif
- Size is 10, 10
- Coordinate System is `'
- Image Structure Metadata:
- INTERLEAVE=BAND
- Corner Coordinates:
- Upper Left ( 0.0, 0.0)
- Lower Left ( 0.0, 10.0)
- Upper Right ( 10.0, 0.0)
- Lower Right ( 10.0, 10.0)
- Center ( 5.0, 5.0)
- Band 1 Block=10x10 Type=Byte, ColorInterp=Gray
- Mask Flags: PER_DATASET
+.. code-block:: python
+
+ >>> src.write_mask(sieved_msk)
+
+The result is a properly masked dataset that allows some 0 value pixels to be
+considered valid.
+
+Numpy masked arrays
+-------------------
+
+If you want, you can read dataset bands as numpy masked arrays.
+
+.. code-block:: python
+
+ $ rio insp tests/data/RGB.byte.tif
+ Rasterio 0.19.0 Interactive Inspector (Python 2.7.9)
+ Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
+ >>> blue = src.read(1, masked=True)
+ >>> blue.mask
+ array([[ True, True, True, ..., True, True, True],
+ [ True, True, True, ..., True, True, True],
+ [ True, True, True, ..., True, True, True],
+ ...,
+ [ True, True, True, ..., True, True, True],
+ [ True, True, True, ..., True, True, True],
+ [ True, True, True, ..., True, True, True]], dtype=bool)
+
+As mentioned earlier, this mask is the inverse of the GDAL band mask. To get
+a mask conforming to GDAL RFC 15, simply do this:
+.. code-block:: python
+
+ >>> msk = (~blue.mask * 255).astype('uint8')
+
+You can rely on this Rasterio identity for any integer value ``N``.
+
+.. code-block:: python
+
+ >>> (~src.read(N, masked=True).mask * 255 == src.read_masks(N)).all()
+ True
diff --git a/examples/reproject.py b/examples/reproject.py
index 9a8a348..89e6985 100644
--- a/examples/reproject.py
+++ b/examples/reproject.py
@@ -5,7 +5,7 @@ import tempfile
import numpy
import rasterio
-from rasterio import Affine as A
+from rasterio import transform
from rasterio.warp import reproject, RESAMPLING
tempdir = '/tmp'
@@ -17,17 +17,15 @@ with rasterio.drivers():
# with each pixel covering 15".
rows, cols = src_shape = (512, 512)
dpp = 1.0/240 # decimal degrees per pixel
- # The following is equivalent to
- # A(dpp, 0, -cols*dpp/2, 0, -dpp, rows*dpp/2).
- src_transform = A.translation(-cols*dpp/2, rows*dpp/2) * A.scale(dpp, -dpp)
+ west, south, east, north = -cols*dpp/2, -rows*dpp/2, cols*dpp/2, rows*dpp/2
+ src_transform = transform.from_bounds(west, south, east, north, cols, rows)
src_crs = {'init': 'EPSG:4326'}
source = numpy.ones(src_shape, numpy.uint8)*255
# Prepare to reproject this rasters to a 1024 x 1024 dataset in
- # Web Mercator (EPSG:3857) with origin at -8928592, 2999585.
+ # Web Mercator (EPSG:3857) with origin at -237481.5, 237536.4.
dst_shape = (1024, 1024)
- dst_transform = A.from_gdal(-237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0)
- dst_transform = dst_transform.to_gdal()
+ dst_transform = transform.from_origin(-237481.5, 237536.4, 425.0, 425.0)
dst_crs = {'init': 'EPSG:3857'}
destination = numpy.zeros(dst_shape, numpy.uint8)
@@ -59,4 +57,3 @@ with rasterio.drivers():
dst.write_band(1, destination)
info = subprocess.call(['open', tiffname])
-
diff --git a/rasterio/__init__.py b/rasterio/__init__.py
index 0ea3dd2..5b0f009 100644
--- a/rasterio/__init__.py
+++ b/rasterio/__init__.py
@@ -12,13 +12,18 @@ from rasterio.dtypes import (
bool_, ubyte, uint8, uint16, int16, uint32, int32, float32, float64,
complex_)
from rasterio.five import string_types
+from rasterio.profiles import default_gtiff_profile
from rasterio.transform import Affine, guard_transform
+# These modules are imported from the Cython extensions, but are also import
+# here to help tools like cx_Freeze find them automatically
+from rasterio import _err, coords, enums
+
# Classes in rasterio._io are imported below just before we need them.
__all__ = [
'band', 'open', 'drivers', 'copy', 'pad']
-__version__ = "0.18"
+__version__ = "0.24.0"
log = logging.getLogger('rasterio')
class NullHandler(logging.Handler):
@@ -84,12 +89,9 @@ def open(
raise TypeError("invalid mode: %r" % mode)
if driver and not isinstance(driver, string_types):
raise TypeError("invalid driver: %r" % driver)
- if mode in ('r', 'r+'):
- if not os.path.exists(path):
- raise IOError("no such file or directory: %r" % path)
if transform:
transform = guard_transform(transform)
-
+
if mode == 'r':
from rasterio._io import RasterReader
s = RasterReader(path)
diff --git a/rasterio/_base.pyx b/rasterio/_base.pyx
index 96c7540..0cd52b2 100644
--- a/rasterio/_base.pyx
+++ b/rasterio/_base.pyx
@@ -276,8 +276,9 @@ cdef class DatasetReader(object):
def get_nodatavals(self):
cdef void *hband = NULL
- cdef object val
+ cdef double nodataval
cdef int success
+
if not self._nodatavals:
if self._hds == NULL:
raise ValueError("can't read closed raster file")
@@ -285,16 +286,24 @@ cdef class DatasetReader(object):
hband = _gdal.GDALGetRasterBand(self._hds, i+1)
if hband == NULL:
raise ValueError("Null band")
- val = _gdal.GDALGetRasterNoDataValue(hband, &success)
+ nodataval = _gdal.GDALGetRasterNoDataValue(hband, &success)
+ val = nodataval
if not success:
val = None
self._nodatavals.append(val)
return self._nodatavals
property nodatavals:
+ """Nodata values for each band."""
+
def __get__(self):
return self.get_nodatavals()
+ property nodata:
+ """The dataset's single nodata value."""
+ def __get__(self):
+ return self.nodatavals[0]
+
def block_windows(self, bidx=0):
"""Returns an iterator over a band's block windows and their
indexes.
@@ -371,31 +380,45 @@ cdef class DatasetReader(object):
def index(self, x, y):
"""Returns the (row, col) index of the pixel containing (x, y)."""
a, b, c, d, e, f, _, _, _ = self.affine
- return int(round((y-f)/e)), int(round((x-c)/a))
+ return int(math.floor((y-f)/e)), int(math.floor((x-c)/a))
+
+ 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."""
- def window(self, left, bottom, right, top):
- """Returns the window corresponding to the world bounding box."""
- ul = self.index(left, top)
- lr = self.index(right, bottom)
- return tuple(zip(ul, lr))
+ window = tuple(zip(self.index(left, top), self.index(right, bottom)))
+ if boundless:
+ return window
+ else:
+ return crop_window(window, self.height, self.width)
def window_transform(self, window):
"""Returns the affine transform for a dataset window."""
(r, _), (c, _) = window
return self.affine * Affine.translation(c or 0, r or 0)
+ def window_bounds(self, window):
+ """Returns the bounds of a window as x_min, y_min, x_max, y_max."""
+ ((row_min, row_max), (col_min, col_max)) = window
+ x_min, y_min = (col_min, row_max) * self.affine
+ x_max, y_max = (col_max, row_min) * self.affine
+ return x_min, y_min, x_max, y_max
+
@property
def meta(self):
m = {
'driver': self.driver,
'dtype': self.dtypes[0],
- 'nodata': self.nodatavals[0],
+ 'nodata': self.nodata,
'width': self.width,
'height': self.height,
'count': self.count,
'crs': self.crs,
'transform': self.affine.to_gdal(),
- 'affine': self.affine }
+ 'affine': self.affine,
+ 'blockxsize': self.block_shapes[0][1],
+ 'blockysize': self.block_shapes[0][0],
+ 'tiled': self.block_shapes[0][1] != self.width }
self._read = True
return m
@@ -563,6 +586,16 @@ cdef class DatasetReader(object):
# A window is a 2D ndarray indexer in the form of a tuple:
# ((row_start, row_stop), (col_start, col_stop))
+cpdef crop_window(object window, int height, int width):
+ """Returns a window cropped to fall within height and width."""
+ cdef int r_start, r_stop, c_start, c_stop
+ (r_start, r_stop), (c_start, c_stop) = window
+ return (
+ (min(max(r_start, 0), height), max(0, min(r_stop, height))),
+ (min(max(c_start, 0), width), max(0, min(c_stop, width)))
+ )
+
+
cpdef eval_window(object window, int height, int width):
"""Evaluates a window tuple that might contain negative values
in the context of a raster height and width."""
@@ -709,3 +742,29 @@ def _transform(src_crs, dst_crs, xs, ys, zs):
_gdal.OSRDestroySpatialReference(src)
_gdal.OSRDestroySpatialReference(dst)
return retval
+
+
+def is_geographic_crs(crs):
+ cdef void *osr_crs = _osr_from_crs(crs)
+ cdef int retval = _gdal.OSRIsGeographic(osr_crs)
+ _gdal.OSRDestroySpatialReference(osr_crs)
+
+ return retval == 1
+
+
+def is_projected_crs(crs):
+ cdef void *osr_crs = _osr_from_crs(crs)
+ cdef int retval = _gdal.OSRIsProjected(osr_crs)
+ _gdal.OSRDestroySpatialReference(osr_crs)
+
+ return retval == 1
+
+
+def is_same_crs(crs1, crs2):
+ cdef void *osr_crs1 = _osr_from_crs(crs1)
+ cdef void *osr_crs2 = _osr_from_crs(crs2)
+ cdef int retval = _gdal.OSRIsSame(osr_crs1, osr_crs2)
+ _gdal.OSRDestroySpatialReference(osr_crs1)
+ _gdal.OSRDestroySpatialReference(osr_crs2)
+
+ return retval == 1
diff --git a/rasterio/_features.pyx b/rasterio/_features.pyx
index c7c63c6..bc83a65 100644
--- a/rasterio/_features.pyx
+++ b/rasterio/_features.pyx
@@ -29,8 +29,8 @@ def _shapes(image, mask, connectivity, transform):
Data type must be one of rasterio.int16, rasterio.int32,
rasterio.uint8, rasterio.uint16, or rasterio.float32.
mask : numpy ndarray or rasterio Band object
- Values of False will be excluded from feature generation
- Must be of type rasterio.bool_
+ Values of False or 0 will be excluded from feature generation
+ Must evaluate to bool (rasterio.bool_ or rasterio.uint8)
connectivity : int
Use 4 or 8 pixel connectivity for grouping pixels into features
transform : Affine transformation
@@ -150,8 +150,8 @@ def _sieve(image, size, output, mask, connectivity):
output : numpy ndarray
Array of same shape and data type as `image` in which to store results.
mask : numpy ndarray or rasterio Band object
- Values of False will be excluded from feature generation.
- Must be of type rasterio.bool_.
+ Values of False or 0 will be excluded from feature generation.
+ Must evaluate to bool (rasterio.bool_ or rasterio.uint8)
connectivity : int
Use 4 or 8 pixel connectivity for grouping pixels into features.
@@ -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/_fill.pyx b/rasterio/_fill.pyx
index 8d7c998..2515564 100644
--- a/rasterio/_fill.pyx
+++ b/rasterio/_fill.pyx
@@ -15,14 +15,16 @@ from rasterio._io cimport InMemoryRaster
def _fillnodata(image, mask, double max_search_distance=100.0,
int smoothing_iterations=0):
cdef void *memdriver = _gdal.GDALGetDriverByName("MEM")
- cdef void *image_dataset
- cdef void *image_band
- cdef void *mask_dataset
- cdef void *mask_band
+ cdef void *image_dataset = NULL
+ cdef void *image_band = NULL
+ cdef void *mask_dataset = NULL
+ cdef void *mask_band = NULL
+ cdef _io.RasterReader rdr
+ cdef _io.RasterReader mrdr
cdef char **alg_options = NULL
if isinstance(image, np.ndarray):
- # copy numpy ndarray into an in-memory dataset
+ # copy numpy ndarray into an in-memory dataset.
image_dataset = _gdal.GDALCreate(
memdriver,
"image",
@@ -33,9 +35,6 @@ def _fillnodata(image, mask, double max_search_distance=100.0,
NULL)
image_band = _gdal.GDALGetRasterBand(image_dataset, 1)
_io.io_auto(image, image_band, True)
- elif isinstance(image, tuple):
- # TODO
- raise NotImplementedError()
else:
raise ValueError("Invalid source image")
@@ -52,8 +51,10 @@ def _fillnodata(image, mask, double max_search_distance=100.0,
mask_band = _gdal.GDALGetRasterBand(mask_dataset, 1)
_io.io_auto(mask_cast, mask_band, True)
elif isinstance(mask, tuple):
- # TODO
- raise NotImplementedError()
+ if mask.shape != image.shape:
+ raise ValueError("Mask must have same shape as image")
+ mrdr = mask.ds
+ mask_band = mrdr.band(mask.bidx)
elif mask is None:
mask_band = NULL
else:
@@ -77,8 +78,10 @@ def _fillnodata(image, mask, double max_search_distance=100.0,
result = np.empty(image.shape, dtype=image.dtype)
_io.io_auto(result, image_band, False)
- _gdal.GDALClose(image_dataset)
- _gdal.GDALClose(mask_dataset)
+ if image_dataset != NULL:
+ _gdal.GDALClose(image_dataset)
+ if mask_dataset != NULL:
+ _gdal.GDALClose(mask_dataset)
_gdal.CSLDestroy(alg_options)
return result
diff --git a/rasterio/_gdal.pxd b/rasterio/_gdal.pxd
index 0ae46f6..bc42d2e 100644
--- a/rasterio/_gdal.pxd
+++ b/rasterio/_gdal.pxd
@@ -17,24 +17,35 @@ cdef extern from "cpl_string.h":
char ** CSLSetNameValue (char **list, char *name, char *val)
void CSLDestroy (char **list)
+cdef extern from "cpl_vsi.h":
+ ctypedef int vsi_l_offset
+ unsigned char * VSIGetMemFileBuffer (const char *filename,
+ vsi_l_offset *data_len,
+ int take_ownership)
+
cdef extern from "ogr_srs_api.h":
+ void * OCTNewCoordinateTransformation (void *source, void *dest)
+ void OCTDestroyCoordinateTransformation (void *source)
+ int OCTTransform (void *ct, int nCount, double *x, double *y, double *z)
+
+ int OSRAutoIdentifyEPSG (void *srs)
void OSRCleanup ()
void * OSRClone (void *srs)
void OSRDestroySpatialReference (void *srs)
int OSRExportToProj4 (void *srs, char **params)
int OSRExportToWkt (void *srs, char **params)
- int OSRImportFromEPSG (void *srs, int code)
- int OSRImportFromProj4 (void *srs, char *proj)
- int OSRSetFromUserInput (void *srs, char *input)
- int OSRAutoIdentifyEPSG (void *srs)
int OSRFixup(void *srs)
const char * OSRGetAuthorityName (void *srs, const char *key)
const char * OSRGetAuthorityCode (void *srs, const char *key)
+ int OSRImportFromEPSG (void *srs, int code)
+ int OSRImportFromProj4 (void *srs, char *proj)
+ int OSRIsGeographic(void *srs)
+ int OSRIsProjected(void *srs)
+ int OSRIsSame(void *srs1, void *srs2)
void * OSRNewSpatialReference (char *wkt)
void OSRRelease (void *srs)
- void * OCTNewCoordinateTransformation (void *source, void *dest)
- void OCTDestroyCoordinateTransformation (void *source)
- int OCTTransform (void *ct, int nCount, double *x, double *y, double *z)
+ int OSRSetFromUserInput (void *srs, char *input)
+
cdef extern from "gdal.h" nogil:
void GDALAllRegister()
@@ -82,8 +93,7 @@ cdef extern from "gdal.h" nogil:
int GDALSetRasterNoDataValue(void *band, double value)
int GDALDatasetRasterIO(void *band, int, int xoff, int yoff, int xsize, int ysize, void *buffer, int width, int height, int, int count, int *bmap, int poff, int loff, int boff)
int GDALRasterIO(void *band, int, int xoff, int yoff, int xsize, int ysize, void *buffer, int width, int height, int, int poff, int loff)
-
- int GDALSetRasterNoDataValue(void *band, double value)
+ int GDALFillRaster(void *band, double rvalue, double ivalue)
void * GDALCreate(void *driver, const char *filename, int width, int height, int nbands, GDALDataType dtype, const char **options)
void * GDALCreateCopy(void *driver, const char *filename, void *ds, int strict, char **options, void *progress_func, void *progress_data)
@@ -111,6 +121,7 @@ cdef extern from "gdal.h" nogil:
int GDALGetRasterColorInterpretation (void *hBand)
int GDALSetRasterColorInterpretation (void *hBand, int)
+ int GDALGetMaskFlags (void *hBand)
void *GDALGetMaskBand (void *hBand)
int GDALCreateMaskBand (void *hDS, int flags)
diff --git a/rasterio/_io.pxd b/rasterio/_io.pxd
index bc866ea..9e5165f 100644
--- a/rasterio/_io.pxd
+++ b/rasterio/_io.pxd
@@ -61,6 +61,8 @@ ctypedef np.int32_t DTYPE_INT32_t
ctypedef np.float32_t DTYPE_FLOAT32_t
ctypedef np.float64_t DTYPE_FLOAT64_t
+cdef bint in_dtype_range(value, dtype)
+
cdef int io_ubyte(
void *hband,
int mode,
diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx
index e9339a5..64ef9ec 100644
--- a/rasterio/_io.pyx
+++ b/rasterio/_io.pyx
@@ -13,7 +13,7 @@ cimport numpy as np
from rasterio cimport _base, _gdal, _ogr, _io
from rasterio._base import (
- eval_window, window_shape, window_index, tastes_like_gdal)
+ crop_window, eval_window, window_shape, window_index, tastes_like_gdal)
from rasterio._drivers import driver_count, GDALEnv
from rasterio._err import cpl_errs
from rasterio import dtypes
@@ -33,7 +33,14 @@ else:
def emit(self, record):
pass
- log.addHandler(NullHandler())
+log.addHandler(NullHandler())
+
+cdef bint in_dtype_range(value, dtype):
+ """Returns True if value is in the range of dtype, else False."""
+ infos = {
+ 'c': np.finfo, 'f': np.finfo, 'i': np.iinfo, 'u': np.iinfo}
+ rng = infos[np.dtype(dtype).kind](dtype)
+ return rng.min <= value <= rng.max
# Single band IO functions.
@@ -316,11 +323,12 @@ cdef int io_multi_cint16(
cdef int i, j, k
cdef np.int16_t real, imag
- buf = np.empty(
+ buf = np.zeros(
(out.shape[0], 2*out.shape[2]*out.shape[1]),
dtype=np.int16)
cdef np.int16_t[:, :] buf_view = buf
+
with nogil:
bandmap = <int *>_gdal.CPLMalloc(count*sizeof(int))
for i in range(count):
@@ -492,47 +500,113 @@ cdef int io_multi_cfloat64(
return retval
+cdef int io_multi_mask(
+ void *hds,
+ int mode,
+ int xoff,
+ int yoff,
+ int width,
+ int height,
+ np.uint8_t[:, :, :] buffer,
+ long[:] indexes,
+ int count):
+ cdef int i, j, retval=0
+ cdef void *hband
+ cdef void *hmask
+
+ for i in range(count):
+ j = indexes[i]
+ hband = _gdal.GDALGetRasterBand(hds, j)
+ if hband == NULL:
+ raise ValueError("Null band")
+ hmask = _gdal.GDALGetMaskBand(hband)
+ if hmask == NULL:
+ raise ValueError("Null mask band")
+ with nogil:
+ retval = _gdal.GDALRasterIO(
+ hmask, mode, xoff, yoff, width, height,
+ &buffer[i, 0, 0], buffer.shape[2], buffer.shape[1], 1, 0, 0)
+ if retval:
+ break
+ return retval
+
+
cdef int io_auto(image, void *hband, bint write):
"""
Convenience function to handle IO with a GDAL band and a 2D numpy image
:param image: a numpy 2D image
:param hband: an instance of GDALGetRasterBand
- :param write: 1 (True) uses write mode, 0 (False) uses read
+ :param write: 1 (True) uses write mode (writes image into band),
+ 0 (False) uses read mode (reads band into image)
:return: the return value from the data-type specific IO function
"""
- if not len(image.shape) == 2:
- raise ValueError("Specified image must have 2 dimensions")
-
- cdef int width = image.shape[1]
- cdef int height = image.shape[0]
-
+ cdef int ndims = len(image.shape)
+ cdef int height = image.shape[-2]
+ cdef int width = image.shape[-1]
+ cdef int count
+ cdef long[:] indexes
dtype_name = image.dtype.name
- if dtype_name == "float32":
- return io_float32(hband, write, 0, 0, width, height, image)
- elif dtype_name == "float64":
- return io_float64(hband, write, 0, 0, width, height, image)
- elif dtype_name == "uint8":
- return io_ubyte(hband, write, 0, 0, width, height, image)
- elif dtype_name == "int16":
- return io_int16(hband, write, 0, 0, width, height, image)
- elif dtype_name == "int32":
- return io_int32(hband, write, 0, 0, width, height, image)
- elif dtype_name == "uint16":
- return io_uint16(hband, write, 0, 0, width, height, image)
- elif dtype_name == "uint32":
- return io_uint32(hband, write, 0, 0, width, height, image)
+ if ndims == 2:
+ if dtype_name == "float32":
+ return io_float32(hband, write, 0, 0, width, height, image)
+ elif dtype_name == "float64":
+ return io_float64(hband, write, 0, 0, width, height, image)
+ elif dtype_name == "uint8":
+ return io_ubyte(hband, write, 0, 0, width, height, image)
+ elif dtype_name == "int16":
+ return io_int16(hband, write, 0, 0, width, height, image)
+ elif dtype_name == "int32":
+ return io_int32(hband, write, 0, 0, width, height, image)
+ elif dtype_name == "uint16":
+ return io_uint16(hband, write, 0, 0, width, height, image)
+ elif dtype_name == "uint32":
+ return io_uint32(hband, write, 0, 0, width, height, image)
+ else:
+ raise ValueError("Image dtype is not supported for this function."
+ "Must be float32, float64, int16, int32, uint8, "
+ "uint16, or uint32")
+ elif ndims == 3:
+ count = image.shape[0]
+ indexes = np.arange(1, count + 1)
+
+ dtype_name = image.dtype.name
+
+ if dtype_name == "float32":
+ return io_multi_float32(hband, write, 0, 0, width, height, image,
+ indexes, count)
+ elif dtype_name == "float64":
+ return io_multi_float64(hband, write, 0, 0, width, height, image,
+ indexes, count)
+ elif dtype_name == "uint8":
+ return io_multi_ubyte(hband, write, 0, 0, width, height, image,
+ indexes, count)
+ elif dtype_name == "int16":
+ return io_multi_int16(hband, write, 0, 0, width, height, image,
+ indexes, count)
+ elif dtype_name == "int32":
+ return io_multi_int32(hband, write, 0, 0, width, height, image,
+ indexes, count)
+ elif dtype_name == "uint16":
+ return io_multi_uint16(hband, write, 0, 0, width, height, image,
+ indexes, count)
+ elif dtype_name == "uint32":
+ return io_multi_uint32(hband, write, 0, 0, width, height, image,
+ indexes, count)
+ else:
+ raise ValueError("Image dtype is not supported for this function."
+ "Must be float32, float64, int16, int32, uint8, "
+ "uint16, or uint32")
+
else:
- raise ValueError("Image dtype is not supported for this function."
- "Must be float32, float64, int16, int32, uint8, "
- "uint16, or uint32")
+ raise ValueError("Specified image must have 2 or 3 dimensions")
cdef class RasterReader(_base.DatasetReader):
- def read_band(self, bidx, out=None, window=None, masked=None):
+ def read_band(self, bidx, out=None, window=None, masked=False):
"""Read the `bidx` band into an `out` array if provided,
otherwise return a new array.
@@ -545,10 +619,15 @@ cdef class RasterReader(_base.DatasetReader):
example, ((0, 2), (0, 2)) defines a 2x2 window at the upper left
of the raster dataset.
"""
+ warnings.warn(
+ "read_band() is deprecated and will be removed by Rasterio 1.0. "
+ "Please use read() instead.",
+ FutureWarning,
+ stacklevel=2)
return self.read(bidx, out=out, window=window, masked=masked)
- def read(self, indexes=None, out=None, window=None, masked=None,
+ def read(self, indexes=None, out=None, window=None, masked=False,
boundless=False):
"""Read raster bands as a multidimensional array
@@ -576,12 +655,10 @@ cdef class RasterReader(_base.DatasetReader):
a 2x2 window at the upper left of the raster dataset.
masked : bool, optional
- The return type will be either a regular NumPy array, or
- a masked NumPy array depending on the `masked` argument. The
- return type is forced if either `True` or `False`, but will
- be chosen if `None`. For `masked=None` (default), the array
- will be the same type as `out` (if used), or will be masked
- if any of the nodatavals are not `None`.
+ If `masked` is `True` the return value will be a masked
+ array. Otherwise (the default) the return value will be a
+ regular array. Masks will be exactly the inverse of the
+ GDAL RFC 15 conforming arrays returned by read_masks().
boundless : bool, optional (default `False`)
If `True`, windows that extend beyond the dataset's extent
@@ -597,6 +674,8 @@ cdef class RasterReader(_base.DatasetReader):
preferentially used by callers.
"""
+ cdef void *hband = NULL
+
return2d = False
if indexes is None:
indexes = self.indexes
@@ -615,8 +694,30 @@ cdef class RasterReader(_base.DatasetReader):
if bidx not in self.indexes:
raise IndexError("band index out of range")
idx = self.indexes.index(bidx)
- check_dtypes.add(self.dtypes[idx])
- nodatavals.append(self._nodatavals[idx])
+
+ dtype = self.dtypes[idx]
+ check_dtypes.add(dtype)
+
+ ndv = self._nodatavals[idx]
+ # Change given nodatavals to the closest value that
+ # can be represented by this band's data type to
+ # match GDAL's strategy.
+ if ndv is not None:
+ if np.dtype(dtype).kind in ('i', 'u'):
+ info = np.iinfo(dtype)
+ dt_min, dt_max = info.min, info.max
+ elif np.dtype(dtype).kind in ('f', 'c'):
+ info = np.finfo(dtype)
+ dt_min, dt_max = info.min, info.max
+ else:
+ dt_min, dt_max = False, True
+ if ndv < dt_min:
+ ndv = dt_min
+ elif ndv > dt_max:
+ ndv = dt_max
+
+ nodatavals.append(ndv)
+
# Mixed dtype reads are not supported at this time.
if len(check_dtypes) > 1:
raise ValueError("more than one 'dtype' found")
@@ -632,13 +733,12 @@ cdef class RasterReader(_base.DatasetReader):
win_shape += (
window[0][1]-window[0][0], window[1][1]-window[1][0])
else:
- w = eval_window(window, self.height, self.width)
- minr = min(max(w[0][0], 0), self.height)
- maxr = max(0, min(w[0][1], self.height))
- minc = min(max(w[1][0], 0), self.width)
- maxc = max(0, min(w[1][1], self.width))
- win_shape += (maxr - minr, maxc - minc)
- window = ((minr, maxr), (minc, maxc))
+ window = crop_window(
+ eval_window(window, self.height, self.width),
+ self.height, self.width
+ )
+ (r_start, r_stop), (c_start, c_stop) = window
+ win_shape += (r_stop - r_start, c_stop - c_start)
else:
win_shape += self.shape
@@ -651,14 +751,33 @@ cdef class RasterReader(_base.DatasetReader):
raise ValueError(
"'out' shape %s does not match window shape %s" %
(out.shape, win_shape))
- if masked is None:
- masked = hasattr(out, 'mask')
- if masked is None:
- masked = any([x is not None for x in nodatavals])
+
+ # Masking
+ # -------
+ #
+ # If masked is True, we check the GDAL mask flags using
+ # GDALGetMaskFlags. If GMF_ALL_VALID for all bands, we do not
+ # call read_masks(), but pass `mask=False` to the masked array
+ # constructor. Else, we read the GDAL mask bands using
+ # read_masks(), invert them and use them in constructing masked
+ # arrays.
+
+ if masked:
+
+ mask_flags = [0]*self.count
+ for i, j in zip(range(self.count), self.indexes):
+ hband = _gdal.GDALGetRasterBand(self._hds, j)
+ mask_flags[i] = _gdal.GDALGetMaskFlags(hband)
+
+ all_valid = all([flag & 0x01 == 1 for flag in mask_flags])
+
+ log.debug("all_valid: %s", all_valid)
+ log.debug("mask_flags: %r", mask_flags)
+
if out is None:
out = np.zeros(win_shape, dtype)
for ndv, arr in zip(
- self.nodatavals, out if len(win_shape) == 3 else [out]):
+ nodatavals, out if len(out.shape) == 3 else [out]):
if ndv is not None:
arr.fill(ndv)
@@ -667,6 +786,23 @@ cdef class RasterReader(_base.DatasetReader):
if not boundless or not window:
out = self._read(indexes, out, window, dtype)
+ if masked:
+ if all_valid:
+ mask = np.ma.nomask
+ else:
+ mask = np.empty(out.shape, 'uint8')
+ mask = ~self._read(
+ indexes, mask, window, 'uint8', masks=True
+ ).astype('bool')
+
+ kwds = {'mask': mask}
+ # Set a fill value only if the read bands share a
+ # single nodata value.
+ if len(set(nodatavals)) == 1:
+ if nodatavals[0] is not None:
+ kwds['fill_value'] = nodatavals[0]
+ out = np.ma.array(out, **kwds)
+
else:
# Compute the overlap between the dataset and the boundless window.
overlap = ((
@@ -682,9 +818,171 @@ cdef class RasterReader(_base.DatasetReader):
overlap_w = overlap[1][1] - overlap[1][0]
scaling_h = float(out.shape[-2:][0])/window_h
scaling_w = float(out.shape[-2:][1])/window_w
- buffer_shape = (int(overlap_h*scaling_h), int(overlap_w*scaling_w))
+ buffer_shape = (
+ int(round(overlap_h*scaling_h)),
+ int(round(overlap_w*scaling_w)))
data = np.empty(win_shape[:-2] + buffer_shape, dtype)
data = self._read(indexes, data, overlap, dtype)
+
+ if masked:
+ mask = np.empty(win_shape[:-2] + buffer_shape, 'uint8')
+ mask = ~self._read(
+ indexes, mask, overlap, 'uint8', masks=True
+ ).astype('bool')
+ kwds = {'mask': mask}
+ if len(set(nodatavals)) == 1:
+ if nodatavals[0] is not None:
+ kwds['fill_value'] = nodatavals[0]
+ data = np.ma.array(data, **kwds)
+
+ else:
+ data = None
+ if masked:
+ kwds = {'mask': True}
+ if len(set(nodatavals)) == 1:
+ if nodatavals[0] is not None:
+ kwds['fill_value'] = nodatavals[0]
+ out = np.ma.array(out, **kwds)
+
+ if data is not None:
+ # Determine where to put the data in the output window.
+ data_h, data_w = buffer_shape
+ roff = 0
+ coff = 0
+ if window[0][0] < 0:
+ roff = int(round(window_h*scaling_h)) - data_h
+ if window[1][0] < 0:
+ coff = int(round(window_w*scaling_w)) - data_w
+
+ for dst, src in zip(
+ out if len(out.shape) == 3 else [out],
+ data if len(data.shape) == 3 else [data]):
+ dst[roff:roff+data_h, coff:coff+data_w] = src
+
+ if masked:
+ if not hasattr(out, 'mask'):
+ kwds = {'mask': True}
+ if len(set(nodatavals)) == 1:
+ if nodatavals[0] is not None:
+ kwds['fill_value'] = nodatavals[0]
+ out = np.ma.array(out, **kwds)
+
+ for dst, src in zip(
+ out.mask if len(out.shape) == 3 else [out.mask],
+ data.mask if len(data.shape) == 3 else [data.mask]):
+ dst[roff:roff+data_h, coff:coff+data_w] = src
+
+ if return2d:
+ out.shape = out.shape[1:]
+
+ return out
+
+
+ def read_masks(self, indexes=None, out=None, window=None, boundless=False):
+ """Read raster band masks as a multidimensional array
+
+ Parameters
+ ----------
+ indexes : list of ints or a single int, optional
+ If `indexes` is a list, the result is a 3D array, but is
+ a 2D array if it is a band index number.
+
+ out: numpy ndarray, optional
+ As with Numpy ufuncs, this is an optional reference to an
+ output array with the same dimensions and shape into which
+ data will be placed.
+
+ *Note*: the method's return value may be a view on this
+ array. In other words, `out` is likely to be an
+ incomplete representation of the method's results.
+
+ window : a pair (tuple) of pairs of ints, optional
+ The optional `window` argument is a 2 item tuple. The first
+ item is a tuple containing the indexes of the rows at which
+ the window starts and stops and the second is a tuple
+ containing the indexes of the columns at which the window
+ starts and stops. For example, ((0, 2), (0, 2)) defines
+ a 2x2 window at the upper left of the raster dataset.
+
+ boundless : bool, optional (default `False`)
+ If `True`, windows that extend beyond the dataset's extent
+ are permitted and partially or completely filled arrays will
+ be returned as appropriate.
+
+ Returns
+ -------
+ Numpy ndarray or a view on a Numpy ndarray
+
+ Note: as with Numpy ufuncs, an object is returned even if you
+ use the optional `out` argument and the return value shall be
+ preferentially used by callers.
+ """
+
+ return2d = False
+ if indexes is None:
+ indexes = self.indexes
+ elif isinstance(indexes, int):
+ indexes = [indexes]
+ return2d = True
+ if out is not None and out.ndim == 2:
+ out.shape = (1,) + out.shape
+ if not indexes:
+ raise ValueError("No indexes to read")
+
+ # Get the natural shape of the read window, boundless or not.
+ win_shape = (len(indexes),)
+ if window:
+ if boundless:
+ win_shape += (
+ window[0][1]-window[0][0], window[1][1]-window[1][0])
+ else:
+ w = eval_window(window, self.height, self.width)
+ minr = min(max(w[0][0], 0), self.height)
+ maxr = max(0, min(w[0][1], self.height))
+ minc = min(max(w[1][0], 0), self.width)
+ maxc = max(0, min(w[1][1], self.width))
+ win_shape += (maxr - minr, maxc - minc)
+ window = ((minr, maxr), (minc, maxc))
+ else:
+ win_shape += self.shape
+
+ dtype = 'uint8'
+
+ if out is not None:
+ if out.dtype != np.dtype(dtype):
+ raise ValueError(
+ "the out array's dtype '%s' does not match '%s'"
+ % (out.dtype, dtype))
+ if out.shape[0] != win_shape[0]:
+ raise ValueError(
+ "'out' shape %s does not match window shape %s" %
+ (out.shape, win_shape))
+ if out is None:
+ out = np.zeros(win_shape, 'uint8')
+
+ # We can jump straight to _read() in some cases. We can ignore
+ # the boundless flag if there's no given window.
+ if not boundless or not window:
+ out = self._read(indexes, out, window, dtype, masks=True)
+
+ 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)))
+
+ if overlap != ((0, 0), (0, 0)):
+ # Prepare a buffer.
+ window_h, window_w = win_shape[-2:]
+ overlap_h = overlap[0][1] - overlap[0][0]
+ overlap_w = overlap[1][1] - overlap[1][0]
+ scaling_h = float(out.shape[-2:][0])/window_h
+ scaling_w = float(out.shape[-2:][1])/window_w
+ buffer_shape = (int(overlap_h*scaling_h), int(overlap_w*scaling_w))
+ data = np.empty(win_shape[:-2] + buffer_shape, 'uint8')
+ data = self._read(indexes, data, overlap, dtype, masks=True)
else:
data = None
@@ -702,32 +1000,13 @@ cdef class RasterReader(_base.DatasetReader):
data if len(data.shape) == 3 else [data]):
dst[roff:roff+data_h, coff:coff+data_w] = src
- # Masking the output. TODO: explain the logic better.
- if masked:
- if len(set(nodatavals)) == 1:
- if nodatavals[0] is None:
- out = np.ma.masked_array(out, copy=False)
- elif np.isnan(nodatavals[0]):
- out = np.ma.masked_where(np.isnan(out), out, copy=False)
- else:
- out = np.ma.masked_equal(out, nodatavals[0], copy=False)
- else:
- out = np.ma.masked_array(out, copy=False)
- for aix in range(len(indexes)):
- if nodatavals[aix] is None:
- band_mask = False
- elif np.isnan(nodatavals[aix]):
- band_mask = np.isnan(out[aix])
- else:
- band_mask = out[aix] == nodatavals[aix]
- out[aix].mask = band_mask
if return2d:
out.shape = out.shape[1:]
return out
- def _read(self, indexes, out, window, dtype):
+ def _read(self, indexes, out, window, dtype, masks=False):
"""Read raster bands as a multidimensional array
If `indexes` is a list, the result is a 3D array, but
@@ -769,7 +1048,11 @@ cdef class RasterReader(_base.DatasetReader):
indexes_count = <int>indexes_arr.shape[0]
gdt = dtypes.dtype_rev[dtype]
- if gdt == 1:
+ if masks:
+ retval = io_multi_mask(
+ self._hds, 0, xoff, yoff, width, height,
+ out, indexes_arr, indexes_count)
+ elif gdt == 1:
retval = io_multi_ubyte(
self._hds, 0, xoff, yoff, width, height,
out, indexes_arr, indexes_count)
@@ -822,7 +1105,7 @@ cdef class RasterReader(_base.DatasetReader):
return out
- def read_mask(self, out=None, window=None):
+ def read_mask(self, indexes=None, out=None, window=None, boundless=False):
"""Read the mask band into an `out` array if provided,
otherwise return a new array containing the dataset's
valid data mask.
@@ -835,6 +1118,13 @@ cdef class RasterReader(_base.DatasetReader):
"""
cdef void *hband
cdef void *hmask
+
+ warnings.warn(
+ "read_mask() is deprecated and will be removed by Rasterio 1.0. "
+ "Please use read_masks() instead.",
+ FutureWarning,
+ stacklevel=2)
+
if self._hds == NULL:
raise ValueError("can't write closed raster file")
hband = _gdal.GDALGetRasterBand(self._hds, 1)
@@ -859,13 +1149,16 @@ cdef class RasterReader(_base.DatasetReader):
xoff = yoff = 0
width = self.width
height = self.height
- retval = io_ubyte(
+
+ io_ubyte(
hmask, 0, xoff, yoff, width, height, out)
return out
def sample(self, xy, indexes=None):
"""Get the values of a dataset at certain positions
+ Values are from the nearest pixel. They are not interpolated.
+
Parameters
----------
xy : iterable, pairs of floats
@@ -889,7 +1182,6 @@ cdef class RasterReader(_base.DatasetReader):
cdef class RasterUpdater(RasterReader):
# Read-write access to raster data and metadata.
- # TODO: the r+ mode.
def __init__(
self, path, mode, driver=None,
@@ -1012,6 +1304,13 @@ cdef class RasterUpdater(RasterReader):
raise ValueError("NULL dataset")
if self._init_nodata is not None:
+
+ if not in_dtype_range(self._init_nodata, self._init_dtype):
+ raise ValueError(
+ "Given nodata value, %s, is beyond the valid "
+ "range of its data type, %s." % (
+ self._init_nodata, self._init_dtype))
+
for i in range(self._count):
hband = _gdal.GDALGetRasterBand(self._hds, i+1)
success = _gdal.GDALSetRasterNoDataValue(
@@ -1141,12 +1440,15 @@ cdef class RasterUpdater(RasterReader):
def set_nodatavals(self, vals):
cdef void *hband = NULL
- cdef double val
+ cdef double nodataval
+ cdef int success
+
for i, val in zip(self.indexes, vals):
hband = _gdal.GDALGetRasterBand(self._hds, i)
- success = _gdal.GDALSetRasterNoDataValue(hband, val)
+ nodataval = val
+ success = _gdal.GDALSetRasterNoDataValue(hband, nodataval)
if success:
- raise ValueError("Invalid nodata values")
+ raise ValueError("Invalid nodata value: %r", val)
self._nodatavals = vals
property nodatavals:
@@ -1157,8 +1459,21 @@ cdef class RasterUpdater(RasterReader):
return self.get_nodatavals()
def __set__(self, value):
+ warnings.warn(
+ "nodatavals.__set__() is deprecated and will be removed by "
+ "Rasterio 1.0. Please use nodata.__set__() instead.",
+ FutureWarning,
+ stacklevel=2)
self.set_nodatavals(value)
+ property nodata:
+ """The dataset's single nodata value."""
+
+ def __get__(self):
+ return self.nodatavals[0]
+
+ def __set__(self, value):
+ self.set_nodatavals([value for old_val in self.nodatavals])
def write(self, src, indexes=None, window=None):
"""Write the src array into indexed bands of the dataset.
@@ -1366,7 +1681,7 @@ cdef class RasterUpdater(RasterReader):
_gdal.GDALSetRasterColorTable(hBand, hTable)
_gdal.GDALDestroyColorTable(hTable)
- def write_mask(self, src, window=None):
+ def write_mask(self, mask, window=None):
"""Write the valid data mask src array into the dataset's band
mask.
@@ -1389,6 +1704,7 @@ cdef class RasterUpdater(RasterReader):
if hmask == NULL:
raise ValueError("NULL band mask")
log.debug("Created mask band")
+
if window:
window = eval_window(window, self.height, self.width)
yoff = window[0][0]
@@ -1399,14 +1715,19 @@ cdef class RasterUpdater(RasterReader):
xoff = yoff = 0
width = self.width
height = self.height
- if src.dtype == np.bool:
- array = 255 * src.astype(np.uint8)
+
+ if mask is True:
+ _gdal.GDALFillRaster(hmask, 255, 0)
+ elif mask is False:
+ _gdal.GDALFillRaster(hmask, 0, 0)
+ elif mask.dtype == np.bool:
+ array = 255 * mask.astype(np.uint8)
retval = io_ubyte(
hmask, 1, xoff, yoff, width, height, array)
else:
retval = io_ubyte(
- hmask, 1, xoff, yoff, width, height, src)
-
+ hmask, 1, xoff, yoff, width, height, mask)
+
cdef class InMemoryRaster:
"""
@@ -1647,14 +1968,30 @@ def writer(path, mode, **kwargs):
fname = name_b
with cpl_errs:
hds = _gdal.GDALOpen(fname, 0)
- if hds == NULL:
- raise ValueError("NULL dataset")
- drv = _gdal.GDALGetDatasetDriver(hds)
- drv_name = _gdal.GDALGetDriverShortName(drv)
- drv_name_b = drv_name
- driver = drv_name_b.decode('utf-8')
- _gdal.GDALClose(hds)
+ if hds == NULL:
+ raise ValueError("NULL dataset")
+ drv = _gdal.GDALGetDatasetDriver(hds)
+ drv_name = _gdal.GDALGetDriverShortName(drv)
+ drv_name_b = drv_name
+ driver = drv_name_b.decode('utf-8')
+ _gdal.GDALClose(hds)
+
if driver == 'GTiff':
return RasterUpdater(path, mode)
else:
return IndirectRasterUpdater(path, mode)
+
+
+def virtual_file_to_buffer(filename):
+ """Read content of a virtual file into a Python bytes buffer."""
+ cdef unsigned char *buff = NULL
+ cdef const char *cfilename = NULL
+ cdef _gdal.vsi_l_offset buff_len = 0
+
+ filename_b = filename if not isinstance(filename, string_types) else filename.encode('utf-8')
+ cfilename = filename_b
+ buff = _gdal.VSIGetMemFileBuffer(cfilename, &buff_len, 0)
+ n = buff_len
+ log.debug("Buffer length: %d bytes", n)
+ cdef np.uint8_t[:] buff_view = <np.uint8_t[:n]>buff
+ return buff_view
diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx
index b84470a..8e13061 100644
--- a/rasterio/_warp.pyx
+++ b/rasterio/_warp.pyx
@@ -156,11 +156,16 @@ def _transform_geom(
def _reproject(
source, destination,
- src_transform=None, src_crs=None,
- dst_transform=None, dst_crs=None,
- resampling=RESAMPLING.nearest,
+ src_transform=None,
+ src_crs=None,
+ src_nodata=None,
+ dst_transform=None,
+ dst_crs=None,
+ dst_nodata=None,
+ resampling=RESAMPLING.nearest,
**kwargs):
- """Reproject a source raster to a destination.
+ """
+ Reproject a source raster to a destination raster.
If the source and destination are ndarrays, coordinate reference
system definitions and affine transformation parameters are required
@@ -169,25 +174,73 @@ def _reproject(
If the source and destination are rasterio Bands, shorthand for
bands of datasets on disk, the coordinate reference systems and
transforms will be read from the appropriate datasets.
+
+ Parameters
+ ------------
+ source: ndarray or rasterio Band
+ Source raster.
+ destination: ndarray or rasterio Band
+ Target raster.
+ src_transform: affine transform object, optional
+ Source affine transformation. Required if source and destination
+ are ndarrays. Will be derived from source if it is a rasterio Band.
+ src_crs: dict, optional
+ Source coordinate reference system, in rasterio dict format.
+ Required if source and destination are ndarrays.
+ Will be derived from source if it is a rasterio Band.
+ Example: {'init': 'EPSG:4326'}
+ src_nodata: int or float, optional
+ The source nodata value. Pixels with this value will not be used
+ for interpolation. If not set, it will be default to the
+ nodata value of the source image if a masked ndarray or rasterio band,
+ if available. Must be provided if dst_nodata is not None.
+ dst_transform: affine transform object, optional
+ Target affine transformation. Required if source and destination
+ are ndarrays. Will be derived from target if it is a rasterio Band.
+ dst_crs: dict, optional
+ Target coordinate reference system. Required if source and destination
+ are ndarrays. Will be derived from target if it is a rasterio Band.
+ dst_nodata: int or float, optional
+ The nodata value used to initialize the destination; it will remain
+ in all areas not covered by the reprojected source. Defaults to the
+ nodata value of the destination image (if set), the value of
+ src_nodata, or 0 (gdal default).
+ resampling: int
+ Resampling method to use. One of the following:
+ RESAMPLING.nearest,
+ RESAMPLING.bilinear,
+ RESAMPLING.cubic,
+ RESAMPLING.cubic_spline,
+ RESAMPLING.lanczos,
+ RESAMPLING.average,
+ RESAMPLING.mode
+ kwargs: dict, optional
+ Additional arguments passed to transformation function.
+
+ Returns
+ ---------
+ out: None
+ Output is written to destination.
"""
- cdef int retval=0, rows, cols
- cdef void *hrdriver
- cdef void *hdsin
- cdef void *hdsout
- cdef void *hbandin
- cdef void *hbandout
+
+ cdef int retval=0, rows, cols, src_count
+ cdef void *hrdriver = NULL
+ cdef void *hdsin = NULL
+ cdef void *hdsout = NULL
+ cdef void *hbandin = NULL
+ cdef void *hbandout = NULL
cdef _io.RasterReader rdr
cdef _io.RasterUpdater udr
cdef _io.GDALAccess GA
cdef double gt[6]
cdef char *srcwkt = NULL
cdef char *dstwkt= NULL
- cdef const char *proj_c
- cdef void *osr
- cdef char **warp_extras
- cdef char *key_c
- cdef char *val_c
- cdef const char* pszWarpThreads
+ cdef const char *proj_c = NULL
+ cdef void *osr = NULL
+ cdef char **warp_extras = NULL
+ cdef char *key_c = NULL
+ cdef char *val_c = NULL
+ cdef const char* pszWarpThread = NULL
# If the source is an ndarray, we copy to a MEM dataset.
# We need a src_transform and src_dst in this case. These will
@@ -200,6 +253,10 @@ def _reproject(
rows = source.shape[1]
cols = source.shape[2]
dtype = np.dtype(source.dtype).name
+ if src_nodata is None and hasattr(source, 'fill_value'):
+ # source is a masked array
+ src_nodata = source.fill_value
+
hrdriver = _gdal.GDALGetDriverByName("MEM")
if hrdriver == NULL:
raise ValueError("NULL driver for 'MEM'")
@@ -216,19 +273,7 @@ def _reproject(
gt[i] = src_transform[i]
retval = _gdal.GDALSetGeoTransform(hdsin, gt)
log.debug("Set transform on temp source dataset: %d", retval)
- osr = _gdal.OSRNewSpatialReference(NULL)
- if osr == NULL:
- raise ValueError("Null spatial reference")
- params = []
- for k, v in src_crs.items():
- if v is True or (k == 'no_defs' and v):
- params.append("+%s" % k)
- else:
- params.append("+%s=%s" % (k, v))
- proj = " ".join(params)
- proj_b = proj.encode()
- proj_c = proj_b
- _gdal.OSRImportFromProj4(osr, proj_c)
+ osr = _base._osr_from_crs(src_crs)
_gdal.OSRExportToWkt(osr, &srcwkt)
_gdal.GDALSetProjection(hdsin, srcwkt)
_gdal.CPLFree(srcwkt)
@@ -236,34 +281,7 @@ def _reproject(
log.debug("Set CRS on temp source dataset: %s", srcwkt)
# Copy arrays to the dataset.
- #hbandin = _gdal.GDALGetRasterBand(hdsin, 1)
- #if hbandin == NULL:
- # raise ValueError("NULL input band")
- #log.debug("Got temp source band")
- indexes = np.array(range(1, src_count+1))
- if dtype == dtypes.ubyte:
- retval = _io.io_multi_ubyte(
- hdsin, 1, 0, 0, cols, rows, source, indexes, src_count)
- elif dtype == dtypes.uint16:
- retval = _io.io_multi_uint16(
- hdsin, 1, 0, 0, cols, rows, source, indexes, src_count)
- elif dtype == dtypes.int16:
- retval = _io.io_multi_int16(
- hdsin, 1, 0, 0, cols, rows, source, indexes, src_count)
- elif dtype == dtypes.uint32:
- retval = _io.io_multi_uint32(
- hdsin, 1, 0, 0, cols, rows, source, indexes, src_count)
- elif dtype == dtypes.int32:
- retval = _io.io_multi_int32(
- hdsin, 1, 0, 0, cols, rows, source, indexes, src_count)
- elif dtype == dtypes.float32:
- retval = _io.io_multi_float32(
- hdsin, 1, 0, 0, cols, rows, source, indexes, src_count)
- elif dtype == dtypes.float64:
- retval = _io.io_multi_float64(
- hdsin, 1, 0, 0, cols, rows, source, indexes, src_count)
- else:
- raise ValueError("Invalid dtype")
+ retval = _io.io_auto(source, hdsin, 1)
# TODO: handle errors (by retval).
log.debug("Wrote array to temp source dataset")
@@ -272,6 +290,8 @@ def _reproject(
rdr = source.ds
hdsin = rdr._hds
src_count = 1
+ if src_nodata is None:
+ src_nodata = rdr.nodata
else:
raise ValueError("Invalid source")
@@ -297,30 +317,22 @@ def _reproject(
gt[i] = dst_transform[i]
retval = _gdal.GDALSetGeoTransform(hdsout, gt)
log.debug("Set transform on temp destination dataset: %d", retval)
- osr = _gdal.OSRNewSpatialReference(NULL)
- if osr == NULL:
- raise ValueError("Null spatial reference")
- params = []
- for k, v in dst_crs.items():
- if v is True or (k == 'no_defs' and v):
- params.append("+%s" % k)
- else:
- params.append("+%s=%s" % (k, v))
- proj = " ".join(params)
- log.debug("Proj4 string: %s", proj)
- proj_b = proj.encode()
- proj_c = proj_b
- _gdal.OSRImportFromProj4(osr, proj_c)
+ osr = _base._osr_from_crs(dst_crs)
_gdal.OSRExportToWkt(osr, &dstwkt)
retval = _gdal.GDALSetProjection(hdsout, dstwkt)
log.debug("Setting Projection: %d", retval)
_gdal.CPLFree(dstwkt)
_gdal.OSRDestroySpatialReference(osr)
log.debug("Set CRS on temp destination dataset: %s", dstwkt)
+ if dst_nodata is None and hasattr(destination, "fill_value"):
+ # destination is a masked array
+ dst_nodata = destination.fill_value
elif isinstance(destination, tuple):
udr = destination.ds
hdsout = udr._hds
+ if dst_nodata is None:
+ dst_nodata = udr.nodata
else:
raise ValueError("Invalid destination")
@@ -338,8 +350,12 @@ def _reproject(
log.debug("Created transformer")
psWOptions = _gdal.GDALCreateWarpOptions()
-
+
+ # Note: warp_extras is pointed to different memory locations on every
+ # call to CSLSetNameValue call below, but needs to be set here to
+ # get the defaults
warp_extras = psWOptions.papszWarpOptions
+
for k, v in kwargs.items():
k, v = k.upper(), str(v).upper()
key_b = k.encode('utf-8')
@@ -347,6 +363,7 @@ def _reproject(
key_c = key_b
val_c = val_b
warp_extras = _gdal.CSLSetNameValue(warp_extras, key_c, val_c)
+ log.debug("Setting warp option %s: %s" % (k, v))
pszWarpThreads = _gdal.CSLFetchNameValue(warp_extras, "NUM_THREADS")
if pszWarpThreads == NULL:
@@ -358,6 +375,54 @@ def _reproject(
log.debug("Created warp options")
psWOptions.eResampleAlg = <_gdal.GDALResampleAlg>resampling
+
+ # Set src_nodata and dst_nodata
+ if src_nodata is None and dst_nodata is not None:
+ raise ValueError("src_nodata must be provided because dst_nodata "
+ "is not None")
+ log.debug("src_nodata: %s" % src_nodata)
+
+ if dst_nodata is None:
+ if src_nodata is not None:
+ dst_nodata = src_nodata
+ else:
+ dst_nodata = 0 # GDAL default
+ log.debug("dst_nodata: %s" % dst_nodata)
+
+ # Validate nodata values
+ if src_nodata is not None:
+ if not _io.in_dtype_range(src_nodata, source.dtype):
+ raise ValueError("src_nodata must be in valid range for "
+ "source dtype")
+
+ psWOptions.padfSrcNoDataReal = <double*>_gdal.CPLMalloc(
+ src_count * sizeof(double))
+ psWOptions.padfSrcNoDataImag = <double*>_gdal.CPLMalloc(
+ src_count * sizeof(double))
+ for i in range(src_count):
+ psWOptions.padfSrcNoDataReal[i] = src_nodata
+ psWOptions.padfSrcNoDataImag[i] = 0.0
+ warp_extras = _gdal.CSLSetNameValue(
+ warp_extras, "UNIFIED_SRC_NODATA", "YES")
+
+
+ if dst_nodata is not None and not _io.in_dtype_range(
+ dst_nodata, destination.dtype):
+ raise ValueError("dst_nodata must be in valid range for "
+ "destination dtype")
+
+ psWOptions.padfDstNoDataReal = <double*>_gdal.CPLMalloc(src_count * sizeof(double))
+ psWOptions.padfDstNoDataImag = <double*>_gdal.CPLMalloc(src_count * sizeof(double))
+ for i in range(src_count):
+ psWOptions.padfDstNoDataReal[i] = dst_nodata
+ psWOptions.padfDstNoDataImag[i] = 0.0
+ warp_extras = _gdal.CSLSetNameValue(
+ warp_extras, "INIT_DEST", "NO_DATA")
+
+ # Important: set back into struct or values set above are lost
+ # This is because CSLSetNameValue returns a new list each time
+ psWOptions.papszWarpOptions = warp_extras
+
# TODO: Approximate transformations.
#if maxerror > 0.0:
# psWOptions.pTransformerArg = _gdal.GDALCreateApproxTransformer(
@@ -385,11 +450,14 @@ def _reproject(
psWOptions.panDstBands[i] = i+1
log.debug("Set transformer options")
- # TODO: Src nodata and alpha band.
+ # TODO: alpha band.
eErr = oWarper.Initialize(psWOptions)
if eErr == 0:
- _, rows, cols = destination.shape
+
+ log.debug("Destination shape: %r", destination.shape)
+
+ rows, cols = destination.shape[-2:]
log.debug(
"Chunk and warp window: %d, %d, %d, %d",
0, 0, cols, rows)
@@ -400,6 +468,7 @@ def _reproject(
except Exception:
log.exception(
"Caught exception in warping. Source not reprojected.")
+ raise
else:
reprojected = True
@@ -416,44 +485,9 @@ def _reproject(
_gdal.GDALClose(hdsin)
if reprojected and isinstance(destination, np.ndarray):
- try:
- dtype = np.dtype(destination.dtype).name
- _, rows, cols = destination.shape
- indexes = np.array(range(1, src_count+1))
- if dtype == dtypes.ubyte:
- retval = _io.io_multi_ubyte(
- hdsout, 0, 0, 0, cols, rows,
- destination, indexes, src_count)
- elif dtype == dtypes.uint16:
- retval = _io.io_multi_uint16(
- hdsout, 0, 0, 0, cols, rows,
- destination, indexes, src_count)
- elif dtype == dtypes.int16:
- retval = _io.io_multi_int16(
- hdsout, 0, 0, 0, cols, rows,
- destination, indexes, src_count)
- elif dtype == dtypes.uint32:
- retval = _io.io_multi_uint32(
- hdsout, 0, 0, 0, cols, rows,
- destination, indexes, src_count)
- elif dtype == dtypes.int32:
- retval = _io.io_multi_int32(
- hdsout, 0, 0, 0, cols, rows,
- destination, indexes, src_count)
- elif dtype == dtypes.float32:
- retval = _io.io_multi_float32(
- hdsout, 0, 0, 0, cols, rows,
- destination, indexes, src_count)
- elif dtype == dtypes.float64:
- retval = _io.io_multi_float64(
- hdsout, 0, 0, 0, cols, rows,
- destination, indexes, src_count)
- else:
- raise ValueError("Invalid dtype")
- # TODO: handle errors (by retval).
- except Exception:
- raise
- finally:
- if hdsout != NULL:
- _gdal.GDALClose(hdsout)
+ retval = _io.io_auto(destination, hdsout, 0)
+ # TODO: handle errors (by retval).
+
+ if hdsout != NULL:
+ _gdal.GDALClose(hdsout)
diff --git a/rasterio/crs.py b/rasterio/crs.py
index eef744c..c70f03e 100644
--- a/rasterio/crs.py
+++ b/rasterio/crs.py
@@ -10,6 +10,7 @@
# {'proj': 'longlat', 'ellps': 'WGS84', 'datum': 'WGS84', 'no_defs': True}
#
+from rasterio._base import is_geographic_crs, is_projected_crs
from rasterio.five import string_types
def to_string(crs):
diff --git a/rasterio/features.py b/rasterio/features.py
index 3673c3a..d52da1a 100644
--- a/rasterio/features.py
+++ b/rasterio/features.py
@@ -22,6 +22,49 @@ class NullHandler(logging.Handler):
log.addHandler(NullHandler())
+def geometry_mask(
+ geometries,
+ out_shape,
+ transform,
+ all_touched=False,
+ invert=False):
+ """Create a mask from shapes. By default, mask is intended for use as a
+ numpy mask, where pixels that overlap shapes are False.
+
+ Parameters
+ ----------
+ geometries : iterable over geometries (GeoJSON-like objects)
+ out_shape : tuple or list
+ Shape of output numpy ndarray.
+ transform : Affine transformation object
+ Transformation from pixel coordinates of `image` to the
+ coordinate system of the input `shapes`. See the `transform`
+ property of dataset objects.
+ all_touched : boolean, optional
+ If True, all pixels touched by geometries will be burned in. If
+ false, only pixels whose center is within the polygon or that
+ are selected by Bresenham's line algorithm will be burned in.
+ invert: boolean, optional
+ If True, mask will be True for pixels that overlap shapes.
+ False by default.
+
+ Returns
+ -------
+ out : numpy ndarray of type 'bool'
+ Result
+ """
+
+ fill, mask_value = (0, 1) if invert else (1, 0)
+
+ return rasterize(
+ geometries,
+ out_shape=out_shape,
+ transform=transform,
+ all_touched=all_touched,
+ fill=fill,
+ default_value=mask_value).astype('bool')
+
+
def shapes(image, mask=None, connectivity=4, transform=IDENTITY):
"""
Return a generator of (polygon, value) for each each set of adjacent pixels
@@ -34,8 +77,8 @@ def shapes(image, mask=None, connectivity=4, transform=IDENTITY):
Data type must be one of rasterio.int16, rasterio.int32,
rasterio.uint8, rasterio.uint16, or rasterio.float32.
mask : numpy ndarray or rasterio Band object, optional
- Values of False will be excluded from feature generation
- Must be of type rasterio.bool_
+ Values of False or 0 will be excluded from feature generation
+ Must evaluate to bool (rasterio.bool_ or rasterio.uint8)
connectivity : int, optional
Use 4 or 8 pixel connectivity for grouping pixels into features
transform : Affine transformation, optional
@@ -67,8 +110,8 @@ def shapes(image, mask=None, connectivity=4, transform=IDENTITY):
raise ValueError('image dtype must be one of: %s'
% (', '.join(valid_dtypes)))
- if mask is not None and np.dtype(mask.dtype) != np.dtype(rasterio.bool_):
- raise ValueError("Mask must be dtype rasterio.bool_")
+ if mask is not None and np.dtype(mask.dtype).name not in ('bool', 'uint8'):
+ raise ValueError("Mask must be dtype rasterio.bool_ or rasterio.uint8")
if connectivity not in (4, 8):
raise ValueError("Connectivity Option must be 4 or 8")
@@ -99,8 +142,8 @@ def sieve(image, size, out=None, output=None, mask=None, connectivity=4):
output : older alias for `out`, will be removed before 1.0.
output : numpy ndarray, optional
mask : numpy ndarray or rasterio Band object, optional
- Values of False will be excluded from feature generation
- Must be of type rasterio.bool_
+ Values of False or 0 will be excluded from feature generation
+ Must evaluate to bool (rasterio.bool_ or rasterio.uint8)
connectivity : int, optional
Use 4 or 8 pixel connectivity for grouping pixels into features
@@ -127,7 +170,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')
@@ -140,8 +183,9 @@ def sieve(image, size, out=None, output=None, mask=None, connectivity=4):
raise ValueError('connectivity must be 4 or 8')
if mask is not None:
- if np.dtype(mask.dtype) != np.dtype(rasterio.bool_):
- raise ValueError('Mask must be dtype rasterio.bool_')
+ if np.dtype(mask.dtype) not in ('bool', 'uint8'):
+ raise ValueError('Mask must be dtype rasterio.bool_ or '
+ 'rasterio.uint8')
elif mask.shape != image.shape:
raise ValueError('mask shape must be same as image shape')
@@ -155,12 +199,15 @@ 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')
+ raise ValueError('out raster must match dtype of image')
elif out.shape != image.shape:
- raise ValueError('mask shape must be same as image shape')
+ raise ValueError('out raster shape must be same as image shape')
with rasterio.drivers():
_sieve(image, size, out, mask, connectivity)
@@ -183,7 +230,7 @@ def rasterize(
Parameters
----------
shapes : iterable of (geometry, value) pairs or iterable over
- geometries `geometry` can either be an object that implements
+ geometries. `geometry` can either be an object that implements
the geo interface or GeoJSON-like object.
out_shape : tuple or list
Shape of output numpy ndarray.
@@ -201,11 +248,11 @@ def rasterize(
all_touched : boolean, optional
If True, all pixels touched by geometries will be burned in. If
false, only pixels whose center is within the polygon or that
- are selected by brezenhams line algorithm will be burned in.
+ are selected by Bresenham's line algorithm will be burned in.
default_value : int or float, optional
Used as value for all geometries, if not provided in `shapes`.
dtype : rasterio or numpy data type, optional
- Used as data type for results, if `output` is not provided.
+ Used as data type for results, if `out` is not provided.
Returns
-------
diff --git a/rasterio/fill.py b/rasterio/fill.py
index 37ec58e..605add7 100644
--- a/rasterio/fill.py
+++ b/rasterio/fill.py
@@ -1,41 +1,46 @@
import rasterio
from rasterio._fill import _fillnodata
+
def fillnodata(
image,
mask=None,
max_search_distance=100.0,
smoothing_iterations=0):
"""
- Fill selected raster regions by interpolation from the edges.
+ Fill holes in a raster dataset by interpolation from the edges.
+
+ This algorithm will interpolate values for all designated nodata
+ pixels (marked by zeros in `mask`). For each pixel a four direction
+ conic search is done to find values to interpolate from (using
+ inverse distance weighting). Once all values are interpolated, zero
+ or more smoothing iterations (3x3 average filters on interpolated
+ pixels) are applied to smooth out artifacts.
+
+ This algorithm is generally suitable for interpolating missing
+ regions of fairly continuously varying rasters (such as elevation
+ models for instance). It is also suitable for filling small holes
+ and cracks in more irregularly varying images (like aerial photos).
+ It is generally not so great for interpolating a raster from sparse
+ point data.
- This algorithm will interpolate values for all designated nodata pixels
- (marked by zeros in `mask`). For each pixel a four direction conic search
- is done to find values to interpolate from (using inverse distance
- weighting). Once all values are interpolated, zero or more smoothing
- iterations (3x3 average filters on interpolated pixels) are applied to
- smooth out artifacts.
-
- This algorithm is generally suitable for interpolating missing regions of
- fairly continuously varying rasters (such as elevation models for
- instance). It is also suitable for filling small holes and cracks in more
- irregularly varying images (like aerial photos). It is generally not so
- great for interpolating a raster from sparse point data.
-
Parameters
----------
image : numpy ndarray
- The source containing nodata holes.
- mask : numpy ndarray
+ The source containing nodata holes.
+ mask : numpy ndarray or None
A mask band indicating which pixels to interpolate. Pixels to
- interpolate into are indicated by the value 0. Values of 1 indicate
- areas to use during interpolation. Must be same shape as image.
+ interpolate into are indicated by the value 0. Values > 0
+ indicate areas to use during interpolation. Must be same shape
+ as image. If `None`, a mask will be diagnosed from the source
+ data.
max_search_distance : float, optional
- The maxmimum number of pixels to search in all directions to find
- values to interpolate from. The default is 100.
+ The maxmimum number of pixels to search in all directions to
+ find values to interpolate from. The default is 100.
smoothing_iterations : integer, optional
- The number of 3x3 smoothing filter passes to run. The default is 0.
-
+ The number of 3x3 smoothing filter passes to run. The default is
+ 0.
+
Returns
-------
out : numpy ndarray
@@ -44,4 +49,5 @@ def fillnodata(
max_search_distance = float(max_search_distance)
smoothing_iterations = int(smoothing_iterations)
with rasterio.drivers():
- return _fillnodata(image, mask, max_search_distance, smoothing_iterations)
+ return _fillnodata(
+ image, mask, max_search_distance, smoothing_iterations)
diff --git a/rasterio/profiles.py b/rasterio/profiles.py
new file mode 100644
index 0000000..3b2bedf
--- /dev/null
+++ b/rasterio/profiles.py
@@ -0,0 +1,48 @@
+"""Raster dataset profiles."""
+
+from rasterio.dtypes import uint8
+
+
+class Profile:
+ """Base class for Rasterio dataset profiles.
+
+ Subclasses will declare a format driver and driver-specific
+ creation options.
+ """
+ driver = None
+ defaults = {}
+
+ def __call__(self, **kwargs):
+ """Returns a mapping of keyword args for writing a new datasets.
+
+ Example:
+
+ profile = SomeProfile()
+ with rasterio.open('foo.tif', 'w', **profile()) as dst:
+ # Write data ...
+
+ """
+ if kwargs.get('driver', self.driver) != self.driver:
+ raise ValueError(
+ "Overriding this profile's driver is not allowed.")
+ profile = self.defaults.copy()
+ profile.update(**kwargs)
+ profile['driver'] = self.driver
+ return profile
+
+
+class DefaultGTiffProfile(Profile):
+ """A tiled, band-interleaved, LZW-compressed, 8-bit GTiff profile."""
+ driver = 'GTiff'
+ defaults = {
+ 'interleave': 'band',
+ 'tiled': True,
+ 'blockxsize': 256,
+ 'blockysize': 256,
+ 'compress': 'lzw',
+ 'nodata': 0,
+ 'dtype': uint8
+ }
+
+
+default_gtiff_profile = DefaultGTiffProfile()
diff --git a/rasterio/rio/bands.py b/rasterio/rio/bands.py
index 0350e27..b167f44 100644
--- a/rasterio/rio/bands.py
+++ b/rasterio/rio/bands.py
@@ -1,5 +1,4 @@
import logging
-import os.path
import sys
import click
@@ -8,7 +7,7 @@ from cligj import files_inout_arg, format_opt
import rasterio
from rasterio.five import zip_longest
-from rasterio.rio.cli import cli
+from rasterio.rio.cli import cli, bidx_mult_opt, output_opt, resolve_inout
PHOTOMETRIC_CHOICES = [val.lower() for val in [
@@ -25,14 +24,14 @@ PHOTOMETRIC_CHOICES = [val.lower() for val in [
# Stack command.
@cli.command(short_help="Stack a number of bands into a multiband dataset.")
@files_inout_arg
+ at output_opt
@format_opt
- at click.option('--bidx', multiple=True,
- help="Indexes of input file bands.")
+ at bidx_mult_opt
@click.option('--photometric', default=None,
type=click.Choice(PHOTOMETRIC_CHOICES),
help="Photometric interpretation")
@click.pass_context
-def stack(ctx, files, driver, bidx, photometric):
+def stack(ctx, files, output, driver, bidx, photometric):
"""Stack a number of bands from one or more input files into a
multiband dataset.
@@ -70,8 +69,7 @@ def stack(ctx, files, driver, bidx, photometric):
logger = logging.getLogger('rio')
try:
with rasterio.drivers(CPL_DEBUG=verbosity>2):
- output = files[-1]
- files = files[:-1]
+ output, files = resolve_inout(files=files, output=output)
output_count = 0
indexes = []
for path, item in zip_longest(files, bidx, fillvalue=None):
@@ -121,7 +119,6 @@ def stack(ctx, files, driver, bidx, photometric):
dst.write(data, range(dst_idx, dst_idx+len(index)))
dst_idx += len(index)
- sys.exit(0)
except Exception:
- logger.exception("Failed. Exception caught")
- sys.exit(1)
+ logger.exception("Exception caught during processing")
+ raise click.Abort()
diff --git a/rasterio/rio/calc.py b/rasterio/rio/calc.py
new file mode 100644
index 0000000..9a8c4bb
--- /dev/null
+++ b/rasterio/rio/calc.py
@@ -0,0 +1,146 @@
+# Calc command.
+
+import logging
+import sys
+import traceback
+
+import click
+import snuggs
+from cligj import files_inout_arg
+
+import rasterio
+from rasterio.fill import fillnodata
+from rasterio.features import sieve
+from rasterio.rio.cli import (
+ cli, dtype_opt, masked_opt, output_opt, resolve_inout)
+
+
+def get_bands(inputs, d, i=None):
+ """Get a rasterio.Band object from calc's inputs"""
+ path = inputs[d] if d in dict(inputs) else inputs[int(d)-1][1]
+ src = rasterio.open(path)
+ return (rasterio.band(src, i) if i else
+ [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
+
+
+ at cli.command(short_help="Raster data calculator.")
+ at click.argument('command')
+ at files_inout_arg
+ at output_opt
+ at click.option('--name', multiple=True,
+ help='Specify an input file with a unique short (alphas only) '
+ 'name for use in commands like '
+ '"a=tests/data/RGB.byte.tif".')
+ at dtype_opt
+ at masked_opt
+ at click.pass_context
+def calc(ctx, command, files, output, name, dtype, masked):
+ """A raster data calculator
+
+ Evaluates an expression using input datasets and writes the result
+ to a new dataset.
+
+ Command syntax is lisp-like. An expression consists of an operator
+ or function name and one or more strings, numbers, or expressions
+ enclosed in parentheses. Functions include ``read`` (gets a raster
+ array) and ``asarray`` (makes a 3-D array from 2-D arrays).
+
+ \b
+ * (read i) evaluates to the i-th input dataset (a 3-D array).
+ * (read i j) evaluates to the j-th band of the i-th dataset (a 2-D
+ array).
+ * (take foo j) evaluates to the j-th band of a dataset named foo (see
+ help on the --name option above).
+ * Standard numpy array operators (+, -, *, /) are available.
+ * When the final result is a list of arrays, a multi band output
+ file is written.
+ * When the final result is a single array, a single band output
+ file is written.
+
+ Example:
+
+ \b
+ $ rio calc "(+ 2 (* 0.95 (read 1)))" tests/data/RGB.byte.tif \\
+ > /tmp/out.tif
+
+ Produces a 3-band GeoTIFF with all values scaled by 0.95 and
+ incremented by 2.
+
+ \b
+ $ rio calc "(asarray (+ 125 (read 1)) (read 1) (read 1))" \\
+ > tests/data/shade.tif /tmp/out.tif
+
+ Produces a 3-band RGB GeoTIFF, with red levels incremented by 125,
+ from the single-band input.
+
+ """
+ import numpy as np
+
+ verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+ logger = logging.getLogger('rio')
+
+ try:
+ with rasterio.drivers(CPL_DEBUG=verbosity > 2):
+ output, files = resolve_inout(files=files, output=output)
+
+ inputs = ([tuple(n.split('=')) for n in name] +
+ [(None, n) for n in files])
+
+ with rasterio.open(inputs[0][1]) as first:
+ kwargs = first.meta
+ kwargs['transform'] = kwargs.pop('affine')
+ dtype = dtype or first.meta['dtype']
+ kwargs['dtype'] = dtype
+
+ ctxkwds = {}
+ for i, (name, path) in enumerate(inputs):
+ with rasterio.open(path) as src:
+ # Using the class method instead of instance
+ # method. Latter raises
+ #
+ # TypeError: astype() got an unexpected keyword
+ # argument 'copy'
+ #
+ # possibly something to do with the instance being
+ # a masked array.
+ ctxkwds[name or '_i%d' % (i+1)] = src.read(masked=masked)
+
+ # 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)
+
+ if len(res.shape) == 3:
+ results = np.ndarray.astype(res, dtype, copy=False)
+ else:
+ results = np.asanyarray(
+ [np.ndarray.astype(res, dtype, copy=False)])
+
+ kwargs['count'] = results.shape[0]
+
+ with rasterio.open(output, 'w', **kwargs) as dst:
+ dst.write(results)
+
+ except snuggs.ExpressionError as err:
+ click.echo("Expression Error:")
+ click.echo(' %s' % err.text)
+ click.echo(' ' + ' ' * err.offset + "^")
+ click.echo(err)
+ raise click.Abort()
+ #except Exception as err:
+ # t, v, tb = sys.exc_info()
+ # for line in traceback.format_exception_only(t, v):
+ # click.echo(line, nl=False)
+ # raise click.Abort()
diff --git a/rasterio/rio/cli.py b/rasterio/rio/cli.py
index 37f2ff2..7bfea62 100644
--- a/rasterio/rio/cli.py
+++ b/rasterio/rio/cli.py
@@ -1,6 +1,10 @@
+"""Rasterio's command line interface core."""
+
import json
import logging
+import os
import sys
+import traceback
import click
from cligj import verbose_opt, quiet_opt
@@ -13,8 +17,52 @@ def configure_logging(verbosity):
logging.basicConfig(stream=sys.stderr, level=log_level)
+class BrokenCommand(click.Command):
+ """A dummy command that provides help for broken plugins."""
+
+ def __init__(self, name):
+ click.Command.__init__(self, name)
+ self.help = (
+ "Warning: entry point could not be loaded. Contact "
+ "its author for help.\n\n\b\n"
+ + traceback.format_exc())
+ self.short_help = (
+ "Warning: could not load plugin. See `rio %s --help`." % self.name)
+
+
+class RioGroup(click.Group):
+ """Custom formatting for the commands of broken plugins."""
+
+ def format_commands(self, ctx, formatter):
+ """Extra format methods for multi methods that adds all the commands
+ after the options.
+ """
+ rows = []
+ for subcommand in self.list_commands(ctx):
+ cmd = self.get_command(ctx, subcommand)
+ # What is this, the tool lied about a command. Ignore it
+ if cmd is None:
+ continue
+
+ help = cmd.short_help or ''
+
+ # Mark broken subcommands with a pile of poop.
+ name = cmd.name
+ if isinstance(cmd, BrokenCommand):
+ if os.environ.get('RIO_HONESTLY'):
+ name += u'\U0001F4A9'
+ else:
+ name += u'\u2020'
+
+ rows.append((name, help))
+
+ if rows:
+ with formatter.section('Commands'):
+ formatter.write_dl(rows)
+
+
# The CLI command group.
- at click.group(help="Rasterio command line interface.")
+ at click.group(help="Rasterio command line interface.", cls=RioGroup)
@verbose_opt
@quiet_opt
@click.version_option(version=rasterio.__version__, message='%(version)s')
@@ -26,6 +74,118 @@ def cli(ctx, verbose, quiet):
ctx.obj['verbosity'] = verbosity
+# Common arguments and options
+
+# TODO: move file_in_arg and file_out_arg to cligj
+
+# Singular input file
+file_in_arg = click.argument(
+ 'INPUT',
+ type=click.Path(exists=True, resolve_path=True))
+
+# Singular output file
+file_out_arg = click.argument(
+ 'OUTPUT',
+ type=click.Path(resolve_path=True))
+
+bidx_opt = click.option(
+ '-b', '--bidx',
+ type=int,
+ default=1,
+ help="Input file band index (default: 1)")
+
+bidx_mult_opt = click.option(
+ '-b', '--bidx',
+ multiple=True,
+ help="Indexes of input file bands.")
+
+# TODO: may be better suited to cligj
+bounds_opt = click.option(
+ '--bounds',
+ nargs=4, type=float, default=None,
+ help='Output bounds: left, bottom, right, top.')
+
+dtype_opt = click.option(
+ '-t', '--dtype',
+ type=click.Choice([
+ 'ubyte', 'uint8', 'uint16', 'int16', 'uint32', 'int32',
+ 'float32', 'float64']),
+ default=None,
+ help="Output data type (default: float64).")
+
+like_file_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.')
+
+masked_opt = click.option(
+ '--masked/--not-masked',
+ default=True,
+ help="Evaluate expressions using masked arrays (the default) or ordinary "
+ "numpy arrays.")
+
+output_opt = click.option(
+ '-o', '--output',
+ default=None,
+ type=click.Path(resolve_path=True),
+ help="Path to output file (optional alternative to a positional arg "
+ "for some commands).")
+
+
+resolution_opt = click.option(
+ '-r', '--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')
+
+"""
+Registry of command line options (also see cligj options):
+-a, --all: Use all pixels touched by features. In rio-mask, rio-rasterize
+--as-mask/--not-as-mask: interpret band as mask or not. In rio-shapes
+--band/--mask: use band or mask. In rio-shapes
+--bbox:
+-b, --bidx: band index(es) (singular or multiple value versions).
+ In rio-info, rio-sample, rio-shapes, rio-stack (different usages)
+--bounds: bounds in world coordinates.
+ In rio-info, rio-rasterize (different usages)
+--count: count of bands. In rio-info
+--crop: Crop raster to extent of features. In rio-mask
+--crs: CRS of input raster. In rio-info
+--default-value: default for rasterized pixels. In rio-rasterize
+--dimensions: Output width, height. In rio-rasterize
+--dst-crs: destination CRS. In rio-transform
+--fill: fill value for pixels not covered by features. In rio-rasterize
+--formats: list available formats. In rio-info
+--height: height of raster. In rio-info
+-i, --invert: Invert mask created from features: In rio-mask
+-j, --geojson-mask: GeoJSON for masking raster. In rio-mask
+--lnglat: geograhpic coordinates of center of raster. In rio-info
+--masked/--not-masked: read masked data from source file.
+ In rio-calc, rio-info
+-m, --mode: output file mode (r, r+). In rio-insp
+--name: input file name alias. In rio-calc
+--nodata: nodata value. In rio-info, rio-merge (different usages)
+--photometric: photometric interpretation. In rio-stack
+--property: GeoJSON property to use as values for rasterize. In rio-rasterize
+-r, --res: output resolution.
+ In rio-info, rio-rasterize (different usages. TODO: try to combine
+ usages, prefer rio-rasterize version)
+--sampling: Inverse of sampling fraction. In rio-shapes
+--shape: shape (width, height) of band. In rio-info
+--src-crs: source CRS.
+ In rio-insp, rio-rasterize (different usages. TODO: consolidate usages)
+--stats: print raster stats. In rio-inf
+-t, --dtype: data type. In rio-calc, rio-info (different usages)
+--width: width of raster. In rio-info
+--with-nodata/--without-nodata: include nodata regions or not. In rio-shapes.
+-v, --tell-me-more, --verbose
+"""
+
+
def coords(obj):
"""Yield all coordinate coordinate tuples from a geometry or feature.
From python-geojson package."""
@@ -67,7 +227,7 @@ def write_features(
'bbox': bbox,
'features': [feat]}, **dump_kwds))
fobj.write('\n')
- # Aggregate all features into a single object expressed as
+ # Aggregate all features into a single object expressed as
# bbox or collection.
else:
features = list(collection())
@@ -78,7 +238,18 @@ def write_features(
else:
fobj.write(json.dumps({
'bbox': collection.bbox,
- 'type': 'FeatureCollection',
+ 'type': 'FeatureCollection',
'features': features},
**dump_kwds))
fobj.write('\n')
+
+
+def resolve_inout(input=None, output=None, files=None):
+ """Resolves inputs and outputs from standard args and options.
+
+ Returns `output_filename, [input_filename0, ...]`."""
+ resolved_output = output or (files[-1] if files else None)
+ resolved_inputs = (
+ [input] if input else [] +
+ list(files[:-1 if not output else None]) if files else [])
+ return resolved_output, resolved_inputs
diff --git a/rasterio/rio/features.py b/rasterio/rio/features.py
index ade9b31..2962096 100644
--- a/rasterio/rio/features.py
+++ b/rasterio/rio/features.py
@@ -1,10 +1,9 @@
-import functools
import json
import logging
-from math import ceil, floor
+from math import ceil
import os
import sys
-import warnings
+import shutil
import click
from cligj import (
@@ -15,16 +14,160 @@ from cligj import (
import rasterio
from rasterio.transform import Affine
-from rasterio.rio.cli import cli, coords, write_features
+from rasterio.rio.cli import (
+ cli, coords, write_features, file_in_arg, file_out_arg, like_file_opt,
+ bounds_opt, resolution_opt, output_opt, resolve_inout)
logger = logging.getLogger('rio')
-warnings.simplefilter('default')
+
+
+# Common options used below
+all_touched_opt = click.option(
+ '-a', '--all', '--all_touched', 'all_touched',
+ is_flag=True,
+ default=False,
+ help='Use all pixels touched by features, otherwise (default) use only '
+ 'pixels whose center is within the polygon or that are selected by '
+ 'Bresenhams line algorithm')
+
+
+# Mask command
+ at cli.command(short_help='Mask in raster using features.')
+ at files_inout_arg
+ at output_opt
+ at click.option('-j', '--geojson-mask', 'geojson_mask',
+ type=click.Path(), default=None,
+ help='GeoJSON file to use for masking raster. Use "-" to read '
+ 'from stdin. If not provided, original raster will be '
+ 'returned')
+ at format_opt
+ at all_touched_opt
+ at click.option('--crop', is_flag=True, default=False,
+ help='Crop output raster to the extent of the geometries. '
+ 'GeoJSON must overlap input raster to use --crop')
+ at click.option('-i', '--invert', is_flag=True, default=False,
+ help='Inverts the mask, so that areas covered by features are'
+ 'masked out and areas not covered are retained. Ignored '
+ 'if using --crop')
+ at click.pass_context
+def mask(
+ ctx,
+ files,
+ output,
+ geojson_mask,
+ driver,
+ all_touched,
+ crop,
+ invert):
+
+ """Masks in raster using GeoJSON features (masks out all areas not covered
+ by features), and optionally crops the output raster to the extent of the
+ features. Features are assumed to be in the same coordinate reference
+ system as the input raster.
+
+ GeoJSON must be the first input file or provided from stdin:
+
+ > rio mask input.tif output.tif --geojson-mask features.json
+
+ > rio mask input.tif output.tif --geojson-mask - < features.json
+
+ If the output raster exists, it will be completely overwritten with the
+ results of this operation.
+
+ The result is always equal to or within the bounds of the input raster.
+
+ --crop and --invert options are mutually exclusive.
+
+ --crop option is not valid if features are completely outside extent of
+ input raster.
+ """
+
+ from rasterio.features import geometry_mask
+ from rasterio.features import bounds as calculate_bounds
+
+ verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+
+ output, files = resolve_inout(files=files, output=output)
+ input = files[0]
+
+ if geojson_mask is None:
+ click.echo('No GeoJSON provided, INPUT will be copied to OUTPUT',
+ err=True)
+ shutil.copy(input, output)
+ return
+
+ if crop and invert:
+ click.echo('Invert option ignored when using --crop', err=True)
+ invert = False
+
+ with rasterio.drivers(CPL_DEBUG=verbosity > 2):
+ try:
+ geojson = json.loads(click.open_file(geojson_mask).read())
+ except ValueError:
+ raise click.BadParameter('GeoJSON could not be read from '
+ '--geojson-mask or stdin',
+ param_hint='--geojson-mask')
+
+ if 'features' in geojson:
+ geometries = (f['geometry'] for f in geojson['features'])
+ elif 'geometry' in geojson:
+ geometries = (geojson['geometry'], )
+ else:
+ raise click.BadParameter('Invalid GeoJSON', param=input,
+ param_hint='input')
+ bounds = geojson.get('bbox', calculate_bounds(geojson))
+
+ with rasterio.open(input) as src:
+ disjoint_bounds = _disjoint_bounds(bounds, src.bounds)
+
+ if crop:
+ if disjoint_bounds:
+ raise click.BadParameter('not allowed for GeoJSON outside '
+ 'the extent of the input raster',
+ param=crop, param_hint='--crop')
+
+ window = src.window(*bounds)
+ transform = src.window_transform(window)
+ (r1, r2), (c1, c2) = window
+ mask_shape = (r2 - r1, c2 - c1)
+ else:
+ if disjoint_bounds:
+ click.echo('GeoJSON outside bounds of existing output '
+ 'raster. Are they in different coordinate '
+ 'reference systems?',
+ err=True)
+
+ window = None
+ transform = src.affine
+ mask_shape = src.shape
+
+ mask = geometry_mask(
+ geometries,
+ out_shape=mask_shape,
+ transform=transform,
+ all_touched=all_touched,
+ invert=invert)
+
+ meta = src.meta.copy()
+ meta.update({
+ 'driver': driver,
+ 'height': mask.shape[0],
+ 'width': mask.shape[1],
+ 'transform': transform
+ })
+
+ with rasterio.open(output, 'w', **meta) as out:
+ for bidx in range(1, src.count + 1):
+ img = src.read(bidx, masked=True, window=window)
+ img.mask = img.mask | mask
+ out.write_band(bidx, img.filled(src.nodatavals[bidx-1]))
# Shapes command.
- at cli.command(short_help="Write the shapes of features.")
+ at cli.command(short_help="Write shapes extracted from bands or masks.")
@click.argument('input', type=click.Path(exists=True))
+ at output_opt
@precision_opt
@indent_opt
@compact_opt
@@ -34,21 +177,49 @@ warnings.simplefilter('default')
@use_rs_opt
@geojson_type_feature_opt(True)
@geojson_type_bbox_opt(False)
- at click.option('--bands/--mask', default=True,
- help="Extract shapes from one of the dataset bands or from "
- "its nodata mask")
- at click.option('--bidx', type=int, default=1,
- help="Index of the source band")
+ at click.option('--band/--mask', default=True,
+ help="Choose to extract from a band (the default) or a mask.")
+ at click.option('--bidx', 'bandidx', type=int, default=None,
+ help="Index of the band or mask that is the source of shapes.")
@click.option('--sampling', type=int, default=1,
- help="Inverse of the sampling fraction")
+ help="Inverse of the sampling fraction; "
+ "a value of 10 decimates.")
@click.option('--with-nodata/--without-nodata', default=False,
help="Include or do not include (the default) nodata regions.")
+ at click.option('--as-mask/--not-as-mask', default=False,
+ help="Interpret a band as a mask and output only one class of "
+ "valid data shapes.")
@click.pass_context
def shapes(
- ctx, input, precision, indent, compact, projection, sequence,
- use_rs, geojson_type, bands, bidx, sampling, with_nodata):
- """Writes features of a dataset out as GeoJSON. It's intended for
- use with single-band rasters and reads from the first band.
+ ctx, input, output, precision, indent, compact, projection, sequence,
+ use_rs, geojson_type, band, bandidx, sampling, with_nodata, as_mask):
+ """Extracts shapes from one band or mask of a dataset and writes
+ them out as GeoJSON. Unless otherwise specified, the shapes will be
+ transformed to WGS 84 coordinates.
+
+ The default action of this command is to extract shapes from the
+ first band of the input dataset. The shapes are polygons bounding
+ contiguous regions (or features) of the same raster value. This
+ command performs poorly for int16 or float type datasets.
+
+ Bands other than the first can be specified using the `--bidx`
+ option:
+
+ $ rio shapes --bidx 3 tests/data/RGB.byte.tif
+
+ The valid data footprint of a dataset's i-th band can be extracted
+ by using the `--mask` and `--bidx` options:
+
+ $ rio shapes --mask --bidx 1 tests/data/RGB.byte.tif
+
+ Omitting the `--bidx` option results in a footprint extracted from
+ the conjunction of all band masks. This is generally smaller than
+ any individual band's footprint.
+
+ A dataset band may be analyzed as though it were a binary mask with
+ the `--as-mask` option:
+
+ $ rio shapes --as-mask --bidx 1 tests/data/RGB.byte.tif
"""
# These import numpy, which we don't want to do unless it's needed.
import numpy
@@ -62,7 +233,11 @@ def shapes(
dump_kwds['indent'] = indent
if compact:
dump_kwds['separators'] = (',', ':')
- stdout = click.get_text_stream('stdout')
+
+ stdout = click.open_file(
+ output, 'w') if output else click.get_text_stream('stdout')
+
+ bidx = 1 if bandidx is None and band else bandidx
# This is the generator for (feature, bbox) pairs.
class Collection(object):
@@ -78,30 +253,58 @@ def shapes(
def __call__(self):
with rasterio.open(input) as src:
img = None
- nodata_mask = None
- if bands:
+ msk = None
+
+ # Adjust transforms.
+ if sampling == 1:
+ transform = src.affine
+ else:
+ transform = src.affine * Affine.scale(float(sampling))
+
+ # Most of the time, we'll use the valid data mask.
+ # We skip reading it if we're extracting every possible
+ # feature (even invalid data features) from a band.
+ if not band or (band and not as_mask and not with_nodata):
if sampling == 1:
- img = src.read_band(bidx)
- transform = src.transform
- # Decimate the band.
+ msk = src.read_masks(bidx)
else:
- img = numpy.zeros(
- (src.height//sampling, src.width//sampling),
- dtype=src.dtypes[src.indexes.index(bidx)])
- img = src.read_band(bidx, img)
- transform = src.affine * Affine.scale(float(sampling))
- if not bands or not with_nodata:
+ msk_shape = (
+ src.height//sampling, src.width//sampling)
+ if bidx is None:
+ msk = numpy.zeros(
+ (src.count,) + msk_shape, 'uint8')
+ else:
+ msk = numpy.zeros(msk_shape, 'uint8')
+ msk = src.read_masks(bidx, msk)
+
+ if bidx is None:
+ msk = numpy.logical_or.reduce(msk).astype('uint8')
+
+ # Possibly overidden below.
+ img = msk
+
+ # Read the band data unless the --mask option is given.
+ if band:
if sampling == 1:
- nodata_mask = src.read_mask()
- transform = src.transform
- # Decimate the mask.
+ img = src.read(bidx, masked=False)
else:
- nodata_mask = numpy.zeros(
+ img = numpy.zeros(
(src.height//sampling, src.width//sampling),
- dtype=numpy.uint8)
- nodata_mask = src.read_mask(nodata_mask)
- transform = src.affine * Affine.scale(float(sampling))
-
+ dtype=src.dtypes[src.indexes.index(bidx)])
+ img = src.read(bidx, img, masked=False)
+
+ # If --as-mask option was given, convert the image
+ # to a binary image. This reduces the number of shape
+ # categories to 2 and likely reduces the number of
+ # shapes.
+ if as_mask:
+ tmp = numpy.ones_like(img, 'uint8') * 255
+ tmp[img == 0] = 0
+ img = tmp
+ if not with_nodata:
+ msk = tmp
+
+ # Transform the raster bounds.
bounds = src.bounds
xs = [bounds[0], bounds[2]]
ys = [bounds[1], bounds[3]]
@@ -114,12 +317,12 @@ def shapes(
self._xs = xs
self._ys = ys
+ # Prepare keyword arguments for shapes().
kwargs = {'transform': transform}
- # Default is to exclude nodata features.
- if nodata_mask is not None:
- kwargs['mask'] = (nodata_mask==255)
- if img is None:
- img = nodata_mask
+ if not with_nodata:
+ kwargs['mask'] = msk
+
+ # Yield GeoJSON features.
for g, i in rasterio.features.shapes(img, **kwargs):
if projection == 'geographic':
g = rasterio.warp.transform_geom(
@@ -129,52 +332,48 @@ def shapes(
yield {
'type': 'Feature',
'id': str(i),
- 'properties': {'val': i},
+ 'properties': {
+ 'val': i, 'filename': os.path.basename(src.name)
+ },
'bbox': [min(xs), min(ys), max(xs), max(ys)],
- 'geometry': g }
+ 'geometry': g
+ }
if not sequence:
geojson_type = 'collection'
try:
- with rasterio.drivers(CPL_DEBUG=verbosity>2):
+ with rasterio.drivers(CPL_DEBUG=(verbosity > 2)):
write_features(
stdout, Collection(), sequence=sequence,
geojson_type=geojson_type, use_rs=use_rs,
**dump_kwds)
- sys.exit(0)
except Exception:
- logger.exception("Failed. Exception caught")
- sys.exit(1)
+ logger.exception("Exception caught during processing")
+ raise click.Abort()
# Rasterize command.
@cli.command(short_help='Rasterize features.')
@files_inout_arg
+ at output_opt
@format_opt
- at 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.')
- at click.option('--bounds', nargs=4, type=float, default=None,
- help='Output bounds: left, bottom, right, top.')
+ at like_file_opt
+ at bounds_opt
@click.option('--dimensions', nargs=2, type=int, default=None,
help='Output dataset width, height in number of pixels.')
- at 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')
- at click.option('--src_crs', default='EPSG:4326',
+ at resolution_opt
+ at click.option('--src-crs', '--src_crs', 'src_crs', default=None,
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')
- at click.option('--all_touched', is_flag=True, default=False)
- at click.option('--default_value', type=float, default=1, help='Default value '
- 'for rasterized pixels')
- at 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')
+ 'codes for now. Used as output coordinate system if output '
+ 'does not exist or --like option is not used. '
+ 'Default: EPSG:4326')
+ at all_touched_opt
+ at click.option('--default-value', '--default_value', 'default_value',
+ type=float, default=1, help='Default value for rasterized pixels')
+ at 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.')
@@ -182,6 +381,7 @@ def shapes(
def rasterize(
ctx,
files,
+ output,
driver,
like,
bounds,
@@ -197,24 +397,24 @@ def rasterize(
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.
+ coordinate reference system as the output unless --src-crs is provided.
--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.
+ transform and data type from that raster will be used to create the output.
+ Only a single band will be output.
+
+ The GeoJSON is assumed to be in the same coordinate reference system unless
+ --src-crs is provided.
--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
-
- The GeoJSON must be within the bounds of the existing output or --like
- raster, or an error will be raised.
+ --driver, --bounds, --dimensions, and --res are ignored when output exists
+ or --like raster is provided
If the output does not exist and --like raster is not provided, the input
@@ -233,16 +433,18 @@ def rasterize(
added in the future.
"""
- import numpy as np
+ from rasterio._base import is_geographic_crs, is_same_crs
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()
+ output, files = resolve_inout(files=files, output=output)
input = click.open_file(files.pop(0) if files else '-')
+ has_src_crs = src_crs is not None
+ src_crs = src_crs or 'EPSG:4326'
+
# 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):
@@ -257,12 +459,6 @@ def rasterize(
return feature['properties'].get(property, default_value)
return default_value
- 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])
-
geojson = json.loads(input.read())
if 'features' in geojson:
geometries = []
@@ -278,10 +474,16 @@ def rasterize(
if os.path.exists(output):
with rasterio.open(output, 'r+') as out:
- if disjoint_bounds(geojson_bounds, out.bounds):
- raise click.BadParameter('GeoJSON outside bounds of '
+ if has_src_crs and not is_same_crs(src_crs, out.crs):
+ raise click.BadParameter('GeoJSON does not match crs of '
'existing output raster',
- param=input, param_hint='input')
+ param='input', param_hint='input')
+
+ if _disjoint_bounds(geojson_bounds, out.bounds):
+ click.echo("GeoJSON outside bounds of existing output "
+ "raster. Are they in different coordinate "
+ "reference systems?",
+ err=True)
meta = out.meta.copy()
@@ -305,10 +507,17 @@ def rasterize(
else:
if like is not None:
template_ds = rasterio.open(like)
- if disjoint_bounds(geojson_bounds, template_ds.bounds):
- raise click.BadParameter('GeoJSON outside bounds of '
- '--like raster', param=input,
- param_hint='input')
+
+ if has_src_crs and not is_same_crs(src_crs, template_ds.crs):
+ raise click.BadParameter('GeoJSON does not match crs of '
+ '--like raster',
+ param='input', param_hint='input')
+
+ if _disjoint_bounds(geojson_bounds, template_ds.bounds):
+ click.echo("GeoJSON outside bounds of --like raster. "
+ "Are they in different coordinate reference "
+ "systems?",
+ err=True)
kwargs = template_ds.meta.copy()
kwargs['count'] = 1
@@ -317,11 +526,12 @@ def rasterize(
else:
bounds = bounds or geojson_bounds
- if src_crs == 'EPSG:4326':
+ if is_geographic_crs(src_crs):
if (bounds[0] < -180 or bounds[2] > 180 or
- bounds[1] < -80 or bounds[3] > 80):
+ bounds[1] < -80 or bounds[3] > 80):
raise click.BadParameter(
- 'Bounds are beyond the valid extent for EPSG:4326.',
+ "Bounds are beyond the valid extent for "
+ "EPSG:4326.",
ctx, param=bounds, param_hint='--bounds')
if dimensions:
@@ -376,3 +586,8 @@ def rasterize(
with rasterio.open(output, 'w', **kwargs) as out:
out.write_band(1, result)
+
+
+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 427f672..b627fdc 100644
--- a/rasterio/rio/info.py
+++ b/rasterio/rio/info.py
@@ -1,16 +1,108 @@
-# Info command.
+"""Fetch and edit raster dataset metadata from the command line."""
import json
import logging
-import os.path
-import pprint
import sys
import click
import rasterio
import rasterio.crs
-from rasterio.rio.cli import cli
+from rasterio.rio.cli import cli, bidx_opt, file_in_arg, masked_opt
+from rasterio.transform import guard_transform
+
+
+ at cli.command('edit-info', short_help="Edit dataset metadata.")
+ at file_in_arg
+ at click.option('--nodata', type=float, default=None,
+ help="New nodata value")
+ at click.option('--crs', help="New coordinate reference system")
+ at click.option('--transform', help="New affine transform matrix")
+ at click.option('--tag', 'tags', multiple=True, metavar='KEY=VAL',
+ help="New tag.")
+ at click.pass_context
+def edit(ctx, input, nodata, crs, transform, tags):
+ """Edit a dataset's metadata: coordinate reference system, affine
+ transformation matrix, nodata value, and tags.
+
+ CRS may be either a PROJ.4 or EPSG:nnnn string, or a JSON-encoded
+ PROJ.4 object.
+
+ Transforms are either JSON-encoded Affine objects (preferred) like
+
+ [300.038, 0.0, 101985.0, 0.0, -300.042, 2826915.0]
+
+ or JSON-encoded GDAL geotransform arrays like
+
+ [101985.0, 300.038, 0.0, 2826915.0, 0.0, -300.042]
+ """
+ import numpy as np
+
+ verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+ logger = logging.getLogger('rio')
+
+ def in_dtype_range(value, dtype):
+ infos = {'c': np.finfo, 'f': np.finfo, 'i': np.iinfo,
+ 'u': np.iinfo}
+ rng = infos[np.dtype(dtype).kind](dtype)
+ return rng.min <= value <= rng.max
+
+ with rasterio.drivers(CPL_DEBUG=(verbosity > 2)) as env:
+ with rasterio.open(input, 'r+') as dst:
+
+ # Update nodata.
+ if nodata is not None:
+
+ dtype = dst.dtypes[0]
+ if not in_dtype_range(nodata, dtype):
+ raise click.BadParameter(
+ "outside the range of the file's "
+ "data type (%s)." % dtype,
+ param=nodata, param_hint='nodata')
+
+ dst.nodata = nodata
+
+ # Update CRS. Value might be a PROJ.4 string or a JSON
+ # encoded dict.
+ if crs:
+ new_crs = crs.strip()
+ try:
+ new_crs = json.loads(crs)
+ except ValueError:
+ pass
+
+ if not (rasterio.crs.is_geographic_crs(new_crs) or
+ rasterio.crs.is_projected_crs(new_crs)):
+ raise click.BadParameter(
+ "'%s' is not a recognized CRS." % crs,
+ param=crs, param_hint='crs')
+
+ dst.crs = new_crs
+
+ # Update transform. Value might be a JSON encoded
+ # Affine object or a GDAL geotransform array.
+ if transform:
+ try:
+ transform_obj = json.loads(transform)
+ except ValueError:
+ raise click.BadParameter(
+ "'%s' is not a JSON array." % transform,
+ param=transform, param_hint='transform')
+
+ try:
+ transform_obj = guard_transform(transform_obj)
+ except:
+ raise click.BadParameter(
+ "'%s' is not recognized as an Affine or GDAL "
+ "geotransform array." % transform,
+ param=transform, param_hint='transform')
+
+ dst.transform = transform_obj
+
+ # Update tags.
+ if tags:
+ tags = dict(p.split('=') for p in tags)
+ dst.update_tags(**tags)
@cli.command(short_help="Print information about the rio environment.")
@@ -32,7 +124,7 @@ def env(ctx, key):
@cli.command(short_help="Print information about a data file.")
- at click.argument('input', type=click.Path(exists=True))
+ at file_in_arg
@click.option('--meta', 'aspect', flag_value='meta', default=True,
help="Show data file structure (default).")
@click.option('--tags', 'aspect', flag_value='tags',
@@ -44,7 +136,7 @@ def env(ctx, key):
# a string.
@click.option('--count', 'meta_member', flag_value='count',
help="Print the count of bands.")
- at click.option('--dtype', 'meta_member', flag_value='dtype',
+ at click.option('-t', '--dtype', 'meta_member', flag_value='dtype',
help="Print the dtype name.")
@click.option('--nodata', 'meta_member', flag_value='nodata',
help="Print the nodata value.")
@@ -61,7 +153,7 @@ def env(ctx, key):
@click.option('--bounds', 'meta_member', flag_value='bounds',
help="Print the boundary coordinates "
"(left, bottom, right, top).")
- at click.option('--res', 'meta_member', flag_value='res',
+ at click.option('-r', '--res', 'meta_member', flag_value='res',
help="Print pixel width and height.")
@click.option('--lnglat', 'meta_member', flag_value='lnglat',
help="Print longitude and latitude at center.")
@@ -70,10 +162,11 @@ def env(ctx, key):
"(use --bidx).")
@click.option('-v', '--tell-me-more', '--verbose', is_flag=True,
help="Output extra information.")
- at click.option('--bidx', type=int, default=1,
- help="Input file band index (default: 1).")
+ at bidx_opt
+ at masked_opt
@click.pass_context
-def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx):
+def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx,
+ masked):
"""Print metadata about the dataset as JSON.
Optionally print a single metadata item as a string.
@@ -81,12 +174,13 @@ def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx):
verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
logger = logging.getLogger('rio')
mode = 'r' if (verbose or meta_member == 'stats') else 'r-'
+
try:
with rasterio.drivers(CPL_DEBUG=(verbosity > 2)):
with rasterio.open(input, mode) as src:
info = src.meta
+ info['transform'] = info['affine'][:6]
del info['affine']
- del info['transform']
info['shape'] = info['height'], info['width']
info['bounds'] = src.bounds
proj4 = rasterio.crs.to_string(src.crs)
@@ -98,11 +192,12 @@ def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx):
if verbose:
stats = [{'min': float(b.min()),
'max': float(b.max()),
- 'mean': float(b.mean())} for b in src.read()]
+ 'mean': float(b.mean())
+ } for b in src.read(masked=masked)]
info['stats'] = stats
if aspect == 'meta':
if meta_member == 'stats':
- band = src.read(bidx)
+ band = src.read(bidx, masked=masked)
click.echo('%f %f %f' % (
float(band.min()),
float(band.max()),
@@ -115,9 +210,8 @@ def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx):
else:
click.echo(json.dumps(info, indent=indent))
elif aspect == 'tags':
- click.echo(json.dumps(src.tags(ns=namespace),
+ click.echo(json.dumps(src.tags(ns=namespace),
indent=indent))
- sys.exit(0)
except Exception:
- logger.exception("Failed. Exception caught")
- sys.exit(1)
+ logger.exception("Exception caught during processing")
+ raise click.Abort()
diff --git a/rasterio/rio/main.py b/rasterio/rio/main.py
index bd17c85..910f893 100644
--- a/rasterio/rio/main.py
+++ b/rasterio/rio/main.py
@@ -1,9 +1,32 @@
-#!/usr/bin/env python
-
-from rasterio.rio.cli import cli
-from rasterio.rio.bands import stack
-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
-from rasterio.rio.sample import sample
+# main: loader of all the command entry points.
+
+import sys
+import traceback
+
+from pkg_resources import iter_entry_points
+
+from rasterio.rio.cli import BrokenCommand, cli
+
+
+# Find and load all entry points in the rasterio.rio_commands group.
+# This includes the standard commands included with Rasterio as well
+# as commands provided by other packages.
+#
+# At a mimimum, commands must use the rasterio.rio.cli.cli command
+# group decorator like so:
+#
+# from rasterio.rio.cli import cli
+#
+# @cli.command()
+# def foo(...):
+# ...
+
+for entry_point in iter_entry_points('rasterio.rio_commands'):
+ try:
+ entry_point.load()
+ except Exception:
+ # Catch this so a busted plugin doesn't take down the CLI.
+ # Handled by registering a dummy command that does nothing
+ # other than explain the error.
+ cli.add_command(
+ BrokenCommand(entry_point.name))
diff --git a/rasterio/rio/merge.py b/rasterio/rio/merge.py
index dd9335b..a4c1afe 100644
--- a/rasterio/rio/merge.py
+++ b/rasterio/rio/merge.py
@@ -10,21 +10,21 @@ import click
from cligj import files_inout_arg, format_opt
import rasterio
-from rasterio.rio.cli import cli
+from rasterio.rio.cli import cli, bounds_opt, output_opt, resolve_inout
from rasterio.transform import Affine
@cli.command(short_help="Merge a stack of raster datasets.")
@files_inout_arg
+ at output_opt
@format_opt
- at click.option('--bounds', nargs=4, type=float, default=None,
- help="Output bounds: left, bottom, right, top.")
- at click.option('--res', nargs=2, type=float, default=None,
+ at bounds_opt
+ at click.option('-r', '--res', nargs=2, type=float, default=None,
help="Output dataset resolution: pixel width, pixel height")
- at click.option('--nodata', '-n', type=float, default=None,
+ at click.option('--nodata', type=float, default=None,
help="Override nodata values defined in input datasets")
@click.pass_context
-def merge(ctx, files, driver, bounds, res, nodata):
+def merge(ctx, files, output, driver, bounds, res, nodata):
"""Copy valid pixels from input files to an output file.
All files must have the same number of bands, data type, and
@@ -45,8 +45,7 @@ def merge(ctx, files, driver, bounds, res, nodata):
try:
with rasterio.drivers(CPL_DEBUG=verbosity>2):
- output = files[-1]
- files = files[:-1]
+ output, files = resolve_inout(files=files, output=output)
with rasterio.open(files[0]) as first:
first_res = first.res
@@ -92,10 +91,16 @@ def merge(ctx, files, driver, bounds, res, nodata):
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
@@ -119,32 +124,58 @@ def merge(ctx, files, driver, bounds, res, nodata):
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.
- window = src.window(*dst.bounds)
- data = np.zeros_like(dest)
- data = src.read(
- out=data,
- window=window,
- boundless=True,
+
+ # 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)
- np.copyto(dest, data,
+
+ # 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(
- dest==nodataval, data.mask==False))
+ region==nodataval, temp.mask==False))
if dst.mode == 'r+':
- data = dst.read(masked=True)
- np.copyto(dest, data,
+ temp = dst.read(masked=True)
+ np.copyto(dest, temp,
where=np.logical_and(
- dest==nodataval, data.mask==False))
+ dest==nodataval, temp.mask==False))
dst.write(dest)
dst.close()
- sys.exit(0)
except Exception:
- logger.exception("Failed. Exception caught")
- sys.exit(1)
+ logger.exception("Exception caught during processing")
+ raise click.Abort()
diff --git a/rasterio/rio/rio.py b/rasterio/rio/rio.py
index 78b1662..b8cd2b9 100644
--- a/rasterio/rio/rio.py
+++ b/rasterio/rio/rio.py
@@ -16,7 +16,7 @@ from cligj import (
geojson_type_feature_opt, geojson_type_bbox_opt)
import rasterio
-from rasterio.rio.cli import cli, write_features
+from rasterio.rio.cli import cli, write_features, file_in_arg
warnings.simplefilter('default')
@@ -30,8 +30,9 @@ warnings.simplefilter('default')
# Insp command.
@cli.command(short_help="Open a data file and start an interpreter.")
- at click.argument('input', type=click.Path(exists=True))
+ at file_in_arg
@click.option(
+ '-m',
'--mode',
type=click.Choice(['r', 'r+']),
default='r',
@@ -51,17 +52,16 @@ def insp(ctx, input, mode):
rasterio.__version__,
'.'.join(map(str, sys.version_info[:3]))),
src)
- sys.exit(0)
except Exception:
- logger.exception("Failed. Exception caught")
- sys.exit(1)
+ logger.exception("Exception caught during processing")
+ raise click.Abort()
# Bounds command.
@cli.command(short_help="Write bounding boxes to stdout as GeoJSON.")
# One or more files, the bounds of each are a feature in the collection
# object or feature sequence.
- at click.argument('input', nargs=-1, type=click.Path(exists=True))
+ at click.argument('INPUT', nargs=-1, type=click.Path(exists=True))
@precision_opt
@indent_opt
@compact_opt
@@ -132,7 +132,9 @@ def bounds(ctx, input, precision, indent, compact, projection, sequence,
[xs[0], ys[1]],
[xs[0], ys[0]] ]]},
'properties': {
- 'id': str(i), 'title': path} }
+ 'id': str(i),
+ 'title': path,
+ 'filename': os.path.basename(path)} }
self._xs.extend(bbox[::2])
self._ys.extend(bbox[1::2])
@@ -146,17 +148,17 @@ def bounds(ctx, input, precision, indent, compact, projection, sequence,
stdout, col, sequence=sequence,
geojson_type=geojson_type, use_rs=use_rs,
**dump_kwds)
- sys.exit(0)
+
except Exception:
- logger.exception("Failed. Exception caught")
- sys.exit(1)
+ logger.exception("Exception caught during processing")
+ raise click.Abort()
# Transform command.
@cli.command(short_help="Transform coordinates.")
- at click.argument('input', default='-', required=False)
- at click.option('--src_crs', default='EPSG:4326', help="Source CRS.")
- at click.option('--dst_crs', default='EPSG:4326', help="Destination CRS.")
+ at click.argument('INPUT', default='-', required=False)
+ at click.option('--src-crs', '--src_crs', default='EPSG:4326', help="Source CRS.")
+ at click.option('--dst-crs', '--dst_crs', default='EPSG:4326', help="Destination CRS.")
@precision_opt
@click.pass_context
def transform(ctx, input, src_crs, dst_crs, precision):
@@ -196,7 +198,6 @@ def transform(ctx, input, src_crs, dst_crs, precision):
result[1::2] = ys
print(json.dumps(result))
- sys.exit(0)
except Exception:
- logger.exception("Failed. Exception caught")
- sys.exit(1)
+ logger.exception("Exception caught during processing")
+ raise click.Abort()
diff --git a/rasterio/rio/sample.py b/rasterio/rio/sample.py
index fd4684a..96eef67 100644
--- a/rasterio/rio/sample.py
+++ b/rasterio/rio/sample.py
@@ -14,7 +14,7 @@ warnings.simplefilter('default')
@cli.command(short_help="Sample a dataset.")
@click.argument('files', nargs=-1, required=True, metavar='FILE "[x, y]"')
- at click.option('--bidx', default=None, help="Indexes of input file bands.")
+ at click.option('-b', '--bidx', default=None, help="Indexes of input file bands.")
@click.pass_context
def sample(ctx, files, bidx):
"""Sample a dataset at one or more points
@@ -88,7 +88,7 @@ def sample(ctx, files, bidx):
(json.loads(line) for line in points),
indexes=indexes):
click.echo(json.dumps(vals.tolist()))
- sys.exit(0)
+
except Exception:
- logger.exception("Failed. Exception caught")
- sys.exit(1)
+ logger.exception("Exception caught during processing")
+ raise click.Abort()
diff --git a/rasterio/tool.py b/rasterio/tool.py
index a7c2d15..6601cfc 100644
--- a/rasterio/tool.py
+++ b/rasterio/tool.py
@@ -29,7 +29,7 @@ def show(source, cmap='gray'):
tuple.
"""
if isinstance(source, tuple):
- arr = source[0].read_band(source[1])
+ arr = source[0].read(source[1])
else:
arr = source
if plt is not None:
@@ -43,7 +43,7 @@ def stats(source):
"""Return a tuple with raster min, max, and mean.
"""
if isinstance(source, tuple):
- arr = source[0].read_band(source[1])
+ arr = source[0].read(source[1])
else:
arr = source
return Stats(numpy.min(arr), numpy.max(arr), numpy.mean(arr))
diff --git a/rasterio/transform.py b/rasterio/transform.py
index 9b29945..c02bc1c 100644
--- a/rasterio/transform.py
+++ b/rasterio/transform.py
@@ -23,3 +23,18 @@ def guard_transform(transform):
else:
transform = Affine(*transform)
return transform
+
+
+def from_origin(west, north, xsize, ysize):
+ """Return an Affine transformation for a georeferenced raster given
+ the coordinates of its upper left corner `west`, `north` and pixel
+ sizes `xsize`, `ysize`."""
+ return Affine.translation(west, north) * Affine.scale(xsize, -ysize)
+
+
+def from_bounds(west, south, east, north, width, height):
+ """Return an Affine transformation for a georeferenced raster given
+ its bounds `west`, `south`, `east`, `north` and its `width` and
+ `height` in number of pixels."""
+ return Affine.translation(west, north) * Affine.scale(
+ (east - west)/width, (south - north)/height)
diff --git a/rasterio/warp.py b/rasterio/warp.py
index dd971d8..219eac7 100644
--- a/rasterio/warp.py
+++ b/rasterio/warp.py
@@ -1,5 +1,9 @@
"""Raster warping and reprojection"""
+from affine import Affine
+from math import ceil
+import numpy as np
+
from rasterio._base import _transform
from rasterio._warp import _transform_geom, _reproject, RESAMPLING
from rasterio.transform import guard_transform
@@ -77,18 +81,88 @@ def transform_geom(
precision)
+def transform_bounds(left, bottom, right, top, src_crs, dst_crs, densify_pts=21):
+ """
+ Transforms bounds from src_crs to dst_crs, optionally densifying the edges
+ (to account for nonlinear transformations along these edges) and extracting
+ the outermost bounds.
+
+ Note: this does not account for the antimeridian.
+
+ Parameters
+ ----------
+ left, bottom, right, top: float
+ Bounding coordinates in src_crs, from the bounds property of a raster.
+ src_crs: dict
+ Source coordinate reference system, in rasterio dict format.
+ Example: {'init': 'EPSG:4326'}
+ dst_crs: dict
+ Target coordinate reference system.
+ densify_pts: uint, optional
+ Number of points to add to each edge to account for nonlinear
+ edges produced by the transform process. Large numbers will produce
+ worse performance. Default: 21 (gdal default).
+
+ Returns
+ -------
+ left, bottom, right, top: float
+ Outermost coordinates in target coordinate reference system.
+ """
+
+ if densify_pts < 0:
+ raise ValueError('densify parameter must be >= 0')
+
+ in_xs = []
+ in_ys = []
+
+ if densify_pts > 0:
+ densify_factor = 1.0 / float(densify_pts + 1)
+
+ # Add points along outer edges.
+ for x in (left, right):
+ in_xs.extend([x] * (densify_pts + 2))
+ in_ys.extend(
+ bottom + np.arange(0, densify_pts + 2, dtype=np.float32)
+ * ((top - bottom) * densify_factor)
+ )
+
+ for y in (bottom, top):
+ in_xs.extend(
+ left + np.arange(1, densify_pts + 1, dtype=np.float32)
+ * ((right - left) * densify_factor)
+ )
+ in_ys.extend([y] * densify_pts)
+
+ else:
+ in_xs = [left, left, right, right]
+ in_ys = [bottom, top, bottom, top]
+
+ xs, ys = transform(src_crs, dst_crs, in_xs, in_ys)
+ return (min(xs), min(ys), max(xs), max(ys))
+
+
def reproject(
source,
destination,
src_transform=None,
src_crs=None,
+ src_nodata=None,
dst_transform=None,
dst_crs=None,
+ dst_nodata=None,
resampling=RESAMPLING.nearest,
**kwargs):
"""
Reproject a source raster to a destination raster.
+ If the source and destination are ndarrays, coordinate reference
+ system definitions and affine transformation parameters are required
+ for reprojection.
+
+ If the source and destination are rasterio Bands, shorthand for
+ bands of datasets on disk, the coordinate reference systems and
+ transforms will be read from the appropriate datasets.
+
Parameters
------------
source: ndarray or rasterio Band
@@ -103,12 +177,22 @@ def reproject(
Required if source and destination are ndarrays.
Will be derived from source if it is a rasterio Band.
Example: {'init': 'EPSG:4326'}
+ src_nodata: int or float, optional
+ The source nodata value. Pixels with this value will not be used
+ for interpolation. If not set, it will be default to the
+ nodata value of the source image if a masked ndarray or rasterio band,
+ if available. Must be provided if dst_nodata is not None.
dst_transform: affine transform object, optional
Target affine transformation. Required if source and destination
are ndarrays. Will be derived from target if it is a rasterio Band.
dst_crs: dict, optional
Target coordinate reference system. Required if source and destination
are ndarrays. Will be derived from target if it is a rasterio Band.
+ dst_nodata: int or float, optional
+ The nodata value used to initialize the destination; it will remain
+ in all areas not covered by the reprojected source. Defaults to the
+ nodata value of the destination image (if set), the value of
+ src_nodata, or 0 (GDAL default).
resampling: int
Resampling method to use. One of the following:
RESAMPLING.nearest,
@@ -137,7 +221,81 @@ def reproject(
destination,
src_transform,
src_crs,
+ src_nodata,
dst_transform,
dst_crs,
+ dst_nodata,
resampling,
**kwargs)
+
+
+def calculate_default_transform(
+ left,
+ bottom,
+ right,
+ top,
+ width,
+ height,
+ src_crs,
+ dst_crs,
+ resolution=None,
+ densify_pts=21):
+ """
+ Transforms bounds to destination coordinate system, calculates resolution
+ if not provided, and returns destination transform and dimensions.
+ Intended to be used to calculate parameters for reproject function.
+
+ Destination transform is anchored from the left, top coordinate.
+
+ Destination width and height are calculated from the number of pixels on
+ each dimension required to fit the destination bounds.
+
+ If resolution is not provided, it is calculated using a weighted average
+ of the relative sizes of source width and height compared to the transformed
+ bounds (pixels are assumed to be square).
+
+
+ Parameters
+ ----------
+ left, bottom, right, top: float
+ Bounding coordinates in src_crs, from the bounds property of a raster.
+ src_crs: dict
+ Source coordinate reference system, in rasterio dict format.
+ Example: {'init': 'EPSG:4326'}
+ dst_crs: dict
+ Target coordinate reference system.
+ resolution: tuple (x resolution, y resolution) or float, optional
+ Target resolution, in units of target coordinate reference system.
+ densify_pts: uint, optional
+ Number of points to add to each edge to account for nonlinear
+ edges produced by the transform process. Large numbers will produce
+ worse performance. Default: 21 (gdal default).
+
+ Returns
+ -------
+ tuple of destination affine transform, width, and height
+ """
+
+ xmin, ymin, xmax, ymax = transform_bounds(
+ left, bottom, right, top, src_crs, dst_crs, densify_pts)
+
+ x_dif = xmax - xmin
+ y_dif = ymax - ymin
+ size = float(width + height)
+
+ if resolution is None:
+ # TODO: compare to gdalwarp default
+ avg_resolution = (
+ (x_dif / float(width)) * (float(width) / size) +
+ (y_dif / float(height)) * (float(height) / size)
+ )
+ resolution = (avg_resolution, avg_resolution)
+
+ elif not isinstance(resolution, (tuple, list)):
+ resolution = (resolution, resolution)
+
+ dst_affine = Affine(resolution[0], 0, xmin, 0, -resolution[1], ymax)
+ dst_width = max(int(ceil(x_dif / resolution[0])), 1)
+ dst_height = max(int(ceil(y_dif / resolution[1])), 1)
+
+ return dst_affine, dst_width, dst_height
\ No newline at end of file
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 4d9ff53..69196b1 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,10 +1,11 @@
affine
cligj
coveralls>=0.4
-cython>=0.20
+cython==0.21.2
delocate
enum34
numpy>=1.8.0
+snuggs>=1.2
pytest
setuptools>=0.9.8
wheel
diff --git a/requirements.txt b/requirements.txt
index 66a3935..0b59530 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,4 +2,5 @@ affine
cligj
enum34
numpy>=1.8.0
+snuggs>=1.2
setuptools
diff --git a/setup.py b/setup.py
index 1788e52..7bed441 100755
--- a/setup.py
+++ b/setup.py
@@ -26,7 +26,27 @@ log = logging.getLogger()
if 'all' in sys.warnoptions:
log.level = logging.DEBUG
-# Parse the version from the fiona module.
+def check_output(cmd):
+ # since subprocess.check_output doesn't exist in 2.6
+ # we wrap it here.
+ try:
+ out = subprocess.check_output(cmd)
+ return out.decode('utf')
+ except AttributeError:
+ # For some reasone check_output doesn't exist
+ # So fall back on Popen
+ p = subprocess.Popen(cmd, stdout=subprocess.PIPE)
+ out, err = p.communicate()
+ return out
+
+def copy_data_tree(datadir, destdir):
+ try:
+ shutil.rmtree(destdir)
+ except OSError:
+ pass
+ shutil.copytree(datadir, destdir)
+
+# Parse the version from the rasterio module.
with open('rasterio/__init__.py') as f:
for line in f:
if line.find("__version__") >= 0:
@@ -50,6 +70,7 @@ include_dirs = []
library_dirs = []
libraries = []
extra_link_args = []
+gdal_output = [None]*3
try:
import numpy
@@ -60,18 +81,13 @@ except ImportError:
try:
gdal_config = os.environ.get('GDAL_CONFIG', 'gdal-config')
- with open("gdal-config.txt", "w") as gcfg:
- subprocess.call([gdal_config, "--cflags"], stdout=gcfg)
- subprocess.call([gdal_config, "--libs"], stdout=gcfg)
- subprocess.call([gdal_config, "--datadir"], stdout=gcfg)
- with open("gdal-config.txt", "r") as gcfg:
- cflags = gcfg.readline().strip()
- libs = gcfg.readline().strip()
- datadir = gcfg.readline().strip()
- for item in cflags.split():
+ for i, flag in enumerate(("--cflags", "--libs", "--datadir")):
+ gdal_output[i] = check_output([gdal_config, flag]).strip()
+
+ for item in gdal_output[0].split():
if item.startswith("-I"):
include_dirs.extend(item[2:].split(":"))
- for item in libs.split():
+ for item in gdal_output[1].split():
if item.startswith("-L"):
library_dirs.extend(item[2:].split(":"))
elif item.startswith("-l"):
@@ -80,27 +96,33 @@ try:
# e.g. -framework GDAL
extra_link_args.append(item)
- # Conditionally copy the GDAL data. To be used in conjunction with
- # the bdist_wheel command to make self-contained binary wheels.
- if os.environ.get('PACKAGE_DATA'):
- try:
- shutil.rmtree('rasterio/gdal_data')
- except OSError:
- pass
- shutil.copytree(datadir, 'rasterio/gdal_data')
-
except Exception as e:
- log.warning("Failed to get options via gdal-config: %s", str(e))
+ if os.name == "nt":
+ log.info(("Building on Windows requires extra options to setup.py to locate needed GDAL files.\n"
+ "More information is available in the README."))
+ else:
+ log.warning("Failed to get options via gdal-config: %s", str(e))
+
-# Conditionally copy PROJ.4 data.
+# Conditionally copy the GDAL data. To be used in conjunction with
+# the bdist_wheel command to make self-contained binary wheels.
if os.environ.get('PACKAGE_DATA'):
+ destdir = 'rasterio/gdal_data'
+ if gdal_output[2]:
+ log.info("Copying gdal data from %s" % gdal_output[2])
+ copy_data_tree(gdal_output[2], destdir)
+ else:
+ # check to see if GDAL_DATA is defined
+ gdal_data = os.environ.get('GDAL_DATA', None)
+ if gdal_data:
+ log.info("Copying gdal_data from %s" % gdal_data)
+ copy_data_tree(gdal_data, destdir)
+
+ # Conditionally copy PROJ.4 data.
projdatadir = os.environ.get('PROJ_LIB', '/usr/local/share/proj')
if os.path.exists(projdatadir):
- try:
- shutil.rmtree('rasterio/proj_data')
- except OSError:
- pass
- shutil.copytree(projdatadir, 'rasterio/proj_data')
+ log.info("Copying proj_data from %s" % projdatadir)
+ copy_data_tree(projdatadir, 'rasterio/proj_data')
ext_options = dict(
include_dirs=include_dirs,
@@ -169,7 +191,8 @@ with open('README.rst') as f:
inst_reqs = [
'affine>=1.0',
'cligj',
- 'Numpy>=1.7' ]
+ 'Numpy>=1.7',
+ 'snuggs>=1.3.1']
if sys.version_info < (3, 4):
inst_reqs.append('enum34')
@@ -202,6 +225,21 @@ setup_args = dict(
entry_points='''
[console_scripts]
rio=rasterio.rio.main:cli
+
+ [rasterio.rio_commands]
+ bounds=rasterio.rio.rio:bounds
+ calc=rasterio.rio.calc:calc
+ edit-info=rasterio.rio.info:edit
+ env=rasterio.rio.info:env
+ info=rasterio.rio.info:info
+ insp=rasterio.rio.rio:insp
+ mask=rasterio.rio.features:mask
+ merge=rasterio.rio.merge:merge
+ rasterize=rasterio.rio.features:rasterize
+ sample=rasterio.rio.sample:sample
+ shapes=rasterio.rio.features:shapes
+ stack=rasterio.rio.bands:stack
+ transform=rasterio.rio.rio:transform
''',
include_package_data=True,
ext_modules=ext_modules,
diff --git a/tests/conftest.py b/tests/conftest.py
index c5a9933..90aed75 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,10 +1,12 @@
import functools
import operator
import os
+import shutil
import sys
-import pytest
from click.testing import CliRunner
+import py
+import pytest
if sys.version_info > (3,):
@@ -20,10 +22,19 @@ def pytest_cmdline_main(config):
if reduce(operator.and_, map(os.path.exists, test_files)):
print("Test data present.")
else:
- print("Test data not present. See download directions in tests/README.txt")
+ print("Test data not present. See download directions in tests/data/README.rst")
sys.exit(1)
@pytest.fixture(scope='function')
def runner():
return CliRunner()
+
+
+ at pytest.fixture(scope='function')
+def data():
+ """A temporary directory containing a copy of the files in data."""
+ tmpdir = py.test.ensuretemp('tests/data')
+ for filename in test_files:
+ shutil.copy(filename, str(tmpdir))
+ return tmpdir
diff --git a/tests/test_band_masks.py b/tests/test_band_masks.py
new file mode 100644
index 0000000..36dd382
--- /dev/null
+++ b/tests/test_band_masks.py
@@ -0,0 +1,100 @@
+import logging
+import sys
+
+import numpy as np
+from pytest import fixture
+
+import rasterio
+
+
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+
+def test_masks():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ rm, gm, bm = src.read_masks()
+ r, g, b = src.read(masked=False)
+ assert not r[rm==0].any()
+ assert not g[gm==0].any()
+ assert not b[bm==0].any()
+
+
+def test_masked_true():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ r, g, b = src.read(masked=True)
+ rm, gm, bm = src.read_masks()
+ assert (r.mask==~rm.astype('bool')).all()
+ assert (g.mask==~gm.astype('bool')).all()
+ assert (b.mask==~bm.astype('bool')).all()
+
+
+def test_masked_none():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ r, g, b = src.read(masked=True)
+ rm, gm, bm = src.read_masks()
+ assert (r.mask==~rm.astype('bool')).all()
+ assert (g.mask==~gm.astype('bool')).all()
+ assert (b.mask==~bm.astype('bool')).all()
+
+
+ at fixture(scope='function')
+def tiffs(tmpdir):
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ kwds = src.meta
+
+ del kwds['nodata']
+ with rasterio.open(
+ str(tmpdir.join('no-nodata.tif')), 'w',
+ **kwds) as dst:
+ dst.write(src.read(masked=False))
+
+ with rasterio.open(
+ str(tmpdir.join('sidecar-masked.tif')), 'w',
+ **kwds) as dst:
+ dst.write(src.read(masked=False))
+ mask = np.zeros(src.shape, dtype='uint8')
+ dst.write_mask(mask)
+
+ return tmpdir
+
+
+def test_masking_no_nodata(tiffs):
+ # if the dataset has no defined nodata values, all data is
+ # considered valid data. The GDAL masks bands are arrays of
+ # 255 values. ``read()`` returns masked arrays where `mask`
+ # is False.
+ filename = str(tiffs.join('no-nodata.tif'))
+ with rasterio.open(filename) as src:
+
+ rgb = src.read(masked=False)
+ assert not hasattr(rgb, 'mask')
+ r = src.read(1, masked=False)
+ assert not hasattr(r, 'mask')
+
+ rgb = src.read(masked=True)
+ assert hasattr(rgb, 'mask')
+ assert not rgb.mask.any()
+ r = src.read(1, masked=True)
+ assert hasattr(r, 'mask')
+ assert not r.mask.any()
+
+ rgb = src.read(masked=True)
+ assert hasattr(rgb, 'mask')
+ assert not r.mask.any()
+ r = src.read(1, masked=True)
+ assert not r.mask.any()
+
+ masks = src.read_masks()
+ assert masks.all()
+
+
+def test_masking_sidecar_mask(tiffs):
+ # If the dataset has a .msk sidecar mask band file, all masks will
+ # be derived from that file.
+ with rasterio.open(str(tiffs.join('sidecar-masked.tif'))) as src:
+ rgb = src.read(masked=True)
+ assert rgb.mask.all()
+ r = src.read(1, masked=True)
+ assert r.mask.all()
+ masks = src.read_masks()
+ assert not masks.any()
diff --git a/tests/test_coords.py b/tests/test_coords.py
index 69eafcf..7a2af5c 100644
--- a/tests/test_coords.py
+++ b/tests/test_coords.py
@@ -1,4 +1,3 @@
-
import rasterio
def test_bounds():
@@ -17,18 +16,3 @@ def test_ul():
def test_res():
with rasterio.open('tests/data/RGB.byte.tif') as src:
assert tuple(round(v, 6) for v in src.res) == (300.037927, 300.041783)
-
-def test_index():
- 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)
- assert src.index(101985.0+400.0, 2826915.0-700.0) == (2, 1)
-
-def test_window():
- with rasterio.open('tests/data/RGB.byte.tif') as src:
- left, bottom, right, top = src.bounds
- assert src.window(left, bottom, right, top) == ((0, src.height),
- (0, src.width))
- assert src.window(left, top-400, left+400, top) == ((0, 1), (0, 1))
- assert src.window(left, top-500, left+500, top) == ((0, 2), (0, 2))
-
diff --git a/tests/test_crs.py b/tests/test_crs.py
index 4bb549b..eb6da86 100644
--- a/tests/test_crs.py
+++ b/tests/test_crs.py
@@ -5,6 +5,7 @@ import sys
import rasterio
from rasterio import crs
+from rasterio._base import is_geographic_crs, is_projected_crs, is_same_crs
logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
@@ -56,6 +57,15 @@ def test_write_3857(tmpdir):
AUTHORITY["EPSG","3857"]]""" in info.decode('utf-8')
+def test_from_epsg():
+ crs_dict = crs.from_epsg(4326)
+ assert crs_dict['init'].lower() == 'epsg:4326'
+
+ # Test with invalid EPSG code
+ with pytest.raises(ValueError):
+ assert crs.from_epsg(0)
+
+
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,
@@ -64,3 +74,47 @@ def test_bare_parameters():
# Example produced by pyproj
crs_dict = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
assert crs_dict.get('no_defs', False) is True
+
+ crs_dict = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=False +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
+ assert crs_dict.get('no_defs', True) is False
+
+
+def test_is_geographic():
+ assert is_geographic_crs({'init': 'EPSG:4326'}) is True
+ assert is_geographic_crs({'init': 'EPSG:3857'}) is False
+
+ wgs84_crs = crs.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
+ assert is_geographic_crs(wgs84_crs) is True
+
+ nad27_crs = crs.from_string('+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs')
+ assert is_geographic_crs(nad27_crs) is True
+
+ lcc_crs = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
+ assert is_geographic_crs(lcc_crs) is False
+
+
+def test_is_projected():
+ assert is_projected_crs({'init': 'EPSG:3857'}) is True
+ assert is_projected_crs({'init': 'EPSG:4326'}) is False
+
+ lcc_crs = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
+ assert is_projected_crs(lcc_crs) is True
+
+ wgs84_crs = crs.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
+ assert is_projected_crs(wgs84_crs) is False
+
+
+def test_is_same_crs():
+ crs1 = {'init': 'EPSG:4326'}
+ crs2 = {'init': 'EPSG:3857'}
+
+ assert is_same_crs(crs1, crs1) is True
+ assert is_same_crs(crs1, crs2) is False
+
+ wgs84_crs = crs.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
+ assert is_same_crs(crs1, wgs84_crs) is True
+
+ # Make sure that same projection with different parameter are not equal
+ lcc_crs1 = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
+ lcc_crs2 = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=45 +lat_0=0')
+ assert is_same_crs(lcc_crs1, lcc_crs2) is False
\ No newline at end of file
diff --git a/tests/test_features_rasterize.py b/tests/test_features_rasterize.py
index 2165f06..8848a8f 100644
--- a/tests/test_features_rasterize.py
+++ b/tests/test_features_rasterize.py
@@ -4,7 +4,7 @@ import numpy
import pytest
import rasterio
-from rasterio.features import shapes, rasterize
+from rasterio.features import shapes, rasterize, geometry_mask
logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
@@ -159,3 +159,23 @@ def test_rasterize_geometries_symmetric():
s = shapes(truth, transform=transform)
result = rasterize(s, out_shape=(rows, cols), transform=transform)
assert numpy.array_equal(result, truth)
+
+
+def test_geometry_mask():
+ rows = cols = 10
+ transform = (1.0, 0.0, 0.0, 0.0, -1.0, 0.0)
+ truth = numpy.zeros((rows, cols), dtype=rasterio.bool_)
+ truth[2:5, 2:5] = True
+ with rasterio.drivers():
+ s = shapes((truth * 10).astype(rasterio.ubyte), transform=transform)
+ # Strip out values returned from shapes, and only keep first shape
+ geoms = [next(s)[0]]
+
+ # Regular mask should be the inverse of truth raster
+ mask = geometry_mask(geoms, out_shape=(rows, cols), transform=transform)
+ assert numpy.array_equal(mask, numpy.invert(truth))
+
+ # Inverted mask should be the same as the truth raster
+ mask = geometry_mask(geoms, out_shape=(rows, cols), transform=transform,
+ invert=True)
+ assert numpy.array_equal(mask, truth)
\ No newline at end of file
diff --git a/tests/test_features_shapes.py b/tests/test_features_shapes.py
index 4ab9d78..2419a15 100644
--- a/tests/test_features_shapes.py
+++ b/tests/test_features_shapes.py
@@ -62,12 +62,21 @@ def test_shapes_connectivity():
image[5:11, 5:11] = 1
image[11, 11] = 1
+ shapes = ftrz.shapes(image, connectivity=4)
+ shape, val = next(shapes)
+ assert len(shape['coordinates'][0]) == 5
+
shapes = ftrz.shapes(image, connectivity=8)
shape, val = next(shapes)
assert len(shape['coordinates'][0]) == 9
# Note: geometry is not technically valid at this point, it has a self
# intersection at 11,11
+ # Test invalid connectivity
+ with pytest.raises(ValueError):
+ shapes = ftrz.shapes(image, connectivity=12)
+ next(shapes)
+
def test_shapes_dtype():
"""Test image data type handling"""
@@ -107,4 +116,29 @@ def test_shapes_dtype():
image = numpy.zeros((rows, cols), dtype=dtype)
image[2:5, 2:5] = test_value
shapes = ftrz.shapes(image)
- shape, value = next(shapes)
+ next(shapes)
+
+ # Test mask types
+ image = numpy.zeros((rows, cols), dtype='uint8')
+ image.fill(255)
+ supported_mask_types = (
+ ('bool', 1),
+ ('uint8', 255)
+ )
+ for dtype, mask_value in supported_mask_types:
+ mask = numpy.zeros((rows, cols), dtype=dtype)
+ mask[2:5, 2:5] = mask_value
+ shapes = ftrz.shapes(image, mask=mask)
+ shape, value = next(shapes)
+ assert value == 255
+
+ unsupported_mask_types = (
+ ('int8', -127),
+ ('int16', -32768)
+ )
+ for dtype, mask_value in unsupported_mask_types:
+ with pytest.raises(ValueError):
+ mask = numpy.zeros((rows, cols), dtype=dtype)
+ mask[2:5, 2:5] = mask_value
+ shapes = ftrz.shapes(image, mask=mask)
+ next(shapes)
\ No newline at end of file
diff --git a/tests/test_features_sieve.py b/tests/test_features_sieve.py
index baea580..e34cf34 100644
--- a/tests/test_features_sieve.py
+++ b/tests/test_features_sieve.py
@@ -27,6 +27,11 @@ def test_sieve():
sieved_image = ftrz.sieve(image, 101)
assert not sieved_image.any()
+ # Invalid size value should fail
+ for invalid_size in (0, 45.1234, image.size + 1):
+ with pytest.raises(ValueError):
+ sieved_image = ftrz.sieve(image, invalid_size)
+
def test_sieve_connectivity():
"""Test proper behavior of connectivity"""
@@ -46,6 +51,10 @@ def test_sieve_connectivity():
sieved_image = ftrz.sieve(image, 54, connectivity=8)
assert numpy.array_equal(sieved_image, image)
+ # Invalid connectivity value should fail
+ with pytest.raises(ValueError):
+ ftrz.sieve(image, 54, connectivity=12)
+
def test_sieve_output():
"""Test proper behavior of output image, if passed into sieve"""
@@ -66,6 +75,11 @@ def test_sieve_output():
with pytest.raises(ValueError):
ftrz.sieve(image, 100, output)
+ # Output of a different shape should fail
+ output = numpy.zeros((21, 21), dtype=rasterio.ubyte)
+ with pytest.raises(ValueError):
+ ftrz.sieve(image, 100, output)
+
def test_sieve_mask():
"""Test proper behavior of mask image, if passed int sieve"""
@@ -92,9 +106,8 @@ def test_sieve_mask():
truth[7:10, 7:10] = 1
assert numpy.array_equal(sieved_image, truth)
- # mask of other type than rasterio.bool_ should fail
- mask = numpy.zeros(shape, dtype=rasterio.uint8)
with pytest.raises(ValueError):
+ mask = numpy.ones((21, 21), dtype=rasterio.bool_)
ftrz.sieve(image, 100, mask=mask)
@@ -138,4 +151,35 @@ def test_dtypes():
with pytest.raises(ValueError):
image = numpy.zeros((rows, cols), dtype=dtype)
image[2:5, 2:5] = test_value
- sieved_image = ftrz.sieve(image, 2)
+ ftrz.sieve(image, 2)
+
+ # Test mask types
+ image = numpy.zeros((rows, cols), dtype='uint8')
+ image.fill(255)
+
+ supported_mask_types = (
+ ('bool', 1),
+ ('uint8', 255)
+ )
+ for dtype, mask_value in supported_mask_types:
+ mask = numpy.zeros((rows, cols), dtype=dtype)
+ mask[2:5, 2:5] = mask_value
+ sieved_image = ftrz.sieve(image, 2, mask=mask)
+ assert numpy.array_equal(image, sieved_image)
+
+ unsupported_mask_types = (
+ ('int8', -127),
+ ('int16', -32768)
+ )
+ for dtype, mask_value in unsupported_mask_types:
+ with pytest.raises(ValueError):
+ mask = numpy.zeros((rows, cols), dtype=dtype)
+ mask[2:5, 2:5] = mask_value
+ ftrz.sieve(image, 2, mask=mask)
+
+
+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_indexing.py b/tests/test_indexing.py
index 317ddb3..f536302 100644
--- a/tests/test_indexing.py
+++ b/tests/test_indexing.py
@@ -1,6 +1,6 @@
-
import rasterio
+
def test_index():
with rasterio.open('tests/data/RGB.byte.tif') as src:
left, bottom, right, top = src.bounds
@@ -9,13 +9,36 @@ def test_index():
assert src.index(right, bottom) == (src.height, src.width)
assert src.index(left, bottom) == (src.height, 0)
+
def test_full_window():
with rasterio.open('tests/data/RGB.byte.tif') as src:
- assert src.window(*src.bounds) == tuple(zip((0, 0), src.shape))
+ left, bottom, right, top = src.bounds
+ assert src.window(left, bottom, right, top) == tuple(zip((0, 0), src.shape))
+
def test_window_no_exception():
with rasterio.open('tests/data/RGB.byte.tif') as src:
left, bottom, right, top = src.bounds
left -= 1000.0
- assert src.window(left, bottom, right, top) == (
- (0, src.height), (-3, src.width))
+ assert src.window(left, bottom, right, top, boundless=True) == (
+ (0, src.height), (-4, src.width))
+
+
+def test_index():
+ 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)
+ assert src.index(101985.0+400.0, 2826915.0-700.0) == (2, 1)
+
+
+def test_window():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ left, bottom, right, top = src.bounds
+ dx, dy = src.res
+ eps = 1.0e-8
+ assert src.window(left+eps, bottom+eps, right-eps, top-eps) == ((0, src.height-1),
+ (0, src.width-1))
+ 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))
diff --git a/tests/test_mask_creation.py b/tests/test_mask_creation.py
new file mode 100644
index 0000000..a6f8c3f
--- /dev/null
+++ b/tests/test_mask_creation.py
@@ -0,0 +1,45 @@
+"""
+Tests of band mask creation, both .msk sidecar and internal.
+
+See https://github.com/mapbox/rasterio/issues/293 for bug report.
+"""
+
+import rasterio
+
+
+def test_create_internal_mask(data):
+ """Write an internal mask to the fixture's RGB.byte.tif."""
+ with rasterio.drivers(GDAL_TIFF_INTERNAL_MASK=True):
+ with rasterio.open(str(data.join('RGB.byte.tif')), 'r+') as dst:
+ blue = dst.read(1, masked=False)
+ mask = 255*(blue == 0).astype('uint8')
+ dst.write_mask(mask)
+
+ # There should be no .msk file
+ assert data.join('RGB.byte.tif').exists()
+ assert not data.join('RGB.byte.tif.msk').exists()
+
+ # Check that the mask was saved correctly.
+ with rasterio.open(str(data.join('RGB.byte.tif'))) as src:
+ assert (mask == src.read_mask()).all()
+
+
+def test_create_sidecar_mask(data):
+ """Write a .msk sidecar mask."""
+ with rasterio.drivers():
+ with rasterio.open(str(data.join('RGB.byte.tif')), 'r+') as dst:
+ blue = dst.read(1, masked=False)
+ mask = 255*(blue == 0).astype('uint8')
+ dst.write_mask(mask)
+
+ # There should be a .msk file in this case.
+ assert data.join('RGB.byte.tif').exists()
+ assert data.join('RGB.byte.tif.msk').exists()
+
+ # Check that the mask was saved correctly.
+ with rasterio.open(str(data.join('RGB.byte.tif'))) as src:
+ assert (mask == src.read_mask()).all()
+
+ # Check the .msk file, too.
+ with rasterio.open(str(data.join('RGB.byte.tif.msk'))) as msk:
+ assert (mask == msk.read(1, masked=False)).all()
diff --git a/tests/test_meta.py b/tests/test_meta.py
new file mode 100644
index 0000000..5e59c39
--- /dev/null
+++ b/tests/test_meta.py
@@ -0,0 +1,28 @@
+# Tests of dataset meta keywords and dataset creation
+
+import rasterio
+
+
+def test_blocksize_rgb(tmpdir):
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ kwds = src.meta
+ assert kwds['blockxsize'] == 791
+ assert kwds['blockysize'] == 3
+ assert kwds['tiled'] is False
+
+def test_blocksize_shade(tmpdir):
+ with rasterio.open('tests/data/shade.tif') as src:
+ kwds = src.meta
+ assert kwds['blockxsize'] == 1024
+ assert kwds['blockysize'] == 8
+ assert kwds['tiled'] is False
+
+def test_copy_meta(tmpdir):
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ kwds = src.meta
+ with rasterio.open(
+ str(tmpdir.join('test_copy_meta.tif')), 'w', **kwds) as dst:
+ assert dst.meta['count'] == 3
+ assert dst.meta['blockxsize'] == 791
+ assert dst.meta['blockysize'] == 3
+ assert dst.meta['tiled'] is False
diff --git a/tests/test_nodata.py b/tests/test_nodata.py
index 1fee860..98d02c6 100644
--- a/tests/test_nodata.py
+++ b/tests/test_nodata.py
@@ -13,6 +13,7 @@ def test_nodata(tmpdir):
with rasterio.drivers():
with rasterio.open('tests/data/RGB.byte.tif') as src:
with rasterio.open(dst_path, 'w', **src.meta) as dst:
+ assert dst.nodata == 0.0
assert dst.meta['nodata'] == 0.0
assert dst.nodatavals == [0.0, 0.0, 0.0]
info = subprocess.check_output([
@@ -31,6 +32,7 @@ def test_set_nodata(tmpdir):
meta = src.meta
meta['nodata'] = 42
with rasterio.open(dst_path, 'w', **meta) as dst:
+ assert dst.nodata == 42
assert dst.meta['nodata'] == 42
assert dst.nodatavals == [42, 42, 42]
info = subprocess.check_output([
@@ -41,6 +43,3 @@ def test_set_nodata(tmpdir):
assert re.search(pattern, info, re.DOTALL) is not None
pattern = b'Band 2.*?NoData Value=42'
assert re.search(pattern, info, re.DOTALL) is not None
-
-
-
diff --git a/tests/test_profile.py b/tests/test_profile.py
new file mode 100644
index 0000000..c277da5
--- /dev/null
+++ b/tests/test_profile.py
@@ -0,0 +1,79 @@
+import pytest
+
+import rasterio
+from rasterio.profiles import Profile, DefaultGTiffProfile
+from rasterio.profiles import default_gtiff_profile
+
+
+def test_base_profile():
+ assert Profile()()['driver'] is None
+
+
+def test_base_profile_kwarg():
+ assert Profile()(foo='bar')['foo'] == 'bar'
+
+
+def test_gtiff_profile_format():
+ assert DefaultGTiffProfile()()['driver'] == 'GTiff'
+
+
+def test_gtiff_profile_interleave():
+ assert DefaultGTiffProfile()()['interleave'] == 'band'
+
+
+def test_gtiff_profile_tiled():
+ assert DefaultGTiffProfile()()['tiled'] == True
+
+
+def test_gtiff_profile_blockxsize():
+ assert DefaultGTiffProfile()()['blockxsize'] == 256
+
+
+def test_gtiff_profile_blockysize():
+ assert DefaultGTiffProfile()()['blockysize'] == 256
+
+
+def test_gtiff_profile_compress():
+ assert DefaultGTiffProfile()()['compress'] == 'lzw'
+
+
+def test_gtiff_profile_nodata():
+ assert DefaultGTiffProfile()()['nodata'] == 0
+
+
+def test_gtiff_profile_dtype():
+ assert DefaultGTiffProfile()()['dtype'] == rasterio.uint8
+
+
+def test_gtiff_profile_other():
+ assert DefaultGTiffProfile()(count=3)['count'] == 3
+
+
+def test_gtiff_profile_dtype_override():
+ assert DefaultGTiffProfile()(dtype='uint16')['dtype'] == rasterio.uint16
+
+
+def test_gtiff_profile_protected_driver():
+ """Overriding the driver is not allowed."""
+ with pytest.raises(ValueError):
+ DefaultGTiffProfile()(driver='PNG')
+
+
+def test_open_with_profile(tmpdir):
+ tiffname = str(tmpdir.join('foo.tif'))
+ with rasterio.open(
+ tiffname,
+ 'w',
+ **default_gtiff_profile(
+ count=1, width=1, height=1)) as dst:
+ data = dst.read()
+ assert data.flatten().tolist() == [0]
+
+
+def test_profile_overlay():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ kwds = src.meta
+ kwds.update(**default_gtiff_profile())
+ assert kwds['tiled']
+ assert kwds['compress'] == 'lzw'
+ assert kwds['count'] == 3
diff --git a/tests/test_read.py b/tests/test_read.py
index 98cf22d..572c01d 100644
--- a/tests/test_read.py
+++ b/tests/test_read.py
@@ -92,7 +92,7 @@ class ReaderContextTest(unittest.TestCase):
def test_read_basic(self):
with rasterio.open('tests/data/shade.tif') as s:
- a = s.read() # Gray
+ a = s.read(masked=True) # Gray
self.assertEqual(a.ndim, 3)
self.assertEqual(a.shape, (1, 1024, 1024))
self.assertTrue(hasattr(a, 'mask'))
@@ -101,7 +101,7 @@ class ReaderContextTest(unittest.TestCase):
self.assertEqual(a.dtype, rasterio.ubyte)
self.assertEqual(a.sum((1, 2)).tolist(), [0])
with rasterio.open('tests/data/RGB.byte.tif') as s:
- a = s.read() # RGB
+ a = s.read(masked=True) # RGB
self.assertEqual(a.ndim, 3)
self.assertEqual(a.shape, (3, 718, 791))
self.assertTrue(hasattr(a, 'mask'))
@@ -113,10 +113,10 @@ class ReaderContextTest(unittest.TestCase):
self.assertEqual(list(set(s.nodatavals)), [0])
self.assertEqual(a.dtype, rasterio.ubyte)
with rasterio.open('tests/data/float.tif') as s:
- a = s.read() # floating point values
+ a = s.read(masked=True) # floating point values
self.assertEqual(a.ndim, 3)
self.assertEqual(a.shape, (1, 2, 3))
- self.assertFalse(hasattr(a, 'mask'))
+ self.assert_(hasattr(a, 'mask'))
self.assertEqual(list(set(s.nodatavals)), [None])
self.assertEqual(a.dtype, rasterio.float64)
@@ -161,7 +161,7 @@ class ReaderContextTest(unittest.TestCase):
# correct format
self.assertRaises(ValueError, s.read, window=(300, 320, 320, 330))
# window with 1 nodata on band 3
- a = s.read(window=((300, 320), (320, 330)))
+ a = s.read(window=((300, 320), (320, 330)), masked=True)
self.assertEqual(a.ndim, 3)
self.assertEqual(a.shape, (3, 20, 10))
self.assertTrue(hasattr(a, 'mask'))
@@ -171,7 +171,7 @@ class ReaderContextTest(unittest.TestCase):
'ec8fb3659f40c4a209027231bef12bdb',
'5a9c12aebc126ec6f27604babd67a4e2'])
# window without any missing data, but still is masked result
- a = s.read(window=((310, 330), (320, 330)))
+ a = s.read(window=((310, 330), (320, 330)), masked=True)
self.assertEqual(a.ndim, 3)
self.assertEqual(a.shape, (3, 20, 10))
self.assertTrue(hasattr(a, 'mask'))
@@ -206,7 +206,6 @@ class ReaderContextTest(unittest.TestCase):
# regular array, without mask
a = numpy.empty((3, 718, 791), numpy.ubyte)
b = s.read(out=a)
- self.assertEqual(id(a), id(b))
self.assertFalse(hasattr(a, 'mask'))
self.assertFalse(hasattr(b, 'mask'))
# with masked array
@@ -237,7 +236,7 @@ class ReaderContextTest(unittest.TestCase):
def test_read_nan_nodata(self):
with rasterio.open('tests/data/float_nan.tif') as s:
- a = s.read()
+ a = s.read(masked=True)
self.assertEqual(a.ndim, 3)
self.assertEqual(a.shape, (1, 2, 3))
self.assertTrue(hasattr(a, 'mask'))
@@ -249,7 +248,7 @@ class ReaderContextTest(unittest.TestCase):
self.assertFalse(hasattr(a, 'mask'))
self.assertTrue(numpy.isnan(a).any())
# Window does not contain a nodatavalue, result is still masked
- a = s.read(window=((0, 2), (0, 2)))
+ a = s.read(window=((0, 2), (0, 2)), masked=True)
self.assertEqual(a.ndim, 3)
self.assertEqual(a.shape, (1, 2, 2))
self.assertTrue(hasattr(a, 'mask'))
diff --git a/tests/test_read_boundless.py b/tests/test_read_boundless.py
index e6c2a1f..dc29c36 100644
--- a/tests/test_read_boundless.py
+++ b/tests/test_read_boundless.py
@@ -16,6 +16,7 @@ def test_read_boundless_natural_extent():
assert abs(data[0].mean() - src.read(1).mean()) < 0.0001
assert data.any()
+
def test_read_boundless_beyond():
with rasterio.open('tests/data/RGB.byte.tif') as src:
data = src.read(window=((-200, -100), (-200, -100)), boundless=True)
@@ -49,3 +50,22 @@ def test_read_boundless_resample():
assert data.shape == (3, 800, 800)
assert data.any()
assert data[0,798,798] == 13
+
+
+def test_read_boundless_masked_no_overlap():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ data = src.read(
+ window=((-200, -100), (-200, -100)), boundless=True, masked=True)
+ assert data.shape == (3, 100, 100)
+ assert data.mask.all()
+
+
+def test_read_boundless_masked_overlap():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ data = src.read(
+ window=((-200, 200), (-200, 200)), boundless=True, masked=True)
+ assert data.shape == (3, 400, 400)
+ assert data.mask.any()
+ assert not data.mask.all()
+ assert data.mask[0,399,399] == False
+ assert data.mask[0,0,0] == True
diff --git a/tests/test_read_resample.py b/tests/test_read_resample.py
index f3eb1bc..d21c0a4 100644
--- a/tests/test_read_resample.py
+++ b/tests/test_read_resample.py
@@ -30,4 +30,4 @@ def test_read_out_shape_resample_up():
out = numpy.zeros((7180, 7910), dtype=rasterio.ubyte)
data = s.read(1, out=out, masked=True)
assert data.shape == (7180, 7910)
- assert data.mean() == s.read(1).mean()
+ assert data.mean() == s.read(1, masked=True).mean()
diff --git a/tests/test_rio_calc.py b/tests/test_rio_calc.py
new file mode 100644
index 0000000..40f0bb4
--- /dev/null
+++ b/tests/test_rio_calc.py
@@ -0,0 +1,180 @@
+import sys
+import logging
+
+from click.testing import CliRunner
+
+import rasterio
+from rasterio.rio.calc import calc
+
+
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+
+def test_err(tmpdir):
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '($ 0.1 (read 1))', 'tests/data/shade.tif', outfile],
+ catch_exceptions=False)
+ assert result.exit_code == 1
+
+
+def test_multiband_calc(tmpdir):
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '(+ 125 (* 0.1 (read 1)))', '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'
+ data = src.read(masked=True)
+ assert data.min() == 125
+ assert data.data[0][0][0] == 255
+ assert data.mask[0][0][0]
+
+
+def test_singleband_calc_byindex(tmpdir):
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '(+ 125 (* 0.1 (read 1 1)))',
+ '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'
+ data = src.read(masked=True)
+ assert data.min() == 125
+
+
+def test_singleband_calc_byname(tmpdir):
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '(+ 125 (* 0.1 (take shade 1)))',
+ '--name', 'shade=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'
+ data = src.read(masked=True)
+ assert data.min() == 125
+
+
+def test_parts_calc(tmpdir):
+ # Producing an RGB output from the hill shade.
+ # Red band has bumped up values. Other bands are unchanged.
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '(asarray (+ (read 1 1) 125) (read 1 1) (read 1 1))',
+ 'tests/data/shade.tif',
+ outfile],
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ with rasterio.open(outfile) as src:
+ assert src.count == 3
+ assert src.meta['dtype'] == 'uint8'
+ data = src.read(masked=True)
+ assert data[0].min() == 125
+ assert data[1].min() == 0
+ assert data[2].min() == 0
+
+
+def test_parts_calc_2(tmpdir):
+ # Produce greyscale output from the RGB file.
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '(+ (+ (/ (read 1 1) 3.0) (/ (read 1 2) 3.0)) (/ (read 1 3) 3.0))',
+ 'tests/data/RGB.byte.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'
+ data = src.read(masked=True)
+ assert round(data.mean(), 1) == 60.3
+
+
+def test_copy_rgb(tmpdir):
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '(read 1)',
+ 'tests/data/RGB.byte.tif',
+ outfile],
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ with rasterio.open(outfile) as src:
+ assert src.count == 3
+ assert src.meta['dtype'] == 'uint8'
+ data = src.read(masked=True)
+ assert round(data.mean(), 1) == 60.6
+
+
+def test_fillnodata(tmpdir):
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '(asarray (fillnodata (read 1 1)) (fillnodata (read 1 2)) (fillnodata (read 1 3)))',
+ 'tests/data/RGB.byte.tif',
+ outfile],
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ with rasterio.open(outfile) as src:
+ assert src.count == 3
+ assert src.meta['dtype'] == 'uint8'
+ data = src.read(masked=True)
+ assert round(data.mean(), 1) == 60.6
+
+
+def test_fillnodata_map(tmpdir):
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '(asarray (map fillnodata (read 1)))',
+ 'tests/data/RGB.byte.tif',
+ outfile],
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ with rasterio.open(outfile) as src:
+ assert src.count == 3
+ assert src.meta['dtype'] == 'uint8'
+ data = src.read(masked=True)
+ assert round(data.mean(), 1) == 60.6
+
+
+def test_sieve_band(tmpdir):
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ '(sieve (band 1 1) 42)',
+ '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'
+
+
+def test_sieve_read(tmpdir):
+ outfile = str(tmpdir.join('out.tif'))
+ runner = CliRunner()
+ result = runner.invoke(calc, [
+ "(sieve (read 1 1 'uint8') 42)",
+ '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'
diff --git a/tests/test_rio_cli.py b/tests/test_rio_cli.py
new file mode 100644
index 0000000..875e3fa
--- /dev/null
+++ b/tests/test_rio_cli.py
@@ -0,0 +1,18 @@
+from rasterio.rio import cli
+
+
+def test_resolve_files_inout__output():
+ assert cli.resolve_inout(input='in', output='out') == ('out', ['in'])
+
+
+def test_resolve_files_inout__input():
+ assert cli.resolve_inout(input='in') == (None, ['in'])
+
+
+def test_resolve_files_inout__inout_files():
+ assert cli.resolve_inout(files=('a', 'b', 'c')) == ('c', ['a', 'b'])
+
+
+def test_resolve_files_inout__inout_files_output_o():
+ assert cli.resolve_inout(
+ files=('a', 'b', 'c'), output='out') == ('out', ['a', 'b', 'c'])
diff --git a/tests/test_rio_features.py b/tests/test_rio_features.py
index e524cd7..dd3da8e 100644
--- a/tests/test_rio_features.py
+++ b/tests/test_rio_features.py
@@ -2,15 +2,13 @@ import logging
import os
import re
import sys
-
-import click
-from click.testing import CliRunner
+import numpy
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": {
@@ -31,6 +29,24 @@ TEST_FEATURES = """{
"type": "Feature"
}"""
+TEST_MERC_FEATURES = """{
+ "geometry": {
+ "coordinates": [
+ [
+ [-11858134, 4808920],
+ [-11868134, 4804143],
+ [-11853357, 4804143],
+ [-11840000, 4812000],
+ [-11858134, 4808920]
+ ]
+ ],
+ "type": "Polygon"
+ },
+ "properties": {
+ "val": 10
+ },
+ "type": "Feature"
+}"""
# > rio shapes tests/data/shade.tif --mask --sampling 500 --projected --precision 0
TEST_MERC_FEATURECOLLECTION = """{
@@ -78,11 +94,165 @@ TEST_MERC_FEATURECOLLECTION = """{
}"""
-def test_err():
- runner = CliRunner()
+def test_mask(runner, tmpdir):
+ output = str(tmpdir.join('test.tif'))
+
+ with rasterio.open('tests/data/shade.tif') as src:
+ src_data = src.read(1, masked=True)
+
+ result = runner.invoke(
+ features.mask,
+ ['tests/data/shade.tif', output, '--geojson-mask', '-'],
+ input=TEST_MERC_FEATURES
+ )
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+
+ masked_count = 0
+ with rasterio.open(output) as out:
+ assert out.count == src.count
+ assert out.shape == src.shape
+ out_data = out.read(1, masked=True)
+
+ # Make sure that pixels with 0 value were converted to mask
+ masked_count = (src_data == 0).sum() - (out_data == 0).sum()
+ assert masked_count == 79743
+ assert out_data.mask.sum() - src_data.mask.sum() == masked_count
+
+ # Test using --all-touched option
+ result = runner.invoke(
+ features.mask,
+ [
+ 'tests/data/shade.tif', output,
+ '--all',
+ '--geojson-mask', '-'
+ ],
+ input=TEST_MERC_FEATURES
+ )
+ assert result.exit_code == 0
+ with rasterio.open(output) as out:
+ out_data = out.read(1, masked=True)
+
+ # Make sure that more pixels with 0 value were converted to mask
+ masked_count2 = (src_data == 0).sum() - (out_data == 0).sum()
+ assert masked_count2 > 0 and masked_count > masked_count2
+ assert out_data.mask.sum() - src_data.mask.sum() == masked_count2
+
+ # Test using --invert option
+ result = runner.invoke(
+ features.mask,
+ [
+ 'tests/data/shade.tif', output,
+ '--invert',
+ '--geojson-mask', '-'
+ ],
+ input=TEST_MERC_FEATURES
+ )
+ assert result.exit_code == 0
+ with rasterio.open(output) as out:
+ out_data = out.read(1, masked=True)
+ # Areas that were masked when not inverted should now be 0
+ assert (out_data == 0).sum() == masked_count
+
+ # Test with feature collection
result = runner.invoke(
- features.shapes, ['tests/data/shade.tif', '--bidx', '4'])
- assert result.exit_code == 1
+ features.mask,
+ ['tests/data/shade.tif', output, '--geojson-mask', '-'],
+ input=TEST_MERC_FEATURECOLLECTION
+ )
+ assert result.exit_code == 0
+
+ # Missing GeoJSON should make copy of input to output
+ output2 = str(tmpdir.join('test2.tif'))
+ result = runner.invoke(
+ features.mask,
+ ['tests/data/shade.tif', output2]
+ )
+ assert result.exit_code == 0
+ assert os.path.exists(output2)
+ with rasterio.open('tests/data/shade.tif') as src:
+ with rasterio.open(output2) as out:
+ src_data = src.read(1, masked=True)
+ out_data = out.read(1, masked=True)
+ assert numpy.array_equal(src_data, out_data)
+
+ # Invalid JSON should fail
+ result = runner.invoke(
+ features.mask,
+ ['tests/data/shade.tif', output, '--geojson-mask', '-'],
+ input='{bogus: value}'
+ )
+ assert result.exit_code == 2
+ assert 'GeoJSON could not be read' in result.output
+
+ result = runner.invoke(
+ features.mask,
+ ['tests/data/shade.tif', output, '--geojson-mask', '-'],
+ input='{"bogus": "value"}'
+ )
+ assert result.exit_code == 2
+ assert 'Invalid GeoJSON' in result.output
+
+
+def test_mask_crop(runner, tmpdir):
+ output = str(tmpdir.join('test.tif'))
+
+ with rasterio.open('tests/data/shade.tif') as src:
+
+ result = runner.invoke(
+ features.mask,
+ [
+ 'tests/data/shade.tif', output,
+ '--crop',
+ '--geojson-mask', '-'
+ ],
+ input=TEST_MERC_FEATURES
+ )
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+ 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
+
+ # Adding invert option after crop should be ignored
+ result = runner.invoke(
+ features.mask,
+ [
+ 'tests/data/shade.tif', output,
+ '--crop',
+ '--invert',
+ '--geojson-mask', '-'
+ ],
+ input=TEST_MERC_FEATURES
+ )
+ assert result.exit_code == 0
+ assert 'Invert option ignored' in result.output
+
+
+def test_mask_out_of_bounds(runner, tmpdir):
+ output = str(tmpdir.join('test.tif'))
+ # Crop with out of bounds raster should
+ result = runner.invoke(
+ features.mask,
+ ['tests/data/shade.tif', output, '--geojson-mask', '-'],
+ input=TEST_FEATURES
+ )
+ assert result.exit_code == 0
+ assert 'outside bounds' in result.output
+
+ # Crop with out of bounds raster should fail
+ result = runner.invoke(
+ features.mask,
+ [
+ 'tests/data/shade.tif', output,
+ '--crop',
+ '--geojson-mask', '-'
+ ],
+ input=TEST_FEATURES
+ )
+ assert result.exit_code == 2
+ assert 'not allowed' in result.output
def test_shapes(runner):
@@ -91,6 +261,11 @@ def test_shapes(runner):
assert result.output.count('"FeatureCollection"') == 1
assert result.output.count('"Feature"') == 232
+ # Invalid band index should fail
+ result = runner.invoke(
+ features.shapes, ['tests/data/shade.tif', '--bidx', '4'])
+ assert result.exit_code == 1
+
def test_shapes_sequence(runner):
result = runner.invoke(features.shapes, ['tests/data/shade.tif', '--sequence'])
@@ -122,7 +297,7 @@ 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
+ assert result.output.count('\n') == 70371
def test_shapes_compact(runner):
@@ -154,6 +329,23 @@ 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"') == 7
+
+
+def test_shapes_mask_decimated(runner):
+ result = runner.invoke(
+ features.shapes,
+ ['tests/data/RGB.byte.tif', '--mask', '--sampling', '10'])
+ assert result.exit_code == 0
+ assert result.output.count('"FeatureCollection"') == 1
+ assert result.output.count('"Feature"') == 1
+
+
+def test_shapes_band1_as_mask(runner):
+ result = runner.invoke(features.shapes,
+ ['tests/data/RGB.byte.tif', '--band', '--bidx', '1', '--as-mask'])
+ assert result.exit_code == 0
+ assert result.output.count('"FeatureCollection"') == 1
assert result.output.count('"Feature"') == 9
@@ -180,7 +372,7 @@ def test_rasterize_err(tmpdir, runner):
# Test invalid CRS value
result = runner.invoke(features.rasterize, [output,
'--res', 1,
- '--src_crs', 'BOGUS'],
+ '--src-crs', 'BOGUS'],
input=TEST_MERC_FEATURECOLLECTION)
assert result.exit_code == 2
@@ -197,7 +389,7 @@ def test_rasterize(tmpdir, runner):
assert out.count == 1
assert out.meta['width'] == 20
assert out.meta['height'] == 10
- data = out.read_band(1, masked=False)
+ data = out.read(1, masked=False)
assert (data == 0).sum() == 55
assert (data == 1).sum() == 145
@@ -214,7 +406,7 @@ def test_rasterize(tmpdir, runner):
assert out.count == 1
assert out.meta['width'] == 40
assert out.meta['height'] == 20
- data = out.read_band(1, masked=False)
+ data = out.read(1, masked=False)
assert (data == 0).sum() == 748
assert (data == 1).sum() == 52
@@ -228,10 +420,22 @@ def test_rasterize(tmpdir, runner):
assert out.count == 1
assert out.meta['width'] == 20
assert out.meta['height'] == 10
- data = out.read_band(1, masked=False)
+ data = out.read(1, masked=False)
assert (data == 0).sum() == 55
assert (data == 1).sum() == 145
+ # Test that src-crs is written into new output
+ output = str(tmpdir.join('test4.tif'))
+ result = runner.invoke(features.rasterize,
+ [output,
+ '--dimensions', 20, 10,
+ '--src-crs', 'EPSG:3857'
+ ],
+ input=TEST_MERC_FEATURECOLLECTION)
+ assert result.exit_code == 0
+ with rasterio.open(output) as out:
+ assert out.crs['init'].lower() == 'epsg:3857'
+
def test_rasterize_existing_output(tmpdir, runner):
output = str(tmpdir.join('test.tif'))
@@ -256,16 +460,24 @@ def test_rasterize_existing_output(tmpdir, runner):
"type": "Feature"
}"""
- result = runner.invoke(features.rasterize, [output, '--default_value', 2],
+ 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)
+ data = out.read(1, masked=False)
assert (data == 0).sum() == 55
assert (data == 1).sum() == 125
assert (data == 2).sum() == 20
+ # Confirm that a different src-crs is rejected, even if a geographic crs
+ result = runner.invoke(features.rasterize,
+ [output,
+ '--res', 0.5,
+ '--src-crs', 'EPSG:4269'
+ ], input=TEST_FEATURES)
+ assert result.exit_code == 2
+
def test_rasterize_like(tmpdir, runner):
output = str(tmpdir.join('test.tif'))
@@ -276,14 +488,23 @@ def test_rasterize_like(tmpdir, runner):
assert os.path.exists(output)
with rasterio.open(output) as out:
assert out.count == 1
- data = out.read_band(1, masked=False)
+ data = out.read(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)
+ [output, '--like', str(tmpdir.join('foo.tif'))], input=TEST_FEATURES)
+ assert result.exit_code == 2
+
+ # Test that src-crs different than --like raster crs breaks
+ output = str(tmpdir.join('test3.tif'))
+ result = runner.invoke(features.rasterize,
+ [output,
+ '--like', 'tests/data/shade.tif',
+ '--src-crs', 'EPSG:4326'],
+ input=TEST_FEATURES)
assert result.exit_code == 2
@@ -294,14 +515,14 @@ def test_rasterize_property_value(tmpdir, runner):
[output,
'--res', 1000,
'--property', 'val',
- '--src_crs', 'EPSG:3857'
+ '--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)
+ data = out.read(1, masked=False)
assert (data == 0).sum() == 50
assert (data == 2).sum() == 25
assert (data == 3).sum() == 25
@@ -315,10 +536,11 @@ def test_rasterize_property_value(tmpdir, runner):
assert os.path.exists(output)
with rasterio.open(output) as out:
assert out.count == 1
- data = out.read_band(1, masked=False)
+ data = out.read(1, masked=False)
assert (data == 0).sum() == 55
assert (data == 15).sum() == 145
+
def test_rasterize_out_of_bounds(tmpdir, runner):
output = str(tmpdir.join('test.tif'))
@@ -326,18 +548,14 @@ def test_rasterize_out_of_bounds(tmpdir, runner):
result = runner.invoke(features.rasterize,
[output, '--like', 'tests/data/shade.tif'],
input=TEST_FEATURES)
- assert result.exit_code == 2
-
- # Test out of bounds of existing output raster (first have to create one)
- result = runner.invoke(features.rasterize,
- [output,
- '--res', 1000,
- '--property', 'val',
- '--src_crs', 'EPSG:3857'
- ],
- input=TEST_MERC_FEATURECOLLECTION)
assert result.exit_code == 0
+ assert 'outside bounds' in result.output
assert os.path.exists(output)
+ with rasterio.open(output) as out:
+ data = out.read_band(1, masked=False)
+ assert data.sum() == 0
+ # Confirm that this does not fail when out of bounds for existing raster
result = runner.invoke(features.rasterize, [output], input=TEST_FEATURES)
- assert result.exit_code == 2
+ assert result.exit_code == 0
+ assert 'outside bounds' in result.output
diff --git a/tests/test_rio_info.py b/tests/test_rio_info.py
index 0dd5c65..c1a5d91 100644
--- a/tests/test_rio_info.py
+++ b/tests/test_rio_info.py
@@ -1,11 +1,116 @@
+import json
+import logging
+import sys
+
import click
from click.testing import CliRunner
-
import rasterio
from rasterio.rio import cli, info
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+
+def test_edit_nodata_err(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ result = runner.invoke(info.edit, [inputfile, '--nodata', '-1'])
+ assert result.exit_code == 2
+
+
+def test_edit_nodata(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ result = runner.invoke(info.edit, [inputfile, '--nodata', '255'])
+ assert result.exit_code == 0
+ with rasterio.open(inputfile) as src:
+ assert src.nodata == 255.0
+
+
+def test_edit_crs_err(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ result = runner.invoke(info.edit, [inputfile, '--crs', 'LOL:WUT'])
+ assert result.exit_code == 2
+
+
+def test_edit_crs_epsg(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ result = runner.invoke(info.edit, [inputfile, '--crs', 'EPSG:32618'])
+ assert result.exit_code == 0
+ with rasterio.open(inputfile) as src:
+ assert src.crs == {'init': 'epsg:32618'}
+
+
+def test_edit_crs_proj4(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ result = runner.invoke(info.edit, [inputfile, '--crs', '+init=epsg:32618'])
+ assert result.exit_code == 0
+ with rasterio.open(inputfile) as src:
+ assert src.crs == {'init': 'epsg:32618'}
+
+
+def test_edit_crs_obj(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ result = runner.invoke(
+ info.edit, [inputfile, '--crs', '{"init": "epsg:32618"}'])
+ assert result.exit_code == 0
+ with rasterio.open(inputfile) as src:
+ assert src.crs == {'init': 'epsg:32618'}
+
+
+def test_edit_transform_err_not_json(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ result = runner.invoke(info.edit, [inputfile, '--transform', 'LOL'])
+ assert result.exit_code == 2
+
+
+def test_edit_transform_err_bad_array(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ result = runner.invoke(info.edit, [inputfile, '--transform', '[1,2]'])
+ assert result.exit_code == 2
+
+
+def test_edit_transform_affine(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ input_t = '[300.038, 0.0, 101985.0, 0.0, -300.042, 2826915.0]'
+ result = runner.invoke(info.edit, [inputfile, '--transform', input_t])
+ assert result.exit_code == 0
+ with rasterio.open(inputfile) as src:
+ for a, b in zip(src.affine, json.loads(input_t)):
+ assert round(a, 6) == round(b, 6)
+
+
+def test_edit_transform_gdal(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ input_t = '[300.038, 0.0, 101985.0, 0.0, -300.042, 2826915.0]'
+ result = runner.invoke(info.edit, [
+ inputfile,
+ '--transform', '[101985.0, 300.038, 0.0, 2826915.0, 0.0, -300.042]'])
+ assert result.exit_code == 0
+ with rasterio.open(inputfile) as src:
+ for a, b in zip(src.affine, json.loads(input_t)):
+ assert round(a, 6) == round(b, 6)
+
+
+def test_edit_tags(data):
+ runner = CliRunner()
+ inputfile = str(data.join('RGB.byte.tif'))
+ result = runner.invoke(info.edit, [
+ inputfile, '--tag', 'lol=1', '--tag', 'wut=2'])
+ assert result.exit_code == 0
+ with rasterio.open(inputfile) as src:
+ assert src.tags()['lol'] == '1'
+ assert src.tags()['wut'] == '2'
+
def test_env():
runner = CliRunner()
result = runner.invoke(info.env, ['--formats'])
diff --git a/tests/test_rio_merge.py b/tests/test_rio_merge.py
index ddcb0a2..5fa5f91 100644
--- a/tests/test_rio_merge.py
+++ b/tests/test_rio_merge.py
@@ -8,6 +8,7 @@ from pytest import fixture
import rasterio
from rasterio.rio.merge import merge
+from rasterio.transform import Affine
logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
@@ -92,7 +93,7 @@ def test_merge_warn(test_data_dir_1):
inputs = [str(x) for x in test_data_dir_1.listdir()]
inputs.sort()
runner = CliRunner()
- result = runner.invoke(merge, inputs + [outputname] + ['-n', '-1'])
+ result = runner.invoke(merge, inputs + [outputname] + ['--nodata', '-1'])
assert result.exit_code == 0
assert "using the --nodata option for better results" in result.output
@@ -235,7 +236,7 @@ def test_merge_float(test_data_dir_float):
inputs = [str(x) for x in test_data_dir_float.listdir()]
inputs.sort()
runner = CliRunner()
- result = runner.invoke(merge, inputs + [outputname] + ['-n', '-1.5'])
+ result = runner.invoke(merge, inputs + [outputname] + ['--nodata', '-1.5'])
assert result.exit_code == 0
assert os.path.exists(outputname)
with rasterio.open(outputname) as out:
@@ -245,3 +246,86 @@ def test_merge_float(test_data_dir_float):
expected[0:6, 0:6] = 255
expected[4:8, 4:8] = 254
assert numpy.all(data == expected)
+
+
+# Test below comes from issue #288. There was an off-by-one error in
+# pasting image data into the canvas array.
+
+ at fixture(scope='function')
+def tiffs(tmpdir):
+
+ data = numpy.ones((1, 1, 1), 'uint8')
+
+ kwargs = {'count': '1',
+ 'driver': 'GTiff',
+ 'dtype': 'uint8',
+ 'height': 1,
+ 'width': 1}
+
+ kwargs['transform'] = Affine( 1, 0, 1,
+ 0,-1, 1)
+ with rasterio.open(str(tmpdir.join('a-sw.tif')), 'w', **kwargs) as r:
+ r.write(data * 40)
+
+ kwargs['transform'] = Affine( 1, 0, 2,
+ 0,-1, 2)
+ with rasterio.open(str(tmpdir.join('b-ct.tif')), 'w', **kwargs) as r:
+ r.write(data * 60)
+
+ kwargs['transform'] = Affine( 2, 0, 3,
+ 0,-2, 4)
+ with rasterio.open(str(tmpdir.join('c-ne.tif')), 'w', **kwargs) as r:
+ r.write(data * 90)
+
+ kwargs['transform'] = Affine( 2, 0, 2,
+ 0,-2, 4)
+ with rasterio.open(str(tmpdir.join('d-ne.tif')), 'w', **kwargs) as r:
+ r.write(data * 120)
+
+ return tmpdir
+
+
+def test_merge_tiny(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])
+ assert result.exit_code == 0
+
+ # Output should be
+ #
+ # [[ 0 120 120 90]
+ # [ 0 120 120 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][2][1] == 60
+ assert data[0][3][0] == 40
+
+
+def test_merge_tiny_output_opt(tiffs):
+ outputname = str(tiffs.join('merged.tif'))
+ inputs = [str(x) for x in tiffs.listdir()]
+ inputs.sort()
+ runner = CliRunner()
+ result = runner.invoke(merge, inputs + ['-o', outputname])
+ assert result.exit_code == 0
+
+ # Output should be
+ #
+ # [[ 0 120 120 90]
+ # [ 0 120 120 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][2][1] == 60
+ assert data[0][3][0] == 40
diff --git a/tests/test_rio_rio.py b/tests/test_rio_rio.py
index 06a7774..b8e7ef3 100644
--- a/tests/test_rio_rio.py
+++ b/tests/test_rio_rio.py
@@ -137,7 +137,7 @@ def test_transform_point():
runner = CliRunner()
result = runner.invoke(
rio.transform,
- ['--dst_crs', 'EPSG:32618', '--precision', '2'],
+ ['--dst-crs', 'EPSG:32618', '--precision', '2'],
"[-78.0, 23.0]", catch_exceptions=False)
assert result.exit_code == 0
assert result.output.strip() == '[192457.13, 2546667.68]'
@@ -147,7 +147,7 @@ def test_transform_point_dst_file():
runner = CliRunner()
result = runner.invoke(
rio.transform,
- ['--dst_crs', 'tests/data/RGB.byte.tif', '--precision', '2'],
+ ['--dst-crs', 'tests/data/RGB.byte.tif', '--precision', '2'],
"[-78.0, 23.0]")
assert result.exit_code == 0
assert result.output.strip() == '[192457.13, 2546667.68]'
@@ -157,7 +157,7 @@ def test_transform_point_src_file():
runner = CliRunner()
result = runner.invoke(
rio.transform,
- ['--src_crs', 'tests/data/RGB.byte.tif', '--precision', '2'],
+ ['--src-crs', 'tests/data/RGB.byte.tif', '--precision', '2'],
"[192457.13, 2546667.68]")
assert result.exit_code == 0
assert result.output.strip() == '[-78.0, 23.0]'
@@ -167,7 +167,7 @@ def test_transform_point_2():
runner = CliRunner()
result = runner.invoke(
rio.transform,
- ['[-78.0, 23.0]', '--dst_crs', 'EPSG:32618', '--precision', '2'])
+ ['[-78.0, 23.0]', '--dst-crs', 'EPSG:32618', '--precision', '2'])
assert result.exit_code == 0
assert result.output.strip() == '[192457.13, 2546667.68]'
@@ -176,7 +176,7 @@ def test_transform_point_multi():
runner = CliRunner()
result = runner.invoke(
rio.transform,
- ['--dst_crs', 'EPSG:32618', '--precision', '2'],
+ ['--dst-crs', 'EPSG:32618', '--precision', '2'],
"[-78.0, 23.0]\n[-78.0, 23.0]", catch_exceptions=False)
assert result.exit_code == 0
assert result.output.strip() == '[192457.13, 2546667.68]\n[192457.13, 2546667.68]'
diff --git a/tests/test_rio_sample.py b/tests/test_rio_sample.py
index 2c6a329..deae77a 100644
--- a/tests/test_rio_sample.py
+++ b/tests/test_rio_sample.py
@@ -28,7 +28,7 @@ def test_sample_stdin():
"[220650.0, 2719200.0]\n[220650.0, 2719200.0]",
catch_exceptions=False)
assert result.exit_code == 0
- assert result.output.strip() == '[28, 29, 27]\n[28, 29, 27]'
+ assert result.output.strip() == '[18, 25, 14]\n[18, 25, 14]'
def test_sample_arg():
@@ -38,7 +38,7 @@ def test_sample_arg():
['tests/data/RGB.byte.tif', "[220650.0, 2719200.0]"],
catch_exceptions=False)
assert result.exit_code == 0
- assert result.output.strip() == '[28, 29, 27]'
+ assert result.output.strip() == '[18, 25, 14]'
def test_sample_bidx():
@@ -48,7 +48,7 @@ def test_sample_bidx():
['tests/data/RGB.byte.tif', '--bidx', '1,2', "[220650.0, 2719200.0]"],
catch_exceptions=False)
assert result.exit_code == 0
- assert result.output.strip() == '[28, 29]'
+ assert result.output.strip() == '[18, 25]'
def test_sample_bidx2():
@@ -58,7 +58,7 @@ def test_sample_bidx2():
['tests/data/RGB.byte.tif', '--bidx', '1..2', "[220650.0, 2719200.0]"],
catch_exceptions=False)
assert result.exit_code == 0
- assert result.output.strip() == '[28, 29]'
+ assert result.output.strip() == '[18, 25]'
def test_sample_bidx3():
@@ -68,7 +68,7 @@ def test_sample_bidx3():
['tests/data/RGB.byte.tif', '--bidx', '..2', "[220650.0, 2719200.0]"],
catch_exceptions=False)
assert result.exit_code == 0
- assert result.output.strip() == '[28, 29]'
+ assert result.output.strip() == '[18, 25]'
def test_sample_bidx4():
@@ -78,4 +78,4 @@ def test_sample_bidx4():
['tests/data/RGB.byte.tif', '--bidx', '3', "[220650.0, 2719200.0]"],
catch_exceptions=False)
assert result.exit_code == 0
- assert result.output.strip() == '[27]'
+ assert result.output.strip() == '[14]'
diff --git a/tests/test_sampling.py b/tests/test_sampling.py
index c4030cc..02eef7e 100644
--- a/tests/test_sampling.py
+++ b/tests/test_sampling.py
@@ -4,7 +4,7 @@ import rasterio
def test_sampling():
with rasterio.open('tests/data/RGB.byte.tif') as src:
data = next(src.sample([(220650.0, 2719200.0)]))
- assert list(data) == [28, 29, 27]
+ assert list(data) == [18, 25, 14]
def test_sampling_beyond_bounds():
with rasterio.open('tests/data/RGB.byte.tif') as src:
@@ -14,4 +14,4 @@ def test_sampling_beyond_bounds():
def test_sampling_indexes():
with rasterio.open('tests/data/RGB.byte.tif') as src:
data = next(src.sample([(220650.0, 2719200.0)], indexes=[2]))
- assert list(data) == [29]
+ assert list(data) == [25]
diff --git a/tests/test_tool.py b/tests/test_tool.py
index 5fa4b03..f579c92 100644
--- a/tests/test_tool.py
+++ b/tests/test_tool.py
@@ -13,11 +13,11 @@ def test_stats():
with rasterio.drivers():
with rasterio.open('tests/data/RGB.byte.tif') as src:
results = stats((src, 1))
- assert results[0] == 1
+ assert results[0] == 0
assert results[1] == 255
- assert np.isclose(results[2], 44.4344)
+ assert np.isclose(results[2], 29.9477)
- results2 = stats(src.read_band(1))
+ results2 = stats(src.read(1))
assert np.allclose(np.array(results), np.array(results2))
@@ -39,6 +39,6 @@ def test_show():
pass
try:
- show(src.read_band(1))
+ show(src.read(1))
except ImportError:
- pass
\ No newline at end of file
+ pass
diff --git a/tests/test_transform.py b/tests/test_transform.py
index acd0d54..cf0607e 100644
--- a/tests/test_transform.py
+++ b/tests/test_transform.py
@@ -1,12 +1,5 @@
import rasterio
-
-def test_window():
- with rasterio.open('tests/data/RGB.byte.tif') as src:
- left, bottom, right, top = src.bounds
- assert src.window(left, bottom, right, top) == ((0, src.height),
- (0, src.width))
- assert src.window(left, top-src.res[1], left+src.res[0], top) == (
- (0, 1), (0, 1))
+from rasterio import transform
def test_window_transform():
@@ -21,3 +14,63 @@ def test_window_transform():
((-1, None), (-1, None))).c == src.bounds.left - src.res[0]
assert src.window_transform(
((-1, None), (-1, None))).f == src.bounds.top + src.res[1]
+
+
+def test_from_origin():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ w, n = src.ul(0, 0)
+ xs, ys = src.res
+ tr = transform.from_origin(w, n, xs, ys)
+ assert [round(v, 7) for v in tr] == [round(v, 7) for v in src.affine]
+
+
+def test_from_bounds():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ w, s, e, n = src.bounds
+ tr = transform.from_bounds(w, s, e, n, src.width, src.height)
+ assert [round(v, 7) for v in tr] == [round(v, 7) for v in src.affine]
+
+
+def test_window_bounds():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+
+ rows = src.height
+ cols = src.width
+
+ # Test window for entire DS and each window in the DS
+ assert src.window_bounds(((0, rows), (0, cols))) == src.bounds
+ for _, window in src.block_windows():
+ ds_x_min, ds_y_min, ds_x_max, ds_y_max = src.bounds
+ w_x_min, w_y_min, w_x_max, w_y_max = src.window_bounds(window)
+ assert ds_x_min <= w_x_min <= w_x_max <= ds_x_max
+ assert ds_y_min <= w_y_min <= w_y_max <= ds_y_max
+
+ # Test a small window in each corner, both in and slightly out of bounds
+ p = 10
+ for window in (
+ # In bounds (UL, UR, LL, LR)
+ ((0, p), (0, p)),
+ ((0, p), (cols - p, p)),
+ ((rows - p, p), (0, p)),
+ ((rows - p, p), (cols - p, p)),
+
+ # Out of bounds (UL, UR, LL, LR)
+ ((-1, p), (-1, p)),
+ ((-1, p), (cols - p, p + 1)),
+ ((rows - p, p + 1), (-1, p)),
+ ((rows - p, p + 1), (cols - p, p + 1))):
+
+ # Alternate formula
+
+ ((row_min, row_max), (col_min, col_max)) = window
+ win_aff = src.window_transform(window)
+
+ x_min, y_max = win_aff.c, win_aff.f
+ x_max = win_aff.c + (src.res[0] * (col_max - col_min))
+ y_min = win_aff.f - (src.res[1] * (row_max - row_min))
+
+ expected = (x_min, y_min, x_max, y_max)
+ actual = src.window_bounds(window)
+
+ for e, a in zip(expected, actual):
+ assert round(e, 7) == round(a, 7)
diff --git a/tests/test_update.py b/tests/test_update.py
index 8d84748..a097a54 100644
--- a/tests/test_update.py
+++ b/tests/test_update.py
@@ -9,30 +9,30 @@ import pytest
import rasterio
-def test_update_tags(tmpdir):
- tiffname = str(tmpdir.join('foo.tif'))
- shutil.copy('tests/data/RGB.byte.tif', tiffname)
- with rasterio.open(tiffname, 'r+') as f:
- f.update_tags(a='1', b='2')
- f.update_tags(1, c=3)
- with pytest.raises(ValueError):
- f.update_tags(4, d=4)
- assert f.tags() == {'AREA_OR_POINT': 'Area', 'a': '1', 'b': '2'}
- assert ('c', '3') in f.tags(1).items()
- info = subprocess.check_output(["gdalinfo", tiffname]).decode('utf-8')
- assert re.search("Metadata:\W+a=1\W+AREA_OR_POINT=Area\W+b=2", info)
-
-def test_update_band(tmpdir):
- tiffname = str(tmpdir.join('foo.tif'))
- shutil.copy('tests/data/RGB.byte.tif', tiffname)
+def test_update_tags(data):
+ tiffname = str(data.join('RGB.byte.tif'))
+ with rasterio.drivers():
+ with rasterio.open(tiffname, 'r+') as f:
+ f.update_tags(a='1', b='2')
+ f.update_tags(1, c=3)
+ with pytest.raises(ValueError):
+ f.update_tags(4, d=4)
+ assert f.tags() == {'AREA_OR_POINT': 'Area', 'a': '1', 'b': '2'}
+ assert ('c', '3') in f.tags(1).items()
+ info = subprocess.check_output(["gdalinfo", tiffname]).decode('utf-8')
+ assert re.search("Metadata:\W+a=1\W+AREA_OR_POINT=Area\W+b=2", info)
+
+
+def test_update_band(data):
+ tiffname = str(data.join('RGB.byte.tif'))
with rasterio.open(tiffname, 'r+') as f:
f.write_band(1, numpy.zeros(f.shape, dtype=f.dtypes[0]))
with rasterio.open(tiffname) as f:
assert not f.read_band(1).any()
-def test_update_spatial(tmpdir):
- tiffname = str(tmpdir.join('foo.tif'))
- shutil.copy('tests/data/RGB.byte.tif', tiffname)
+
+def test_update_spatial(data):
+ tiffname = str(data.join('RGB.byte.tif'))
with rasterio.open(tiffname, 'r+') as f:
f.transform = affine.Affine.from_gdal(1.0, 1.0, 0.0, 0.0, 0.0, -1.0)
f.crs = {'init': 'epsg:4326'}
@@ -41,9 +41,9 @@ def test_update_spatial(tmpdir):
assert list(f.affine.to_gdal()) == [1.0, 1.0, 0.0, 0.0, 0.0, -1.0]
assert f.crs == {'init': 'epsg:4326'}
-def test_update_spatial_epsg(tmpdir):
- tiffname = str(tmpdir.join('foo.tif'))
- shutil.copy('tests/data/RGB.byte.tif', tiffname)
+
+def test_update_spatial_epsg(data):
+ tiffname = str(data.join('RGB.byte.tif'))
with rasterio.open(tiffname, 'r+') as f:
f.transform = affine.Affine.from_gdal(1.0, 1.0, 0.0, 0.0, 0.0, -1.0)
f.crs = 'EPSG:4326'
@@ -52,10 +52,38 @@ def test_update_spatial_epsg(tmpdir):
assert list(f.affine.to_gdal()) == [1.0, 1.0, 0.0, 0.0, 0.0, -1.0]
assert f.crs == {'init': 'epsg:4326'}
-def test_update_nodatavals(tmpdir):
- tiffname = str(tmpdir.join('foo.tif'))
- shutil.copy('tests/data/RGB.byte.tif', tiffname)
+
+def test_update_nodatavals(data):
+ tiffname = str(data.join('RGB.byte.tif'))
+ with rasterio.open(tiffname, 'r+') as f:
+ f.nodata = 255
+ with rasterio.open(tiffname) as f:
+ assert f.nodatavals == [255, 255, 255]
+
+
+def test_update_nodatavals_error(data):
+ """GDAL doesn't support un-setting nodata values."""
+ tiffname = str(data.join('RGB.byte.tif'))
+ with rasterio.open(tiffname, 'r+') as f:
+ try:
+ f.nodata = None
+ except TypeError:
+ pass
+
+
+def test_update_mask_true(data):
+ """Provide an option to set a uniformly valid mask."""
+ tiffname = str(data.join('RGB.byte.tif'))
+ with rasterio.open(tiffname, 'r+') as f:
+ f.write_mask(True)
+ with rasterio.open(tiffname) as f:
+ assert f.read_masks().all()
+
+
+def test_update_mask_false(data):
+ """Provide an option to set a uniformly invalid mask."""
+ tiffname = str(data.join('RGB.byte.tif'))
with rasterio.open(tiffname, 'r+') as f:
- f.nodatavals = [-1, -1, -1]
+ f.write_mask(False)
with rasterio.open(tiffname) as f:
- assert f.nodatavals == [-1, -1, -1]
+ assert not f.read_masks().any()
diff --git a/tests/test_warp.py b/tests/test_warp.py
index fcf33be..7ba8cb0 100644
--- a/tests/test_warp.py
+++ b/tests/test_warp.py
@@ -1,18 +1,53 @@
import logging
-import subprocess
import sys
-
-import affine
+import pytest
+from affine import Affine
import numpy
import rasterio
-from rasterio.warp import reproject, RESAMPLING, transform_geom, transform
+from rasterio.warp import (
+ reproject, RESAMPLING, transform_geom, transform, transform_bounds,
+ calculate_default_transform)
logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+DST_TRANSFORM = Affine.from_gdal(-8789636.708, 300.0, 0.0, 2943560.235, 0.0, -300.0)
+
+
+class ReprojectParams(object):
+ """ Class to assist testing reprojection by encapsulating parameters """
+ def __init__(self, left, bottom, right, top, width, height, src_crs,
+ dst_crs):
+ self.width = width
+ self.height = height
+ src_res = float(right - left) / float(width)
+ self.src_transform = Affine(src_res, 0, left, 0, -src_res, top)
+ self.src_crs = src_crs
+ self.dst_crs = dst_crs
+
+ with rasterio.drivers():
+ dt, dw, dh = calculate_default_transform(
+ left, bottom, right, top, width, height, src_crs, dst_crs)
+ self.dst_transform = dt
+ self.dst_width = dw
+ self.dst_height = dh
+
+
+def default_reproject_params():
+ return ReprojectParams(
+ left=-120,
+ bottom=30,
+ right=-80,
+ top=70,
+ width=80,
+ height=80,
+ src_crs={'init': 'EPSG:4326'},
+ dst_crs={'init': 'EPSG:32610'})
+
+
def test_transform():
"""2D and 3D"""
WGS84_crs = {'init': 'EPSG:4326'}
@@ -28,13 +63,133 @@ def test_transform():
assert numpy.allclose(numpy.array(UTM33_result), numpy.array(UTM33_points))
-def test_reproject():
- """Ndarry to ndarray"""
+def test_transform_bounds():
+ with rasterio.drivers():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ l, b, r, t = src.bounds
+ assert numpy.allclose(
+ transform_bounds(l, b, r, t, src.crs, {'init': 'EPSG:4326'}),
+ (
+ -78.95864996545055, 23.564991210854686,
+ -76.57492370013823, 25.550873767433984
+ )
+ )
+
+
+def test_transform_bounds_densify():
+ # This transform is non-linear along the edges, so densification produces
+ # a different result than otherwise
+ src_crs = {'init': 'EPSG:4326'}
+ dst_crs = {'init': 'EPSG:32610'}
+ assert numpy.allclose(
+ transform_bounds(
+ -120, 40, -80, 64,
+ src_crs,
+ dst_crs,
+ densify_pts=0
+ ),
+ (
+ 646695.227266598, 4432069.056898901,
+ 4201818.984205882, 7807592.187464975
+ )
+ )
+
+ assert numpy.allclose(
+ transform_bounds(
+ -120, 40, -80, 64,
+ src_crs,
+ dst_crs,
+ densify_pts=100
+ ),
+ (
+ 646695.2272665979, 4432069.056898901,
+ 4201818.984205882, 7807592.187464977
+ )
+ )
+
+
+def test_transform_bounds_no_change():
+ """ Make sure that going from and to the same crs causes no change """
+ with rasterio.drivers():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ l, b, r, t = src.bounds
+ assert numpy.allclose(
+ transform_bounds(l, b, r, t, src.crs, src.crs),
+ src.bounds
+ )
+
+
+def test_transform_bounds_densify_out_of_bounds():
+ with pytest.raises(ValueError):
+ transform_bounds(
+ -120, 40, -80, 64,
+ {'init': 'EPSG:4326'},
+ {'init': 'EPSG:32610'},
+ densify_pts=-10
+ )
+
+
+def test_calculate_default_transform():
+ target_transform = Affine(
+ 0.0028956983577810586, 0.0, -78.95864996545055,
+ 0.0, -0.0028956983577810586, 25.550873767433984
+ )
+ 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(
+ l, b, r, t, src.width, src.height, src.crs, wgs84_crs)
+
+ assert dst_transform.almost_equals(target_transform)
+ assert width == 824
+ assert height == 686
+
+
+def test_calculate_default_transform_single_resolution():
+ with rasterio.drivers():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ l, b, r, t = src.bounds
+ target_resolution = 0.1
+ target_transform = Affine(
+ target_resolution, 0.0, -78.95864996545055,
+ 0.0, -target_resolution, 25.550873767433984
+ )
+ dst_transform, width, height = calculate_default_transform(
+ l, b, r, t, src.width, src.height, src.crs,
+ {'init': 'EPSG:4326'}, resolution=target_resolution
+ )
+
+ assert dst_transform.almost_equals(target_transform)
+ assert width == 24
+ assert height == 20
+
+
+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,
+ 0.0, -target_resolution[1], 25.550873767433984
+ )
+
+ dst_transform, width, height = calculate_default_transform(
+ l, b, r, t, src.width, src.height, src.crs,
+ {'init': 'EPSG:4326'}, resolution=target_resolution
+ )
+
+ assert dst_transform.almost_equals(target_transform)
+ assert width == 12
+ assert height == 20
+
+
+def test_reproject_ndarray():
with rasterio.drivers():
with rasterio.open('tests/data/RGB.byte.tif') as src:
source = src.read_band(1)
- dst_transform = affine.Affine.from_gdal(
- -8789636.708, 300.0, 0.0, 2943560.235, 0.0, -300.0)
+
dst_crs = dict(
proj='merc',
a=6378137,
@@ -48,23 +203,169 @@ def test_reproject():
nadgrids='@null',
wktext=True,
no_defs=True)
- destin = numpy.empty(src.shape, dtype=numpy.uint8)
+ out = numpy.empty(src.shape, dtype=numpy.uint8)
reproject(
source,
- destin,
+ out,
src_transform=src.transform,
src_crs=src.crs,
- dst_transform=dst_transform,
+ dst_transform=DST_TRANSFORM,
dst_crs=dst_crs,
resampling=RESAMPLING.nearest)
- assert destin.any()
- try:
- import matplotlib.pyplot as plt
- plt.imshow(destin)
- plt.gray()
- plt.savefig('test_reproject.png')
- except:
- pass
+ assert (out > 0).sum() == 438146
+
+
+def test_reproject_epsg():
+ with rasterio.drivers():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ source = src.read_band(1)
+
+ dst_crs = {'init': 'EPSG:3857'}
+ out = numpy.empty(src.shape, dtype=numpy.uint8)
+ reproject(
+ source,
+ out,
+ src_transform=src.transform,
+ src_crs=src.crs,
+ dst_transform=DST_TRANSFORM,
+ dst_crs=dst_crs,
+ resampling=RESAMPLING.nearest)
+ assert (out > 0).sum() == 438146
+
+
+def test_reproject_out_of_bounds():
+ # using EPSG code not appropriate for the transform should return blank image
+ with rasterio.drivers():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ source = src.read_band(1)
+
+ dst_crs = {'init': 'EPSG:32619'}
+ out = numpy.empty(src.shape, dtype=numpy.uint8)
+ reproject(
+ source,
+ out,
+ src_transform=src.transform,
+ src_crs=src.crs,
+ dst_transform=DST_TRANSFORM,
+ dst_crs=dst_crs,
+ resampling=RESAMPLING.nearest)
+ assert not out.any()
+
+
+def test_reproject_nodata():
+ params = default_reproject_params()
+ nodata = 215
+
+ with rasterio.drivers():
+ source = numpy.ones((params.width, params.height), dtype=numpy.uint8)
+ out = numpy.zeros((params.dst_width, params.dst_height),
+ dtype=source.dtype)
+ out.fill(120) # Fill with arbitrary value
+
+ reproject(
+ source,
+ out,
+ src_transform=params.src_transform,
+ src_crs=params.src_crs,
+ src_nodata=nodata,
+ dst_transform=params.dst_transform,
+ dst_crs=params.dst_crs,
+ dst_nodata=nodata
+ )
+
+ assert (out == 1).sum() == 4461
+ assert (out == nodata).sum() == (params.dst_width *
+ params.dst_height - 4461)
+
+
+def test_reproject_dst_nodata_default():
+ """
+ If nodata is not provided, destination will be filled with 0
+ instead of nodata
+ """
+
+ params = default_reproject_params()
+
+ with rasterio.drivers():
+ source = numpy.ones((params.width, params.height), dtype=numpy.uint8)
+ out = numpy.zeros((params.dst_width, params.dst_height),
+ dtype=source.dtype)
+ out.fill(120) # Fill with arbitrary value
+
+ reproject(
+ source,
+ out,
+ src_transform=params.src_transform,
+ src_crs=params.src_crs,
+ dst_transform=params.dst_transform,
+ dst_crs=params.dst_crs
+ )
+
+ assert (out == 1).sum() == 4461
+ assert (out == 0).sum() == (params.dst_width *
+ params.dst_height - 4461)
+
+
+def test_reproject_invalid_dst_nodata():
+ """ dst_nodata must be in value range of data type """
+ params = default_reproject_params()
+
+ with rasterio.drivers():
+ source = numpy.ones((params.width, params.height), dtype=numpy.uint8)
+ out = source.copy()
+
+ with pytest.raises(ValueError):
+ reproject(
+ source,
+ out,
+ src_transform=params.src_transform,
+ src_crs=params.src_crs,
+ src_nodata=0,
+ dst_transform=params.dst_transform,
+ dst_crs=params.dst_crs,
+ dst_nodata=999999999
+ )
+
+
+def test_reproject_missing_src_nodata():
+ """ src_nodata is required if dst_nodata is not None """
+ params = default_reproject_params()
+
+ with rasterio.drivers():
+ source = numpy.ones((params.width, params.height), dtype=numpy.uint8)
+ out = source.copy()
+
+ with pytest.raises(ValueError):
+ reproject(
+ source,
+ out,
+ src_transform=params.src_transform,
+ src_crs=params.src_crs,
+ dst_transform=params.dst_transform,
+ dst_crs=params.dst_crs,
+ dst_nodata=215
+ )
+
+
+def test_reproject_invalid_src_nodata():
+ """ src_nodata must be in range for data type """
+ params = default_reproject_params()
+
+ with rasterio.drivers():
+ source = numpy.ones((params.width, params.height), dtype=numpy.uint8)
+ out = source.copy()
+
+ with pytest.raises(ValueError):
+ reproject(
+ source,
+ out,
+ src_transform=params.src_transform,
+ src_crs=params.src_crs,
+ src_nodata=999999999,
+ dst_transform=params.dst_transform,
+ dst_crs=params.dst_crs,
+ dst_nodata=215
+ )
def test_reproject_multi():
@@ -72,8 +373,6 @@ def test_reproject_multi():
with rasterio.drivers():
with rasterio.open('tests/data/RGB.byte.tif') as src:
source = src.read()
- dst_transform = affine.Affine.from_gdal(
- -8789636.708, 300.0, 0.0, 2943560.235, 0.0, -300.0)
dst_crs = dict(
proj='merc',
a=6378137,
@@ -93,24 +392,15 @@ def test_reproject_multi():
destin,
src_transform=src.transform,
src_crs=src.crs,
- dst_transform=dst_transform,
+ dst_transform=DST_TRANSFORM,
dst_crs=dst_crs,
resampling=RESAMPLING.nearest)
assert destin.any()
- try:
- import matplotlib.pyplot as plt
- plt.imshow(destin)
- plt.gray()
- plt.savefig('test_reproject.png')
- except:
- pass
def test_warp_from_file():
"""File to ndarray"""
with rasterio.open('tests/data/RGB.byte.tif') as src:
- dst_transform = affine.Affine.from_gdal(
- -8789636.708, 300.0, 0.0, 2943560.235, 0.0, -300.0)
dst_crs = dict(
proj='merc',
a=6378137,
@@ -128,24 +418,15 @@ def test_warp_from_file():
reproject(
rasterio.band(src, 1),
destin,
- dst_transform=dst_transform,
+ dst_transform=DST_TRANSFORM,
dst_crs=dst_crs)
assert destin.any()
- try:
- import matplotlib.pyplot as plt
- plt.imshow(destin)
- plt.gray()
- plt.savefig('test_warp_from_filereproject.png')
- except:
- pass
def test_warp_from_to_file(tmpdir):
"""File to file"""
tiffname = str(tmpdir.join('foo.tif'))
with rasterio.open('tests/data/RGB.byte.tif') as src:
- dst_transform = affine.Affine.from_gdal(
- -8789636.708, 300.0, 0.0, 2943560.235, 0.0, -300.0)
dst_crs = dict(
proj='merc',
a=6378137,
@@ -161,20 +442,17 @@ def test_warp_from_to_file(tmpdir):
no_defs=True)
kwargs = src.meta.copy()
kwargs.update(
- transform=dst_transform,
+ transform=DST_TRANSFORM,
crs=dst_crs)
with rasterio.open(tiffname, 'w', **kwargs) as dst:
for i in (1, 2, 3):
reproject(rasterio.band(src, i), rasterio.band(dst, i))
- # subprocess.call(['open', tiffname])
def test_warp_from_to_file_multi(tmpdir):
"""File to file"""
tiffname = str(tmpdir.join('foo.tif'))
with rasterio.open('tests/data/RGB.byte.tif') as src:
- dst_transform = affine.Affine.from_gdal(
- -8789636.708, 300.0, 0.0, 2943560.235, 0.0, -300.0)
dst_crs = dict(
proj='merc',
a=6378137,
@@ -190,7 +468,7 @@ def test_warp_from_to_file_multi(tmpdir):
no_defs=True)
kwargs = src.meta.copy()
kwargs.update(
- transform=dst_transform,
+ transform=DST_TRANSFORM,
crs=dst_crs)
with rasterio.open(tiffname, 'w', **kwargs) as dst:
for i in (1, 2, 3):
@@ -198,7 +476,6 @@ def test_warp_from_to_file_multi(tmpdir):
rasterio.band(src, i),
rasterio.band(dst, i),
num_threads=2)
- # subprocess.call(['open', tiffname])
def test_transform_geom():
diff --git a/tests/test_write.py b/tests/test_write.py
index 121e478..25f5178 100644
--- a/tests/test_write.py
+++ b/tests/test_write.py
@@ -1,10 +1,14 @@
import logging
+import re
import subprocess
import sys
-import re
+
import numpy
+import pytest
+
import rasterio
+
logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
@@ -206,6 +210,18 @@ def test_write_nodata(tmpdir):
info = subprocess.check_output(["gdalinfo", "-stats", name]).decode('utf-8')
assert "NoData Value=0" in info
+
+def test_guard_nodata(tmpdir):
+ name = str(tmpdir.join("test_guard_nodata.tif"))
+ a = numpy.ones((100, 100), dtype=rasterio.ubyte) * 127
+ with pytest.raises(ValueError):
+ with rasterio.open(
+ name, 'w',
+ driver='GTiff', width=100, height=100, count=2,
+ dtype=a.dtype, nodata=-1) as s:
+ pass
+
+
def test_write_lzw(tmpdir):
name = str(tmpdir.join("test_write_lzw.tif"))
a = numpy.ones((100, 100), dtype=rasterio.ubyte) * 127
--
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