[rasterio] 01/06: Imported Upstream version 0.32.0
Sebastiaan Couwenberg
sebastic at moszumanska.debian.org
Thu Mar 24 23:45:54 UTC 2016
This is an automated email from the git hooks/post-receive script.
sebastic pushed a commit to branch master
in repository rasterio.
commit e86356d8e7ba07f9d95aedc0c56b105ab65b91b1
Author: Bas Couwenberg <sebastic at xs4all.nl>
Date: Fri Mar 25 00:24:35 2016 +0100
Imported Upstream version 0.32.0
---
AUTHORS.txt | 37 ++--
CHANGES.txt | 37 +++-
README.rst | 2 +-
docs/cli.rst | 4 +-
docs/colormaps.rst | 2 +-
docs/concurrency.rst | 2 -
docs/datasets.rst | 8 +-
docs/features.rst | 6 +-
docs/georeferencing.rst | 4 +-
docs/masks.rst | 8 +-
docs/tags.rst | 2 +-
docs/windowed-rw.rst | 12 +-
examples/concurrent-cpu-bound.py | 2 -
examples/decimate.py | 2 +-
examples/rasterio_polygonize.py | 2 +-
examples/sieve.py | 2 +-
rasterio/__init__.py | 4 +-
rasterio/_base.pyx | 19 +-
rasterio/_err.pyx | 3 +-
rasterio/_gdal.pxd | 5 +
rasterio/_io.pxd | 1 +
rasterio/_io.pyx | 40 ++--
rasterio/_warp.pyx | 386 +++++++++++++++++++++++----------------
rasterio/crs.py | 2 +-
rasterio/enums.py | 1 +
rasterio/errors.py | 13 +-
rasterio/features.py | 4 +-
rasterio/rio/bands.py | 7 +-
rasterio/rio/calc.py | 7 +-
rasterio/rio/features.py | 135 ++++++--------
rasterio/rio/helpers.py | 26 ++-
rasterio/rio/info.py | 55 +-----
rasterio/rio/merge.py | 18 +-
rasterio/rio/options.py | 63 ++++++-
rasterio/rio/warp.py | 174 ++++++++++++------
rasterio/tool.py | 54 +++++-
rasterio/tools/mask.py | 92 ++++++++++
rasterio/transform.py | 1 +
rasterio/warp.py | 115 +++++++-----
tests/conftest.py | 36 ++++
tests/data/no_band.h5 | Bin 0 -> 9836 bytes
tests/data/world.byte.tif | Bin 0 -> 54885 bytes
tests/test_crs.py | 20 +-
tests/test_read.py | 7 +
tests/test_rio_features.py | 4 +-
tests/test_rio_helpers.py | 31 +++-
tests/test_rio_info.py | 42 +----
tests/test_rio_merge.py | 6 +-
tests/test_rio_options.py | 33 +++-
tests/test_rio_warp.py | 93 +++++++++-
tests/test_tool.py | 45 ++++-
tests/test_tools_mask.py | 51 ++++++
tests/test_transform.py | 21 +++
tests/test_warp.py | 78 ++++++--
tests/test_warp_transform.py | 63 +++++++
55 files changed, 1304 insertions(+), 583 deletions(-)
diff --git a/AUTHORS.txt b/AUTHORS.txt
index 5050878..3e35b11 100644
--- a/AUTHORS.txt
+++ b/AUTHORS.txt
@@ -1,25 +1,34 @@
Authors
=======
-Sean Gillies <sean at mapbox.com>
-Brendan Ward <bcward at consbio.org>
+Aldo Culquicondor <alculquicondor at gmail.com>
+Alessandro Amici <alexamici at gmail.com>
+Alexander <spatial.hast at gmail.com>
Amit Kapadia <amit at planet.com>
+AsgerPetersen <asgerpetersen at gmail.com>
+Bas Couwenberg <sebastic at xs4all.nl>
+Brendan Ward <bcward at consbio.org>
+Etienne B. Racine <etiennebr at gmail.com>
+Even Rouault <even.rouault at spatialys.com>
+Jacques Tardie <hi at jacquestardie.org>
+James McBride <jmcbride at berkeley.edu>
+James Seppi <james.seppi at gmail.com>
+Jeffrey Gerard <jgerard at climate.com>
+Johan Van de Wauw <johan.vandewauw at gmail.com>
+Joshua Arnott <josh at snorfalorpagus.net>
Kelsey Jordahl <kjordahl at alum.mit.edu>
Kevin Wurster <wursterk at gmail.com>
+Martijn Visser <mgvisser at gmail.com>
+Matt Savoie <github at flamingbear.com>
+Matthew Perry <perrygeo at gmail.com>
Maxim Dubinin <sim at gis-lab.info>
-Ryan Grout <rgrout at continuum.io>
Mike Toews <mwtoews at gmail.com>
-AsgerPetersen <asgerpetersen at gmail.com>
-Joshua Arnott <josh at snorfalorpagus.net>
-Alessandro Amici <alexamici at gmail.com>
-Johan Van de Wauw <johan.vandewauw at gmail.com>
-James Seppi <james.seppi at gmail.com>
-Jacques Tardie <hi at jacquestardie.org>
-Etienne B. Racine <etiennebr at gmail.com>
-Christoph Gohlke <cgohlke at uci.edu>
-Martijn Visser <mgvisser at gmail.com>
-Aldo Culquicondor <alculquicondor at gmail.com>
-Robin Wilson <robin at rtwilson.com>
+Nat Wilson <njwilson23 at gmail.com>
Patrick Young <patrick.young at digitalglobe.com>
+Robin Wilson <robin at rtwilson.com>
+Ryan Grout <rgrout at continuum.io>
+Sean Gillies <sean at mapbox.com>
+Trevor R.H. Clarke <tclarke at ball.com>
+cgohlke <cgohlke at uci.edu>
See also https://github.com/mapbox/rasterio/graphs/contributors.
diff --git a/CHANGES.txt b/CHANGES.txt
index 19f4033..74411ee 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,6 +1,41 @@
Changes
=======
+0.32.0 (2016-03-22)
+-------------------
+- Bug fix: geometry factories and warp operations are properly deallocated
+ in normal and error situations (#494, #568).
+- Bug fix: a code block in rio-merge's help has been better formatted (#535).
+- Bug fix: the rasterio.vfs module is imported in __init__.py to assist
+ cx_Freeze (#536).
+- Bug fix: old usage of `read_band()` has been replaced by `read()` throughout
+ the docs (#537).
+- Bug fix: accidental overwriting of existing files is now prevented by the
+ `resolve_inout()` function in `rasterio.rio.helpers`. Commands that take
+ one or more input files plus an output file should use this helper and force
+ overwrite either by using a `--force-overwrite` option or by using the
+ `-o output` option, which implicitly forces overwriting (#539, #540).
+- Bug fix: missing support for NaN nodata value in rio-warp added (#542, #544).
+- Bug fix: missing documentation of `rasterize()`'s `fill` parameter added
+ (#543).
+- Bug fix: raster dataset bounds are densified before transforming so that
+ the projected output of rio-bounds is correct (#556, #557).
+- Bug fix: add 'line' to the `Interleaving` enum (#560).
+- Bug fix: convert `matplotlib` import errors to a `RuntimeWarning` (#562).
+- Bug fix: deallocate CPL strings in error cases (#573).
+- Bug fix: non-invertable affine transforms are prevented using
+ `__future__.division` *#580).
+- Bug fix: rio-warp clips output regions to the limits of the destination
+ CRS unless disabled with `--no-check-invert-proj` (#597).
+- New feature: the functionality previously available only in rio-mask is now
+ available as `rasterio.tools.mask.mask()` (#552).
+- New feature: raster bounds are used to label axes in `rasterio.tool.show()`
+ (#553).
+- New feature: GDAL's suggested warp bounds algorithm is wrapped and exposed
+ for use in `warp()` and rio-warp (#574).
+- Breaking change: align rio-warp's `--bounds` option with rio-merge's: these
+ are in destination CRS units (#541, #545).
+
0.31.0 (2015-12-18)
-------------------
- Warn when rasters have no georeferencing and when the default identity
@@ -19,7 +54,7 @@ Changes
window_intersection(), windows_intersect() (#496, #506).
- Warn when an alpha band that might provide a dataset mask is shadowed by a
nodata attribute (#508, #523).
-- IPython is not the default interpeter for rio-insp and the documentation
+- IPython is not the default interpreter for rio-insp and the documentation
saying it is has been corrected (#518).
- Guard against creating datasets with block sizes larger than the dataset
width and height. Such datasets are semi-broken and are likely to be
diff --git a/README.rst b/README.rst
index 51ddcc0..8cc6fc7 100644
--- a/README.rst
+++ b/README.rst
@@ -117,7 +117,7 @@ using Python.
$ rio insp tests/data/RGB.byte.tif
Rasterio 0.10 Interactive Inspector (Python 3.4.1)
- Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
+ Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>> src.name
'tests/data/RGB.byte.tif'
>>> src.closed
diff --git a/docs/cli.rst b/docs/cli.rst
index 86e591e..a1366ed 100644
--- a/docs/cli.rst
+++ b/docs/cli.rst
@@ -317,7 +317,7 @@ option.
.. code-block:: console
- $ rio info tests/data/RGB.byte.tif --indent 2
+ $ rio info tests/data/RGB.byte.tif --indent 2 --verbose
{
"count": 3,
"crs": "EPSG:32618",
@@ -373,7 +373,7 @@ The ``insp`` command opens a dataset and an interpreter.
$ rio insp tests/data/RGB.byte.tif
Rasterio 0.18 Interactive Inspector (Python 2.7.9)
- Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
+ Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>> print src.name
tests/data/RGB.byte.tif
>>> print src.bounds
diff --git a/docs/colormaps.rst b/docs/colormaps.rst
index 45af887..a7e8cdd 100644
--- a/docs/colormaps.rst
+++ b/docs/colormaps.rst
@@ -14,7 +14,7 @@ to bands using the ``write_colormap()`` method.
with rasterio.drivers():
with rasterio.open('tests/data/shade.tif') as src:
- shade = src.read_band(1)
+ shade = src.read(1)
meta = src.meta
with rasterio.open('/tmp/colormap.tif', 'w', **meta) as dst:
diff --git a/docs/concurrency.rst b/docs/concurrency.rst
index 3d43859..14b4963 100644
--- a/docs/concurrency.rst
+++ b/docs/concurrency.rst
@@ -76,8 +76,6 @@ Here is the program in examples/concurrent-cpu-bound.py.
with rasterio.open(outfile, 'w', **meta) as dst:
# Define a generator for data, window pairs.
- # We use the new read() method here to a 3D array with all
- # bands, but could also use read_band().
def jobs():
for ij, window in dst.block_windows():
data = src.read(window=window)
diff --git a/docs/datasets.rst b/docs/datasets.rst
index e9d097a..12a3870 100644
--- a/docs/datasets.rst
+++ b/docs/datasets.rst
@@ -89,7 +89,7 @@ a file can be read like this:
.. code-block:: pycon
- >>> dataset.read_band(1)
+ >>> dataset.read(1)
array([[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
@@ -106,7 +106,7 @@ elsewhere.
.. code-block::
>>> from matplotlib import pyplot
- >>> pyplot.imshow(dataset.read_band(1), cmap='pink')
+ >>> pyplot.imshow(dataset.read(1), cmap='pink')
<matplotlib.image.AxesImage object at 0x111195c10>
>>> pyplot.show()
@@ -136,7 +136,7 @@ After it's closed, data can no longer be read.
.. code-block:: pycon
- >>> dataset.read_band(1)
+ >>> dataset.read(1)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ValueError: can't read closed raster file
@@ -179,7 +179,7 @@ data types, and the specific format must be specified.
.. code-block:: pycon
- >>> with rasterio.oepn
+ >>> with rasterio.open
Writing data mostly works as with a Python file. There are a few format-
specific differences. TODO: details.
diff --git a/docs/features.rst b/docs/features.rst
index ef51d65..2bd51eb 100644
--- a/docs/features.rst
+++ b/docs/features.rst
@@ -23,7 +23,7 @@ The shapes of the foreground features can be extracted like this:
from rasterio import features
with rasterio.open('13547682814_f2e459f7a5_o_d.png') as src:
- blue = src.read_band(3)
+ blue = src.read(3)
mask = blue != 255
shapes = features.shapes(blue, mask=mask)
@@ -72,7 +72,9 @@ image to be created.
transform=src.transform)
The values for the input shapes are replaced with ``255`` in a generator
-expression. The resulting image, written to disk like this,
+expression. Areas not covered by input geometries are replaced with an
+optional ``fill`` value, which defaults to ``0``. The resulting image,
+written to disk like this,
.. code-block:: python
diff --git a/docs/georeferencing.rst b/docs/georeferencing.rst
index f8c2321..ae34617 100644
--- a/docs/georeferencing.rst
+++ b/docs/georeferencing.rst
@@ -16,7 +16,7 @@ Rasterio distribution root to see.
.. code-block:: pycon
Rasterio 0.9 Interactive Inspector (Python 3.4.1)
- Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
+ Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>> src
<open RasterReader name='tests/data/RGB.byte.tif' mode='r'>
>>> src.crs
@@ -42,7 +42,7 @@ argument.
Coordinate Transformation
-------------------------
-A dataset's pixel coordinate system has its orgin at the "upper left" (imagine
+A dataset's pixel coordinate system has its origin at the "upper left" (imagine
it displayed on your screen). Column index increases to the right, and row
index increases downward. The mapping of these coordinates to "world"
coordinates in the dataset's reference system is done with an affine
diff --git a/docs/masks.rst b/docs/masks.rst
index fb4da60..0cea876 100644
--- a/docs/masks.rst
+++ b/docs/masks.rst
@@ -29,7 +29,7 @@ inverse relationship in the context of RGB.byte.tif.
$ 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.
+ Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>> src.shape
(718, 791)
>>> src.count
@@ -99,7 +99,7 @@ copy of the test data opened using rio-insp in "r+" (update) mode.
$ 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.
+ Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>>
To mark that all pixels of all bands are valid (i.e., to override nodata
@@ -132,7 +132,7 @@ as an RGB image (with the help of `numpy.dstack()`):
$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.
+ Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>> msk = src.read_masks()
>>> show(np.dstack(msk))
@@ -179,7 +179,7 @@ If you want, you can read dataset bands as numpy masked arrays.
$ 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.
+ Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>> blue = src.read(1, masked=True)
>>> blue.mask
array([[ True, True, True, ..., True, True, True],
diff --git a/docs/tags.rst b/docs/tags.rst
index e08f6cd..541709d 100644
--- a/docs/tags.rst
+++ b/docs/tags.rst
@@ -15,7 +15,7 @@ I'm going to use the rasterio interactive inspector in these examples below.
$ rasterio.insp tests/data/RGB.byte.tif
Rasterio 0.6 Interactive Inspector (Python 2.7.5)
- Type "src.name", "src.read_band(1)", or "help(src)" for more information.
+ Type "src.name", "src.read(1)", or "help(src)" for more information.
>>>
Tags belong to namespaces. To get a copy of a dataset's tags from the default
diff --git a/docs/windowed-rw.rst b/docs/windowed-rw.rst
index fe5d497..0393cef 100644
--- a/docs/windowed-rw.rst
+++ b/docs/windowed-rw.rst
@@ -67,7 +67,7 @@ test file.
>>> import rasterio
>>> with rasterio.open('tests/data/RGB.byte.tif') as src:
- ... w = src.read_band(1, window=((0, 100), (0, 100)))
+ ... w = src.read(1, window=((0, 100), (0, 100)))
...
>>> print(w.shape)
(100, 100)
@@ -104,7 +104,7 @@ Below, the window is scaled to one third of the source image.
.. code-block:: python
with rasterio.open('tests/data/RGB.byte.tif') as src:
- b, g, r = (src.read_band(k) for k in (1, 2, 3))
+ b, g, r = (src.read(k) for k in (1, 2, 3))
write_window = (30, 269), (50, 313)
@@ -144,7 +144,7 @@ destination dataset.
read_window = (350, 410), (350, 450)
with rasterio.open('tests/data/RGB.byte.tif') as src:
- b, g, r = (src.read_band(k, window=read_window) for k in (1, 2, 3))
+ b, g, r = (src.read(k, window=read_window) for k in (1, 2, 3))
write_window = (-240, None), (-400, None)
@@ -242,7 +242,7 @@ The block windows themselves can be had from the block_windows function.
...
This function returns an iterator that yields a pair of values. The second is
-a window tuple that can be used in calls to read_band or write_band. The first
+a window tuple that can be used in calls to `read` or `write_band`. The first
is the pair of row and column indexes of this block within all blocks of the
dataset.
@@ -252,7 +252,7 @@ You may read windows of data from a file block-by-block like this.
>>> with rasterio.open('tests/data/RGB.byte.tif') as src:
... for ji, window in src.block_windows(1):
- ... r = src.read_band(1, window=window)
+ ... r = src.read(1, window=window)
... print(r.shape)
... break
...
@@ -266,7 +266,7 @@ it's a good idea to test this assumption in your code.
>>> with rasterio.open('tests/data/RGB.byte.tif') as src:
... assert len(set(src.block_shapes)) == 1
... for ji, window in src.block_windows(1):
- ... b, g, r = (src.read_band(k, window=window) for k in (1, 2, 3))
+ ... b, g, r = (src.read(k, window=window) for k in (1, 2, 3))
... print(ji, r.shape, g.shape, b.shape)
... break
...
diff --git a/examples/concurrent-cpu-bound.py b/examples/concurrent-cpu-bound.py
index 1ceeb7c..98de30c 100644
--- a/examples/concurrent-cpu-bound.py
+++ b/examples/concurrent-cpu-bound.py
@@ -32,8 +32,6 @@ def main(infile, outfile, num_workers=4):
with rasterio.open(outfile, 'w', **meta) as dst:
# Define a generator for data, window pairs.
- # We use the new read() method here to a 3D array with all
- # bands, but could also use read_band().
def jobs():
for ij, window in dst.block_windows():
data = src.read(window=window)
diff --git a/examples/decimate.py b/examples/decimate.py
index 12b2d1a..cd0c440 100644
--- a/examples/decimate.py
+++ b/examples/decimate.py
@@ -7,7 +7,7 @@ import rasterio
with rasterio.drivers():
with rasterio.open('tests/data/RGB.byte.tif') as src:
- b, g, r = (src.read_band(k) for k in (1, 2, 3))
+ b, g, r = (src.read(k) for k in (1, 2, 3))
meta = src.meta
tmpfilename = os.path.join(tempfile.mkdtemp(), 'decimate.tif')
diff --git a/examples/rasterio_polygonize.py b/examples/rasterio_polygonize.py
index 54dd19a..fbd4089 100644
--- a/examples/rasterio_polygonize.py
+++ b/examples/rasterio_polygonize.py
@@ -20,7 +20,7 @@ def main(raster_file, vector_file, driver, mask_value):
with rasterio.drivers():
with rasterio.open(raster_file) as src:
- image = src.read_band(1)
+ image = src.read(1)
if mask_value is not None:
mask = image == mask_value
diff --git a/examples/sieve.py b/examples/sieve.py
index b3aab7f..d074745 100644
--- a/examples/sieve.py
+++ b/examples/sieve.py
@@ -14,7 +14,7 @@ with rasterio.drivers():
# Read a raster to be sieved.
with rasterio.open('tests/data/shade.tif') as src:
- shade = src.read_band(1)
+ shade = src.read(1)
# Print the number of shapes in the source raster.
print("Slope shapes: %d" % len(list(shapes(shade))))
diff --git a/rasterio/__init__.py b/rasterio/__init__.py
index 142ce75..1ae488f 100644
--- a/rasterio/__init__.py
+++ b/rasterio/__init__.py
@@ -17,13 +17,13 @@ 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
+from rasterio import _err, coords, enums, vfs
# Classes in rasterio._io are imported below just before we need them.
__all__ = [
'band', 'open', 'drivers', 'copy', 'pad']
-__version__ = "0.31.0"
+__version__ = "0.32.0"
log = logging.getLogger('rasterio')
class NullHandler(logging.Handler):
diff --git a/rasterio/_base.pyx b/rasterio/_base.pyx
index ea59ac0..d41b11e 100644
--- a/rasterio/_base.pyx
+++ b/rasterio/_base.pyx
@@ -329,6 +329,8 @@ cdef class DatasetReader(object):
property nodata:
"""The dataset's single nodata value."""
def __get__(self):
+ if self.count == 0:
+ return None
return self.nodatavals[0]
property mask_flags:
@@ -358,7 +360,7 @@ cdef class DatasetReader(object):
indexes.
The primary use of this function is to obtain windows to pass to
- read_band() for highly efficient access to raster block data.
+ read() for highly efficient access to raster block data.
"""
cdef int i, j
block_shapes = self.block_shapes
@@ -443,9 +445,13 @@ cdef class DatasetReader(object):
@property
def meta(self):
"""The basic metadata of this dataset."""
+ if self.count == 0:
+ dtype = 'float_'
+ else:
+ dtype = self.dtypes[0]
m = {
'driver': self.driver,
- 'dtype': self.dtypes[0],
+ 'dtype': dtype,
'nodata': self.nodata,
'width': self.width,
'height': self.height,
@@ -475,6 +481,8 @@ cdef class DatasetReader(object):
@property
def is_tiled(self):
+ if len(self.block_shapes) == 0:
+ return False
return self.block_shapes[0][1] != self.width
property profile:
@@ -902,6 +910,13 @@ def _transform(src_crs, dst_crs, xs, ys, zs):
z[i] = zs[i]
transform = _gdal.OCTNewCoordinateTransformation(src, dst)
+ if transform == NULL:
+ _gdal.CPLFree(x)
+ _gdal.CPLFree(y)
+ _gdal.CPLFree(z)
+ _gdal.OSRDestroySpatialReference(src)
+ _gdal.OSRDestroySpatialReference(dst)
+ raise ValueError("Cannot create coordinate transformer")
res = _gdal.OCTTransform(transform, n, x, y, z)
#if res:
# raise ValueError("Failed coordinate transformation")
diff --git a/rasterio/_err.pyx b/rasterio/_err.pyx
index bae5c16..5d02c86 100644
--- a/rasterio/_err.pyx
+++ b/rasterio/_err.pyx
@@ -69,11 +69,12 @@ cdef class GDALErrCtxManager:
if err_type >= 3:
raise exception_map[err_no](msg)
+
cpl_errs = GDALErrCtxManager()
class GDALError(IntEnum):
- none = 0, # CE_None
+ success = 0, # CE_None
debug = 1, # CE_Debug
warning= 2, # CE_Warning
failure = 3, # CE_Failure
diff --git a/rasterio/_gdal.pxd b/rasterio/_gdal.pxd
index 75b937c..3b21e6b 100644
--- a/rasterio/_gdal.pxd
+++ b/rasterio/_gdal.pxd
@@ -224,6 +224,9 @@ cdef extern from "gdal_alg.h":
void* hDstDS, const char *pszDstWKT,
int bGCPUseOK, double dfGCPErrorThreshold,
int nOrder )
+ void *GDALCreateGenImgProjTransformer3(const char *pszSrcWKT,
+ const double *padfSrcGeoTransform, const char *pszDstWKT,
+ const double *padfDstGeoTransform)
int GDALGenImgProjTransform(void *pTransformArg, int bDstToSrc, int nPointCount, double *x, double *y, double *z, int *panSuccess )
void GDALDestroyGenImgProjTransformer( void * )
@@ -233,3 +236,5 @@ cdef extern from "gdal_alg.h":
int GDALFillNodata(void *dst_band, void *mask_band, double max_search_distance, int deprecated, int smoothing_iterations, char **options, void *progress_func, void *progress_data)
int GDALChecksumImage(void *band, int xoff, int yoff, int width, int height)
+
+ int GDALSuggestedWarpOutput2(void * hSrcDS, GDALTransformerFunc pfnRawTransformer, void * pTransformArg, double * padfGeoTransformOut, int * pnPixels, int * pnLines, double * padfExtent, int nOptions)
diff --git a/rasterio/_io.pxd b/rasterio/_io.pxd
index 9e5165f..63534b6 100644
--- a/rasterio/_io.pxd
+++ b/rasterio/_io.pxd
@@ -51,6 +51,7 @@ cdef class InMemoryRaster:
cdef double transform[6]
cdef int band_ids[1]
cdef np.ndarray _image
+ cdef object crs
ctypedef np.uint8_t DTYPE_UBYTE_t
diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx
index 9b003d6..1b5b2e7 100644
--- a/rasterio/_io.pyx
+++ b/rasterio/_io.pyx
@@ -20,6 +20,7 @@ from rasterio._drivers import driver_count, GDALEnv
from rasterio._err import cpl_errs, GDALError
from rasterio import dtypes
from rasterio.coords import BoundingBox
+from rasterio.errors import DriverRegistrationError
from rasterio.five import text_type, string_types
from rasterio.transform import Affine
from rasterio.enums import ColorInterp, MaskFlags, Resampling
@@ -29,16 +30,6 @@ from rasterio.warnings import NodataShadowWarning
log = logging.getLogger('rasterio')
-if 'all' in sys.warnoptions:
- # show messages in console with: python -W all
- logging.basicConfig()
-else:
- # no handler messages shown
- class NullHandler(logging.Handler):
- def emit(self, record):
- pass
-
-log.addHandler(NullHandler())
cdef bint in_dtype_range(value, dtype):
@@ -54,7 +45,11 @@ cdef bint in_dtype_range(value, dtype):
105: np.iinfo,
117: np.iinfo
}
- rng = infos[np.dtype(dtype).kind](dtype)
+ key = np.dtype(dtype).kind
+ if np.isnan(value):
+ return key in ('c', 'f', 99, 102)
+
+ rng = infos[key](dtype)
return rng.min <= value <= rng.max
# Single band IO functions.
@@ -1454,7 +1449,7 @@ cdef class RasterUpdater(RasterReader):
if [abs(v) for v in transform] == [0, 1, 0, 0, 0, 1]:
warnings.warn(
"Dataset uses default geotransform (Affine.identity). "
- "No tranform will be written to the output by GDAL.",
+ "No transform will be written to the output by GDAL.",
UserWarning
)
@@ -1834,7 +1829,7 @@ cdef class InMemoryRaster:
IO with GDAL. Other memory based operations should use numpy arrays.
"""
- def __cinit__(self, image, transform=None):
+ def __cinit__(self, image, transform=None, crs=None):
"""
Create in-memory raster dataset, and populate its initial values with
the values in image.
@@ -1847,12 +1842,18 @@ cdef class InMemoryRaster:
self._image = image
self.dataset = NULL
- cdef void *memdriver = _gdal.GDALGetDriverByName("MEM")
cdef int i = 0 # avoids Cython warning in for loop below
+ cdef const char *srcwkt = NULL
+ cdef void *osr = NULL
# Several GDAL operations require the array of band IDs as input
self.band_ids[0] = 1
+ cdef void *memdriver = _gdal.GDALGetDriverByName("MEM")
+ if memdriver == NULL:
+ raise DriverRegistrationError(
+ "MEM driver is not registered.")
+
self.dataset = _gdal.GDALCreate(
memdriver,
"output",
@@ -1873,6 +1874,17 @@ cdef class InMemoryRaster:
if err:
raise ValueError("transform not set: %s" % transform)
+ # Set projection if specified (for use with
+ # GDALSuggestedWarpOutput2()).
+ if crs:
+ osr = _base._osr_from_crs(crs)
+ _gdal.OSRExportToWkt(osr, &srcwkt)
+ _gdal.GDALSetProjection(self.dataset, srcwkt)
+ log.debug("Set CRS on temp source dataset: %s", srcwkt)
+ _gdal.CPLFree(srcwkt)
+ _gdal.OSRDestroySpatialReference(osr)
+
+
self.band = _gdal.GDALGetRasterBand(self.dataset, 1)
if self.band == NULL:
raise ValueError("NULL output band: {0}".format(i))
diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx
index 7d4e4f4..e2ad940 100644
--- a/rasterio/_warp.pyx
+++ b/rasterio/_warp.pyx
@@ -1,6 +1,6 @@
# distutils: language = c++
-from collections import namedtuple
+from enum import IntEnum
import logging
import numpy as np
@@ -8,7 +8,10 @@ cimport numpy as np
from rasterio cimport _base, _gdal, _ogr, _io, _features
from rasterio import dtypes
-from rasterio.errors import RasterioDriverRegistrationError
+from rasterio._err import cpl_errs, GDALError
+from rasterio._io cimport InMemoryRaster
+from rasterio.errors import DriverRegistrationError
+from rasterio.transform import Affine, from_bounds
cdef extern from "gdalwarper.h" nogil:
@@ -38,14 +41,19 @@ cdef extern from "gdalwarper.h" nogil:
double dfProgressScale=1.0)
-RESAMPLING = namedtuple('RESAMPLING', [
- 'nearest',
- 'bilinear',
- 'cubic',
- 'cubic_spline',
- 'lanczos',
- 'average',
- 'mode'] )(*list(range(7)))
+class Resampling(IntEnum):
+ nearest=0
+ bilinear=1
+ cubic=2
+ cubic_spline=3
+ lanczos=4
+ average=5
+ mode=6
+ max=8
+ min=9
+ med=10
+ q1=11
+ q3=12
cdef extern from "ogr_geometry.h" nogil:
@@ -55,9 +63,6 @@ cdef extern from "ogr_geometry.h" nogil:
cdef cppclass OGRGeometryFactory:
void * transformWithOptions(void *geom, void *ct, char **options)
-# const OGRGeometry* poSrcGeom,
-# OGRCoordinateTransformation *poCT,
-# char** papszOptions
cdef extern from "ogr_spatialref.h":
@@ -111,6 +116,7 @@ def _transform_geom(
<const OGRGeometry *>src_ogr_geom,
<OGRCoordinateTransformation *>transform,
options)
+ del factory
g = _features.GeomBuilder().build(dst_ogr_geom)
_ogr.OGR_G_DestroyGeometry(dst_ogr_geom)
@@ -163,7 +169,8 @@ def _reproject(
dst_transform=None,
dst_crs=None,
dst_nodata=None,
- resampling=RESAMPLING.nearest,
+ resampling=Resampling.nearest,
+ num_threads=1,
**kwargs):
"""
Reproject a source raster to a destination raster.
@@ -208,13 +215,15 @@ def _reproject(
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
+ Resampling.nearest,
+ Resampling.bilinear,
+ Resampling.cubic,
+ Resampling.cubic_spline,
+ Resampling.lanczos,
+ Resampling.average,
+ Resampling.mode
+ num_threads: int
+ Number of worker threads.
kwargs: dict, optional
Additional arguments passed to transformation function.
@@ -260,7 +269,7 @@ def _reproject(
hrdriver = _gdal.GDALGetDriverByName("MEM")
if hrdriver == NULL:
- raise RasterioDriverRegistrationError(
+ raise DriverRegistrationError(
"'MEM' driver not found. Check that this call is contained "
"in a `with rasterio.drivers()` or `with rasterio.open()` "
"block.")
@@ -280,9 +289,9 @@ def _reproject(
osr = _base._osr_from_crs(src_crs)
_gdal.OSRExportToWkt(osr, &srcwkt)
_gdal.GDALSetProjection(hdsin, srcwkt)
+ log.debug("Set CRS on temp source dataset: %s", srcwkt)
_gdal.CPLFree(srcwkt)
_gdal.OSRDestroySpatialReference(osr)
- log.debug("Set CRS on temp source dataset: %s", srcwkt)
# Copy arrays to the dataset.
retval = _io.io_auto(source, hdsin, 1)
@@ -308,7 +317,7 @@ def _reproject(
hrdriver = _gdal.GDALGetDriverByName("MEM")
if hrdriver == NULL:
- raise RasterioDriverRegistrationError(
+ raise DriverRegistrationError(
"'MEM' driver not found. Check that this call is contained "
"in a `with rasterio.drivers()` or `with rasterio.open()` "
"block.")
@@ -318,21 +327,29 @@ def _reproject(
hrdriver, "output", cols, rows, src_count,
dtypes.dtype_rev[np.dtype(destination.dtype).name], NULL)
if hdsout == NULL:
- raise ValueError("NULL output datasource")
+ raise ValueError("Failed to create temp destination dataset.")
_gdal.GDALSetDescription(
hdsout, "Temporary destination dataset for _reproject()")
- log.debug("Created temp destination dataset")
+ log.debug("Created temp destination dataset.")
+
for i in range(6):
gt[i] = dst_transform[i]
- retval = _gdal.GDALSetGeoTransform(hdsout, gt)
- log.debug("Set transform on temp destination dataset: %d", retval)
+
+ if not GDALError.success == _gdal.GDALSetGeoTransform(hdsout, gt):
+ raise ValueError(
+ "Failed to set transform on temp destination dataset.")
+
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)
+ log.debug("CRS for temp destination dataset: %s.", dstwkt)
+
+ if not GDALError.success == _gdal.GDALSetProjection(hdsout, dstwkt):
+ raise ValueError(
+ "Failed to set projection on temp destination dataset.")
+
+ _gdal.CPLFree(dstwkt)
+
if dst_nodata is None and hasattr(destination, "fill_value"):
# destination is a masked array
dst_nodata = destination.fill_value
@@ -347,156 +364,201 @@ def _reproject(
cdef void *hTransformArg = NULL
cdef _gdal.GDALWarpOptions *psWOptions = NULL
- cdef GDALWarpOperation *oWarper = new GDALWarpOperation()
- reprojected = False
- try:
- hTransformArg = _gdal.GDALCreateGenImgProjTransformer(
- hdsin, NULL, hdsout, NULL,
- 1, 1000.0, 0)
- if hTransformArg == NULL:
- raise ValueError("NULL transformer")
- log.debug("Created transformer")
+ hTransformArg = _gdal.GDALCreateGenImgProjTransformer(
+ hdsin, NULL, hdsout, NULL,
+ 1, 1000.0, 0)
+ if hTransformArg == NULL:
+ raise ValueError("NULL transformer")
+ 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
- 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')
- val_b = v.encode('utf-8')
- 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))
+ val_b = str(num_threads).encode('utf-8')
+ warp_extras = _gdal.CSLSetNameValue(warp_extras, "NUM_THREADS", val_b)
+ log.debug("Setting NUM_THREADS option: %s", val_b)
- pszWarpThreads = _gdal.CSLFetchNameValue(warp_extras, "NUM_THREADS")
- if pszWarpThreads == NULL:
- pszWarpThreads = _gdal.CPLGetConfigOption(
- "GDAL_NUM_THREADS", "1")
- warp_extras = _gdal.CSLSetNameValue(
- warp_extras, "NUM_THREADS", pszWarpThreads)
-
- log.debug("Created warp options")
-
- psWOptions.eResampleAlg = <_gdal.GDALResampleAlg>resampling
+ for k, v in kwargs.items():
+ k, v = k.upper(), str(v).upper()
+ key_b = k.encode('utf-8')
+ val_b = v.encode('utf-8')
+ 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))
- # 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)
+ log.debug("Created warp options")
- 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)
+ psWOptions.eResampleAlg = <_gdal.GDALResampleAlg>resampling
- # Validate nodata values
+ # Set src_nodata and dst_nodata
+ if src_nodata is None and dst_nodata is not None:
+ psWOptions.papszWarpOptions = warp_extras
+ _gdal.GDALDestroyGenImgProjTransformer(hTransformArg)
+ _gdal.GDALDestroyWarpOptions(psWOptions)
+ 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:
- 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))
+ 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):
+ psWOptions.papszWarpOptions = warp_extras
+ _gdal.GDALDestroyGenImgProjTransformer(hTransformArg)
+ _gdal.GDALDestroyWarpOptions(psWOptions)
+ 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.padfDstNoDataReal[i] = dst_nodata
- psWOptions.padfDstNoDataImag[i] = 0.0
+ psWOptions.padfSrcNoDataReal[i] = src_nodata
+ psWOptions.padfSrcNoDataImag[i] = 0.0
warp_extras = _gdal.CSLSetNameValue(
- warp_extras, "INIT_DEST", "NO_DATA")
+ warp_extras, "UNIFIED_SRC_NODATA", "YES")
- # 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(
- # _gdal.GDALGenImgProjTransform,
- # hTransformArg,
- # maxerror )
- # psWOptions.pfnTransformer = _gdal.GDALApproxTransform
- #else:
- psWOptions.pfnTransformer = _gdal.GDALGenImgProjTransform
- psWOptions.pTransformerArg = hTransformArg
- psWOptions.hSrcDS = hdsin
- psWOptions.hDstDS = hdsout
- psWOptions.nBandCount = src_count
- psWOptions.panSrcBands = <int *>_gdal.CPLMalloc(src_count*sizeof(int))
- psWOptions.panDstBands = <int *>_gdal.CPLMalloc(src_count*sizeof(int))
- if isinstance(source, tuple):
- psWOptions.panSrcBands[0] = source.bidx
- else:
- for i in range(src_count):
- psWOptions.panSrcBands[i] = i+1
- if isinstance(destination, tuple):
- psWOptions.panDstBands[0] = destination.bidx
- else:
- for i in range(src_count):
- psWOptions.panDstBands[i] = i+1
- log.debug("Set transformer options")
+ if dst_nodata is not None and not _io.in_dtype_range(
+ dst_nodata, destination.dtype):
+ psWOptions.papszWarpOptions = warp_extras
+ _gdal.GDALDestroyGenImgProjTransformer(hTransformArg)
+ _gdal.GDALDestroyWarpOptions(psWOptions)
+ 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
+
+ psWOptions.pfnTransformer = _gdal.GDALGenImgProjTransform
+ psWOptions.pTransformerArg = hTransformArg
+ psWOptions.hSrcDS = hdsin
+ psWOptions.hDstDS = hdsout
+ psWOptions.nBandCount = src_count
+ psWOptions.panSrcBands = <int *>_gdal.CPLMalloc(src_count*sizeof(int))
+ psWOptions.panDstBands = <int *>_gdal.CPLMalloc(src_count*sizeof(int))
+ if isinstance(source, tuple):
+ psWOptions.panSrcBands[0] = source.bidx
+ else:
+ for i in range(src_count):
+ psWOptions.panSrcBands[i] = i+1
+ if isinstance(destination, tuple):
+ psWOptions.panDstBands[0] = destination.bidx
+ else:
+ for i in range(src_count):
+ psWOptions.panDstBands[i] = i+1
+ log.debug("Set transformer options")
- # TODO: alpha band.
+ # TODO: alpha band.
- eErr = oWarper.Initialize(psWOptions)
- if eErr == 0:
+ # Now that the transformer and warp options are set up, we init
+ # and run the warper.
+ cdef GDALWarpOperation *oWarper = new GDALWarpOperation()
+ try:
+ with cpl_errs:
+ oWarper.Initialize(psWOptions)
+
+ rows, cols = destination.shape[-2:]
+ log.debug(
+ "Chunk and warp window: %d, %d, %d, %d.",
+ 0, 0, cols, rows)
+
+ with cpl_errs:
+ if num_threads > 1:
+ log.debug("Executing multi warp with num_threads: %d", num_threads)
+ oWarper.ChunkAndWarpMulti(0, 0, cols, rows)
+ else:
+ oWarper.ChunkAndWarpImage(0, 0, cols, rows)
- log.debug("Destination shape: %r", destination.shape)
+ if dtypes.is_ndarray(destination):
+ retval = _io.io_auto(destination, hdsout, 0)
+ # TODO: handle errors (by retval).
- rows, cols = destination.shape[-2:]
- log.debug(
- "Chunk and warp window: %d, %d, %d, %d",
- 0, 0, cols, rows)
- with nogil:
- eErr = oWarper.ChunkAndWarpMulti(0, 0, cols, rows)
- log.debug("Chunked and warped: %d", retval)
-
- except Exception:
- log.exception(
- "Caught exception in warping. Source not reprojected.")
- raise
-
- else:
- reprojected = True
+ if hdsout != NULL:
+ _gdal.GDALClose(hdsout)
+ # Clean up transformer, warp options, and dataset handles.
finally:
- if hTransformArg != NULL:
- _gdal.GDALDestroyGenImgProjTransformer(hTransformArg)
- #if maxerror > 0.0:
- # _gdal.GDALDestroyApproxTransformer(psWOptions.pTransformerArg)
- if psWOptions != NULL:
- _gdal.GDALDestroyWarpOptions(psWOptions)
+ _gdal.GDALDestroyGenImgProjTransformer(hTransformArg)
+ _gdal.GDALDestroyWarpOptions(psWOptions)
if dtypes.is_ndarray(source):
if hdsin != NULL:
_gdal.GDALClose(hdsin)
- if reprojected and dtypes.is_ndarray(destination):
- retval = _io.io_auto(destination, hdsout, 0)
- # TODO: handle errors (by retval).
- if hdsout != NULL:
- _gdal.GDALClose(hdsout)
+def _calculate_default_transform(
+ src_crs, dst_crs, width, height, left, bottom, right, top, **kwargs):
+ """Wraps GDAL's algorithm."""
+
+ cdef void *hTransformArg = NULL
+ cdef int npixels = 0
+ cdef int nlines = 0
+ cdef double extent[4]
+ cdef double geotransform[6]
+ cdef void *osr = NULL
+ cdef char *wkt = NULL
+ cdef InMemoryRaster temp = None
+
+ extent[:] = [0.0, 0.0, 0.0, 0.0]
+ geotransform[:] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+ # Make an in-memory raster dataset we can pass to
+ # GDALCreateGenImgProjTransformer().
+ transform = from_bounds(left, bottom, right, top, width, height)
+ img = np.empty((height, width))
+
+ osr = _base._osr_from_crs(dst_crs)
+ _gdal.OSRExportToWkt(osr, &wkt)
+ _gdal.OSRDestroySpatialReference(osr)
+
+ with InMemoryRaster(
+ img, transform=transform.to_gdal(), crs=src_crs) as temp:
+ hTransformArg = _gdal.GDALCreateGenImgProjTransformer(
+ temp.dataset, NULL, NULL, wkt, 1, 1000.0, 0)
+ if hTransformArg == NULL:
+ if wkt != NULL:
+ _gdal.CPLFree(wkt)
+ raise ValueError("NULL transformer")
+ log.debug("Created transformer")
+ # geotransform, npixels, and nlines are modified by the
+ # function called below.
+ try:
+ if not GDALError.success == _gdal.GDALSuggestedWarpOutput2(
+ temp.dataset, _gdal.GDALGenImgProjTransform, hTransformArg,
+ geotransform, &npixels, &nlines, extent, 0):
+ raise RuntimeError(
+ "Failed to compute a suggested warp output.")
+ finally:
+ if wkt != NULL:
+ _gdal.CPLFree(wkt)
+ if hTransformArg != NULL:
+ _gdal.GDALDestroyGenImgProjTransformer(hTransformArg)
+
+ # Convert those modified arguments to Python values.
+ dst_affine = Affine.from_gdal(*[geotransform[i] for i in range(6)])
+ dst_width = npixels
+ dst_height = nlines
+
+ return dst_affine, dst_width, dst_height
diff --git a/rasterio/crs.py b/rasterio/crs.py
index 2ee4d45..5a04efe 100644
--- a/rasterio/crs.py
+++ b/rasterio/crs.py
@@ -97,7 +97,7 @@ def from_epsg(code):
# Below is the big list of PROJ4 parameters from
# http://trac.osgeo.org/proj/wiki/GenParms.
-# It is parsed into a list of paramter keys ``all_proj_keys``.
+# It is parsed into a list of parameter keys ``all_proj_keys``.
_param_data = """
+a Semimajor radius of the ellipsoid axis
diff --git a/rasterio/enums.py b/rasterio/enums.py
index c611052..68563ee 100644
--- a/rasterio/enums.py
+++ b/rasterio/enums.py
@@ -44,6 +44,7 @@ class Compression(Enum):
class Interleaving(Enum):
pixel='PIXEL'
+ line='LINE'
band='BAND'
diff --git a/rasterio/errors.py b/rasterio/errors.py
index 288fb0e..3a4f34e 100644
--- a/rasterio/errors.py
+++ b/rasterio/errors.py
@@ -1,9 +1,18 @@
"""A module of errors."""
+from click import FileError
+
class RasterioIOError(IOError):
"""A failure to open a dataset using the presently registered drivers."""
-class RasterioDriverRegistrationError(ValueError):
- """To be raised when, eg, _gdal.GDALGetDriverByName("MEM") returns NULL"""
+class DriverRegistrationError(ValueError):
+ """To be raised when, eg, _gdal.GDALGetDriverByName("MEM") returns NULL."""
+
+
+class FileOverwriteError(FileError):
+ """Rasterio's CLI refuses to implicitly clobber output files."""
+
+ def __init__(self, message):
+ super(FileOverwriteError, self).__init__('', hint=message)
diff --git a/rasterio/features.py b/rasterio/features.py
index e49aa76..59e9c9b 100644
--- a/rasterio/features.py
+++ b/rasterio/features.py
@@ -158,7 +158,7 @@ def sieve(image, size, out=None, output=None, mask=None, connectivity=4):
# Start moving users over to 'out'.
if output is not None:
warnings.warn(
- "The 'output' keyword arg has been superceded by 'out' "
+ "The 'output' keyword arg has been superseded by 'out' "
"and will be removed before Rasterio 1.0.",
FutureWarning,
stacklevel=2) # pragma: no cover
@@ -296,7 +296,7 @@ def rasterize(
if output is not None:
warnings.warn(
- "The 'output' keyword arg has been superceded by 'out' "
+ "The 'output' keyword arg has been superseded by 'out' "
"and will be removed before Rasterio 1.0.",
FutureWarning,
stacklevel=2) # pragma: no cover
diff --git a/rasterio/rio/bands.py b/rasterio/rio/bands.py
index f692644..53c3dae 100644
--- a/rasterio/rio/bands.py
+++ b/rasterio/rio/bands.py
@@ -16,9 +16,11 @@ from rasterio.five import zip_longest
@format_opt
@options.bidx_mult_opt
@options.rgb_opt
+ at options.force_overwrite_opt
@options.creation_options
@click.pass_context
-def stack(ctx, files, output, driver, bidx, photometric, creation_options):
+def stack(ctx, files, output, driver, bidx, photometric, force_overwrite,
+ creation_options):
"""Stack a number of bands from one or more input files into a
multiband dataset.
@@ -55,7 +57,8 @@ def stack(ctx, files, output, driver, bidx, photometric, creation_options):
logger = logging.getLogger('rio')
try:
with rasterio.drivers(CPL_DEBUG=verbosity>2):
- output, files = resolve_inout(files=files, output=output)
+ output, files = resolve_inout(files=files, output=output,
+ force_overwrite=force_overwrite)
output_count = 0
indexes = []
for path, item in zip_longest(files, bidx, fillvalue=None):
diff --git a/rasterio/rio/calc.py b/rasterio/rio/calc.py
index b9df4ff..4bc1227 100644
--- a/rasterio/rio/calc.py
+++ b/rasterio/rio/calc.py
@@ -40,9 +40,11 @@ def read_array(ix, subix=None, dtype=None):
'"a=tests/data/RGB.byte.tif".')
@options.dtype_opt
@options.masked_opt
+ at options.force_overwrite_opt
@options.creation_options
@click.pass_context
-def calc(ctx, command, files, output, name, dtype, masked, creation_options):
+def calc(ctx, command, files, output, name, dtype, masked, force_overwrite,
+ creation_options):
"""A raster data calculator
Evaluates an expression using input datasets and writes the result
@@ -89,7 +91,8 @@ def calc(ctx, command, files, output, name, dtype, masked, creation_options):
try:
with rasterio.drivers(CPL_DEBUG=verbosity > 2):
- output, files = resolve_inout(files=files, output=output)
+ output, files = resolve_inout(files=files, output=output,
+ force_overwrite=force_overwrite)
inputs = ([tuple(n.split('=')) for n in name] +
[(None, n) for n in files])
diff --git a/rasterio/rio/features.py b/rasterio/rio/features.py
index ebfcf19..e35cd66 100644
--- a/rasterio/rio/features.py
+++ b/rasterio/rio/features.py
@@ -2,8 +2,8 @@ import json
import logging
from math import ceil
import os
-import shutil
import re
+import shutil
import click
import cligj
@@ -11,18 +11,27 @@ from cligj import (
precision_opt, indent_opt, compact_opt, projection_geographic_opt,
projection_mercator_opt, projection_projected_opt, sequence_opt,
use_rs_opt, geojson_type_feature_opt, geojson_type_bbox_opt,
- files_inout_arg, format_opt, geojson_type_collection_opt)
+ format_opt, geojson_type_collection_opt)
from .helpers import coords, resolve_inout, write_features, to_lower
from . import options
import rasterio
from rasterio.transform import Affine
from rasterio.coords import disjoint_bounds
+from rasterio.warp import transform_bounds
logger = logging.getLogger('rio')
# Common options used below
+
+# Unlike the version in cligj, this one doesn't require values.
+files_inout_arg = click.argument(
+ 'files',
+ nargs=-1,
+ type=click.Path(resolve_path=True),
+ metavar="INPUTS... OUTPUT")
+
all_touched_opt = click.option(
'-a', '--all', '--all_touched', 'all_touched',
is_flag=True,
@@ -50,6 +59,7 @@ all_touched_opt = click.option(
help='Inverts the mask, so that areas covered by features are'
'masked out and areas not covered are retained. Ignored '
'if using --crop')
+ at options.force_overwrite_opt
@options.creation_options
@click.pass_context
def mask(
@@ -61,6 +71,7 @@ def mask(
all_touched,
crop,
invert,
+ force_overwrite,
creation_options):
"""Masks in raster using GeoJSON features (masks out all areas not covered
@@ -85,12 +96,14 @@ def mask(
input raster.
"""
- from rasterio.features import geometry_mask
+ from rasterio.tools.mask import mask as mask_tool
from rasterio.features import bounds as calculate_bounds
verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+ logger = logging.getLogger('rio')
- output, files = resolve_inout(files=files, output=output)
+ output, files = resolve_inout(
+ files=files, output=output, force_overwrite=force_overwrite)
input = files[0]
if geojson_mask is None:
@@ -113,7 +126,7 @@ def mask(
param_hint='--geojson-mask')
if 'features' in geojson:
- geometries = (f['geometry'] for f in geojson['features'])
+ geometries = [f['geometry'] for f in geojson['features']]
elif 'geometry' in geojson:
geometries = (geojson['geometry'], )
else:
@@ -122,62 +135,29 @@ def mask(
bounds = geojson.get('bbox', calculate_bounds(geojson))
with rasterio.open(input) as src:
- # If y pixel value is positive, then invert y dimension in bounds
- invert_y = src.affine.e > 0
-
- src_bounds = src.bounds
- if invert_y:
- src_bounds = [src.bounds[0], src.bounds[3],
- src.bounds[2], src.bounds[1]]
-
- has_disjoint_bounds = disjoint_bounds(bounds, src_bounds)
-
- if crop:
- if has_disjoint_bounds:
-
- raise click.BadParameter('not allowed for GeoJSON outside '
- 'the extent of the input raster',
- param=crop, param_hint='--crop')
-
- if invert_y:
- bounds = (bounds[0], bounds[3], bounds[2], bounds[1])
-
- window = src.window(*bounds)
- transform = src.window_transform(window)
- (r1, r2), (c1, c2) = window
- mask_shape = (r2 - r1, c2 - c1)
- else:
- if has_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)
+ try:
+ out_image, out_transform = mask_tool(src, geometries,
+ crop=crop, invert=invert,
+ all_touched=all_touched)
+ except ValueError as e:
+ if e.args[0] == 'Input shapes do not overlap raster.':
+ if crop:
+ raise click.BadParameter('not allowed for GeoJSON '
+ 'outside the extent of the '
+ 'input raster', param=crop,
+ param_hint='--crop')
meta = src.meta.copy()
meta.update(**creation_options)
meta.update({
'driver': driver,
- 'height': mask.shape[0],
- 'width': mask.shape[1],
- 'transform': transform
+ 'height': out_image.shape[1],
+ 'width': out_image.shape[2],
+ 'transform': out_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]))
+ out.write(out_image)
# Shapes command.
@@ -303,7 +283,7 @@ def shapes(
if bidx is None:
msk = numpy.logical_or.reduce(msk).astype('uint8')
- # Possibly overidden below.
+ # Possibly overridden below.
img = msk
# Read the band data unless the --mask option is given.
@@ -399,9 +379,10 @@ def shapes(
@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')
- at click.option('--property', type=str, default=None, help='Property in '
+ at click.option('--property', 'prop', 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.')
+ at options.force_overwrite_opt
@options.creation_options
@click.pass_context
def rasterize(
@@ -417,7 +398,8 @@ def rasterize(
all_touched,
default_value,
fill,
- property,
+ prop,
+ force_overwrite,
creation_options):
"""Rasterize GeoJSON into a new or existing raster.
@@ -466,7 +448,8 @@ def rasterize(
verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
- output, files = resolve_inout(files=files, output=output)
+ output, files = resolve_inout(
+ files=files, output=output, force_overwrite=force_overwrite)
has_src_crs = src_crs is not None
src_crs = src_crs or 'EPSG:4326'
@@ -481,8 +464,8 @@ def rasterize(
with rasterio.drivers(CPL_DEBUG=verbosity > 2):
def feature_value(feature):
- if property and 'properties' in feature:
- return feature['properties'].get(property, default_value)
+ if prop and 'properties' in feature:
+ return feature['properties'].get(prop, default_value)
return default_value
with click.open_file(files.pop(0) if files else '-') as gj_f:
@@ -524,7 +507,7 @@ def rasterize(
fill = fill)
for bidx in range(1, meta['count'] + 1):
- data = out.read_band(bidx, masked=True)
+ data = out.read(bidx, masked=True)
# Burn in any non-fill pixels, and update mask accordingly
ne = result != fill
data[ne] = result[ne]
@@ -627,7 +610,7 @@ def rasterize(
@click.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), required=True)
@precision_opt
@indent_opt
@compact_opt
@@ -678,22 +661,20 @@ def bounds(ctx, input, precision, indent, compact, projection, dst_crs,
for i, path in enumerate(input):
with rasterio.open(path) as src:
bounds = src.bounds
- xs = [bounds[0], bounds[2]]
- ys = [bounds[1], bounds[3]]
if dst_crs:
- xs, ys = rasterio.warp.transform(
- src.crs, {'init': dst_crs}, xs, ys)
+ bbox = transform_bounds(src.crs,
+ dst_crs, *bounds)
elif projection == 'mercator':
- xs, ys = rasterio.warp.transform(
- src.crs, {'init': 'epsg:3857'}, xs, ys)
+ bbox = transform_bounds(src.crs,
+ {'init': 'epsg:3857'}, *bounds)
elif projection == 'geographic':
- xs, ys = rasterio.warp.transform(
- src.crs, {'init': 'epsg:4326'}, xs, ys)
+ bbox = transform_bounds(src.crs,
+ {'init': 'epsg:4326'}, *bounds)
+ else:
+ bbox = bounds
if precision >= 0:
- xs = [round(v, precision) for v in xs]
- ys = [round(v, precision) for v in ys]
- bbox = [min(xs), min(ys), max(xs), max(ys)]
+ bbox = [round(b, precision) for b in bbox]
yield {
'type': 'Feature',
@@ -701,11 +682,11 @@ def bounds(ctx, input, precision, indent, compact, projection, dst_crs,
'geometry': {
'type': 'Polygon',
'coordinates': [[
- [xs[0], ys[0]],
- [xs[1], ys[0]],
- [xs[1], ys[1]],
- [xs[0], ys[1]],
- [xs[0], ys[0]] ]]},
+ [bbox[0], bbox[1]],
+ [bbox[2], bbox[1]],
+ [bbox[2], bbox[3]],
+ [bbox[0], bbox[3]],
+ [bbox[0], bbox[1]]]]},
'properties': {
'id': str(i),
'title': path,
diff --git a/rasterio/rio/helpers.py b/rasterio/rio/helpers.py
index 1300b4e..b9c9b02 100644
--- a/rasterio/rio/helpers.py
+++ b/rasterio/rio/helpers.py
@@ -4,6 +4,9 @@ Helper objects used by multiple CLI commands.
import json
+import os
+
+from rasterio.errors import FileOverwriteError
def coords(obj):
@@ -64,11 +67,29 @@ def write_features(
fobj.write('\n')
-def resolve_inout(input=None, output=None, files=None):
+def resolve_inout(input=None, output=None, files=None, force_overwrite=False):
"""Resolves inputs and outputs from standard args and options.
- Returns `output_filename, [input_filename0, ...]`."""
+ :param input: a single input filename, optional.
+ :param output: a single output filename, optional.
+ :param files: a sequence of filenames in which the last is the
+ output filename.
+ :param force_overwrite: whether to force overwriting the output
+ file, bool.
+ :return: the resolved output filename and input filenames as a
+ tuple of length 2.
+
+ If provided, the :param:`output` file may be overwritten. An output
+ file extracted from :param:`files` will not be overwritten unless
+ :param:`force_overwrite` is `True`.
+ """
resolved_output = output or (files[-1] if files else None)
+ force_overwrite = output is not None or force_overwrite
+ if not force_overwrite and resolved_output and os.path.exists(
+ resolved_output):
+ raise FileOverwriteError(
+ "file exists and won't be overwritten without use of the "
+ "`-f` or `-o` options.")
resolved_inputs = (
[input] if input else [] +
list(files[:-1 if not output else None]) if files else [])
@@ -76,4 +97,5 @@ def resolve_inout(input=None, output=None, files=None):
def to_lower(ctx, param, value):
+ """Click callback, converts values to lowercase."""
return value.lower()
diff --git a/rasterio/rio/info.py b/rasterio/rio/info.py
index fa390dd..e7acea9 100644
--- a/rasterio/rio/info.py
+++ b/rasterio/rio/info.py
@@ -17,15 +17,6 @@ from rasterio.transform import guard_transform
# Handlers for info module options.
-def from_like_context(ctx, param, value):
- """Return the value for an option from the context if the option
- or `--all` is given, else return None."""
- if ctx.obj and ctx.obj.get('like') and (
- value == 'like' or ctx.obj.get('all_like')):
- return ctx.obj['like'][param.name]
- else:
- return None
-
def all_handler(ctx, param, value):
"""Get tags from a template file or command line."""
@@ -37,7 +28,7 @@ def all_handler(ctx, param, value):
def crs_handler(ctx, param, value):
"""Get crs value from a template file or command line."""
- retval = from_like_context(ctx, param, value)
+ retval = options.from_like_context(ctx, param, value)
if retval is None and value:
try:
retval = json.loads(value)
@@ -50,35 +41,9 @@ def crs_handler(ctx, param, value):
return retval
-def like_handler(ctx, param, value):
- """Copy a dataset's meta property to the command context for access
- from other callbacks."""
- if ctx.obj is None:
- ctx.obj = {}
- if value:
- with rasterio.open(value) as src:
- metadata = src.meta
- ctx.obj['like'] = metadata
- ctx.obj['like']['transform'] = metadata['affine']
- ctx.obj['like']['tags'] = src.tags()
-
-
-def nodata_handler(ctx, param, value):
- """Get nodata value from a template file or command line."""
- retval = from_like_context(ctx, param, value)
- if retval is None and value is not None:
- try:
- retval = float(value)
- except:
- raise click.BadParameter(
- "%s is not a number." % repr(value),
- param=param, param_hint='nodata')
- return retval
-
-
def tags_handler(ctx, param, value):
"""Get tags from a template file or command line."""
- retval = from_like_context(ctx, param, value)
+ retval = options.from_like_context(ctx, param, value)
if retval is None and value:
try:
retval = dict(p.split('=') for p in value)
@@ -91,7 +56,7 @@ def tags_handler(ctx, param, value):
def transform_handler(ctx, param, value):
"""Get transform value from a template file or command line."""
- retval = from_like_context(ctx, param, value)
+ retval = options.from_like_context(ctx, param, value)
if retval is None and value:
try:
value = json.loads(value)
@@ -111,8 +76,7 @@ def transform_handler(ctx, param, value):
@click.command('edit-info', short_help="Edit dataset metadata.")
@options.file_in_arg
- at click.option('--nodata', callback=nodata_handler, default=None,
- help="New nodata value")
+ at options.nodata_opt
@click.option('--crs', callback=crs_handler, default=None,
help="New coordinate reference system")
@click.option('--transform', callback=transform_handler,
@@ -122,13 +86,7 @@ def transform_handler(ctx, param, value):
@click.option('--all', 'allmd', callback=all_handler, flag_value='like',
is_eager=True, default=False,
help="Copy all metadata items from the template file.")
- at click.option(
- '--like',
- type=click.Path(exists=True),
- callback=like_handler,
- is_eager=True,
- help="Raster dataset to use as a template for obtaining affine "
- "transform (bounds and resolution), crs, and nodata values.")
+ at options.like_opt
@click.pass_context
def edit(ctx, input, nodata, crs, transform, tags, allmd, like):
"""Edit a dataset's metadata: coordinate reference system, affine
@@ -287,7 +245,8 @@ def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx,
proj4 = proj4.split('=')[1].upper()
info['crs'] = proj4
info['res'] = src.res
- info['lnglat'] = src.lnglat()
+ if proj4 != '':
+ info['lnglat'] = src.lnglat()
if verbose:
stats = [{'min': float(b.min()),
'max': float(b.max()),
diff --git a/rasterio/rio/merge.py b/rasterio/rio/merge.py
index 6e37193..ee81dec 100644
--- a/rasterio/rio/merge.py
+++ b/rasterio/rio/merge.py
@@ -19,12 +19,8 @@ from rasterio.transform import Affine
@format_opt
@options.bounds_opt
@options.resolution_opt
- at click.option('--nodata', type=float, default=None,
- help="Override nodata values defined in input datasets")
- at click.option('--force-overwrite', '-f', 'force_overwrite', is_flag=True,
- type=bool, default=False,
- help="Do not prompt for confirmation before overwriting output "
- "file")
+ at options.nodata_opt
+ at options.force_overwrite_opt
@click.option('--precision', type=int, default=7,
help="Number of decimal places of precision in alignment of "
"pixels")
@@ -46,6 +42,8 @@ def merge(ctx, files, output, driver, bounds, res, nodata, force_overwrite,
and are otherwise taken from the first input file.
Note: --res changed from 2 parameters in 0.25.
+
+ \b
--res 0.1 0.1 => --res 0.1 (square)
--res 0.1 0.2 => --res 0.1 --res 0.2 (rectangular)
"""
@@ -55,12 +53,8 @@ def merge(ctx, files, output, driver, bounds, res, nodata, force_overwrite,
verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
logger = logging.getLogger('rio')
- output, files = resolve_inout(files=files, output=output)
-
- if os.path.exists(output) and not force_overwrite:
- raise click.ClickException(
- "Output exists and won't be overwritten without the "
- "`-f` option")
+ output, files = resolve_inout(
+ files=files, output=output, force_overwrite=force_overwrite)
sources = [rasterio.open(f) for f in files]
dest, output_transform = merge_tool(sources, bounds=bounds, res=res,
diff --git a/rasterio/rio/options.py b/rasterio/rio/options.py
index b5da2ae..1911c2d 100644
--- a/rasterio/rio/options.py
+++ b/rasterio/rio/options.py
@@ -50,6 +50,7 @@ import os.path
import click
+import rasterio
from rasterio.vfs import parse_path
@@ -101,6 +102,42 @@ def file_in_handler(ctx, param, value):
return path
+def from_like_context(ctx, param, value):
+ """Return the value for an option from the context if the option
+ or `--all` is given, else return None."""
+ if ctx.obj and ctx.obj.get('like') and (
+ value == 'like' or ctx.obj.get('all_like')):
+ return ctx.obj['like'][param.name]
+ else:
+ return None
+
+
+def like_handler(ctx, param, value):
+ """Copy a dataset's meta property to the command context for access
+ from other callbacks."""
+ if ctx.obj is None:
+ ctx.obj = {}
+ if value:
+ with rasterio.open(value) as src:
+ metadata = src.meta
+ ctx.obj['like'] = metadata
+ ctx.obj['like']['transform'] = metadata['affine']
+ ctx.obj['like']['tags'] = src.tags()
+
+
+def nodata_handler(ctx, param, value):
+ """Get nodata value from a template file or command line."""
+ retval = from_like_context(ctx, param, value)
+ if retval is None and value is not None:
+ try:
+ retval = float(value)
+ except:
+ raise click.BadParameter(
+ "%s is not a number." % repr(value),
+ param=param, param_hint='nodata')
+ return retval
+
+
# Singular input file
file_in_arg = click.argument('INPUT', callback=file_in_handler)
@@ -113,7 +150,7 @@ bidx_opt = click.option(
'-b', '--bidx',
type=int,
default=1,
- help="Input file band index (default: 1)")
+ help="Input file band index (default: 1).")
bidx_mult_opt = click.option(
'-b', '--bidx',
@@ -156,8 +193,9 @@ 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).")
+ help="Path to output file (optional alternative to a positional arg). "
+ "Existing files will be overwritten (`--force-overwrite` is "
+ "implied).")
resolution_opt = click.option(
'-r', '--res',
@@ -165,7 +203,7 @@ resolution_opt = click.option(
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')
+ '--res pixel_width --res pixel_height.')
creation_options = click.option(
'--co', 'creation_options',
@@ -181,3 +219,20 @@ rgb_opt = click.option(
flag_value='rgb',
default=False,
help="Set RGB photometric interpretation.")
+
+force_overwrite_opt = click.option(
+ '--force-overwrite', 'force_overwrite',
+ is_flag=True, type=bool, default=False,
+ help="Always overwrite an existing output file.")
+
+nodata_opt = click.option(
+ '--nodata', callback=nodata_handler, default=None,
+ help="New nodata value.")
+
+like_opt = click.option(
+ '--like',
+ type=click.Path(exists=True),
+ callback=like_handler,
+ is_eager=True,
+ help="Raster dataset to use as a template for obtaining affine "
+ "transform (bounds and resolution), crs, and nodata values.")
diff --git a/rasterio/rio/warp.py b/rasterio/rio/warp.py
index 80a060d..b65f8e8 100644
--- a/rasterio/rio/warp.py
+++ b/rasterio/rio/warp.py
@@ -1,66 +1,96 @@
import logging
from math import ceil
+import warnings
+
import click
-from cligj import format_opt
+from cligj import files_inout_arg, format_opt
+from .helpers import resolve_inout
from . import options
import rasterio
from rasterio import crs
from rasterio.transform import Affine
-from rasterio.warp import (reproject, RESAMPLING, calculate_default_transform,
+from rasterio.warp import (reproject, Resampling, calculate_default_transform,
transform_bounds)
-logger = logging.getLogger('rio')
+# Improper usage of rio-warp can lead to accidental creation of
+# extremely large datasets. We'll put a hard limit on the size of
+# datasets and raise a usage error if the limits are exceeded.
+MAX_OUTPUT_WIDTH = 100000
+MAX_OUTPUT_HEIGHT = 100000
+
+
+def bounds_handler(ctx, param, value):
+ """Warn about future usage changes."""
+ if value:
+ click.echo("Future Warning: "
+ "the semantics of the `--bounds` option will change in Rasterio "
+ "version 1.0 from bounds of the source dataset to bounds of the "
+ "destination dataset.", err=True)
+ return value
+
+
+def x_dst_bounds_handler(ctx, param, value):
+ """Warn about future usage changes."""
+ if value:
+ click.echo("Future Warning: "
+ "the `--x-dst-bounds` option will be removed in Rasterio version "
+ "1.0 in favor of `--bounds`.", err=True)
+ return value
@click.command(short_help='Warp a raster dataset.')
- at options.file_in_arg
- at options.file_out_arg
+ at files_inout_arg
+ at options.output_opt
@format_opt
@options.like_file_opt
@click.option('--dst-crs', default=None,
help='Target coordinate reference system.')
@options.dimensions_opt
- at options.bounds_opt
+ at click.option(
+ '--src-bounds',
+ nargs=4, type=float, default=None,
+ help="Determine output extent from source bounds: left bottom right top "
+ "(note: for future backwards compatibility in 1.0).")
+ at click.option(
+ '--x-dst-bounds',
+ nargs=4, type=float, default=None, callback=x_dst_bounds_handler,
+ help="Set output extent from bounding values: left bottom right top "
+ "(note: this option will be removed in 1.0).")
+ at click.option(
+ '--bounds',
+ nargs=4, type=float, default=None, callback=bounds_handler,
+ help="Determine output extent from source bounds: left bottom right top "
+ "(note: the semantics of this option will change to those of "
+ "`--x-dst-bounds` in version 1.0).")
@options.resolution_opt
- at click.option('--resampling', type=click.Choice(['nearest', 'bilinear', 'cubic',
- 'cubic_spline','lanczos', 'average', 'mode']),
- default='nearest', help='Resampling method (default: nearest).')
+ at click.option('--resampling', type=click.Choice([r.name for r in Resampling]),
+ default='nearest', help="Resampling method.",
+ show_default=True)
@click.option('--threads', type=int, default=1,
help='Number of processing threads.')
+ at click.option('--check-invert-proj', type=bool, default=True,
+ help='Constrain output to valid coordinate region in dst-crs')
+ at options.force_overwrite_opt
@options.creation_options
@click.pass_context
-# TODO: add NODATA options and support for existing output rasters
-def warp(
- ctx,
- input,
- output,
- driver,
- like,
- dst_crs,
- dimensions,
- bounds,
- res,
- resampling,
- threads,
- creation_options):
+def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds,
+ x_dst_bounds, bounds, res, resampling, threads, check_invert_proj,
+ force_overwrite, creation_options):
"""
Warp a raster dataset.
- Currently, the output is always overwritten. This will be changed in a
- later version.
-
- If a template raster is provided using the --like option, the coordinate
- reference system, affine transform, and dimensions of that raster will
- be used for the output. In this case --dst-crs, --bounds, --res, and
- --dimensions options are ignored.
+ If a template raster is provided using the --like option, the
+ coordinate reference system, affine transform, and dimensions of
+ that raster will be used for the output. In this case --dst-crs,
+ --bounds, --res, and --dimensions options are ignored.
\b
$ rio warp input.tif output.tif --like template.tif
- The output coordinate reference system may be either a PROJ.4 or EPSG:nnnn
- string,
+ The output coordinate reference system may be either a PROJ.4 or
+ EPSG:nnnn string,
\b
--dst-crs EPSG:4326
@@ -71,25 +101,31 @@ def warp(
\b
--dst-crs '{"proj": "utm", "zone": 18, ...}'
- If --dimensions are provided, --res and --bounds are ignored. Resolution
- is calculated based on the relationship between the raster bounds in the
- target coordinate system and the dimensions, and may produce rectangular
- rather than square pixels.
+ If --dimensions are provided, --res and --bounds are ignored.
+ Resolution is calculated based on the relationship between the
+ raster bounds in the target coordinate system and the dimensions,
+ and may produce rectangular rather than square pixels.
\b
- $ rio warp input.tif output.tif --dimensions 100 200 --dst-crs EPSG:4326
+ $ rio warp input.tif output.tif --dimensions 100 200 \\
+ > --dst-crs EPSG:4326
If --bounds are provided, --res is required if --dst-crs is provided
- (defaults to source raster resolution otherwise). Bounds are in the source
- coordinate reference system.
+ (defaults to source raster resolution otherwise).
\b
- $ rio warp input.tif output.tif --bounds -78 22 -76 24 --dst-crs \\
- EPSG:4326 --res 0.1
+ $ rio warp input.tif output.tif \\
+ > --bounds -78 22 -76 24 --res 0.1 --dst-crs EPSG:4326
+
"""
verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
- resampling = getattr(RESAMPLING, resampling) # get integer code for method
+ logger = logging.getLogger('rio')
+
+ output, files = resolve_inout(
+ files=files, output=output, force_overwrite=force_overwrite)
+
+ resampling = Resampling[resampling] # get integer code for method
if not len(res):
# Click sets this as an empty tuple if not provided
@@ -98,12 +134,21 @@ def warp(
# Expand one value to two if needed
res = (res[0], res[0]) if len(res) == 1 else res
- with rasterio.drivers(CPL_DEBUG=verbosity > 2):
- with rasterio.open(input) as src:
+ with rasterio.drivers(CPL_DEBUG=verbosity > 2,
+ CHECK_WITH_INVERT_PROJ=check_invert_proj):
+ with rasterio.open(files[0]) as src:
l, b, r, t = src.bounds
out_kwargs = src.meta.copy()
out_kwargs['driver'] = driver
+ # Sort out the bounds options.
+ src_bounds = bounds or src_bounds
+ dst_bounds = x_dst_bounds
+ if src_bounds and dst_bounds:
+ raise click.BadParameter(
+ "Source and destination bounds may not be specified "
+ "simultaneously.")
+
if like:
with rasterio.open(like) as template_ds:
dst_crs = template_ds.crs
@@ -115,11 +160,12 @@ def warp(
try:
dst_crs = crs.from_string(dst_crs)
except ValueError:
- raise click.BadParameter('invalid crs format',
+ raise click.BadParameter("invalid crs format.",
param=dst_crs, param_hint=dst_crs)
if dimensions:
- # Calculate resolution appropriate for dimensions in target
+ # Calculate resolution appropriate for dimensions
+ # in target.
dst_width, dst_height = dimensions
xmin, ymin, xmax, ymax = transform_bounds(src.crs, dst_crs,
*src.bounds)
@@ -130,13 +176,18 @@ def warp(
ymax
)
- elif bounds:
+ elif src_bounds or dst_bounds:
if not res:
- raise click.BadParameter('Required when using --bounds',
+ raise click.BadParameter(
+ "Required when using --bounds.",
param='res', param_hint='res')
- xmin, ymin, xmax, ymax = transform_bounds(src.crs, dst_crs,
- *bounds)
+ if src_bounds:
+ xmin, ymin, xmax, ymax = transform_bounds(
+ src.crs, dst_crs, *src_bounds)
+ else:
+ xmin, ymin, xmax, ymax = dst_bounds
+
dst_transform = Affine(res[0], 0, xmin, 0, -res[1], ymax)
dst_width = max(int(ceil((xmax - xmin) / res[0])), 1)
dst_height = max(int(ceil((ymax - ymin) / res[1])), 1)
@@ -147,7 +198,7 @@ def warp(
resolution=res)
elif dimensions:
- # Same projection, different dimensions, calculate resolution
+ # Same projection, different dimensions, calculate resolution.
dst_crs = src.crs
dst_width, dst_height = dimensions
dst_transform = Affine(
@@ -157,20 +208,20 @@ def warp(
t
)
- elif bounds:
- # Same projection, different dimensions and possibly different
- # resolution
+ elif src_bounds or dst_bounds:
+ # Same projection, different dimensions and possibly
+ # different resolution.
if not res:
res = (src.affine.a, -src.affine.e)
dst_crs = src.crs
- xmin, ymin, xmax, ymax = bounds
+ xmin, ymin, xmax, ymax = (src_bounds or dst_bounds)
dst_transform = Affine(res[0], 0, xmin, 0, -res[1], ymax)
dst_width = max(int(ceil((xmax - xmin) / res[0])), 1)
dst_height = max(int(ceil((ymax - ymin) / res[1])), 1)
elif res:
- # Same projection, different resolution
+ # Same projection, different resolution.
dst_crs = src.crs
dst_transform = Affine(res[0], 0, l, 0, -res[1], t)
dst_width = max(int(ceil((r - l) / res[0])), 1)
@@ -182,6 +233,15 @@ def warp(
dst_width = src.width
dst_height = src.height
+ # When the bounds option is misused, extreme values of
+ # destination width and height may result.
+ if (dst_width < 0 or dst_height < 0 or
+ dst_width > MAX_OUTPUT_WIDTH or
+ dst_height > MAX_OUTPUT_HEIGHT):
+ raise click.BadParameter(
+ "Invalid output dimensions: {0}.".format(
+ (dst_width, dst_height)))
+
out_kwargs.update({
'crs': dst_crs,
'transform': dst_transform,
@@ -205,4 +265,4 @@ def warp(
dst_crs=out_kwargs['crs'],
# dst_nodata=#TODO
resampling=resampling,
- num_threads=threads)
\ No newline at end of file
+ num_threads=threads)
diff --git a/rasterio/tool.py b/rasterio/tool.py
index 7804871..ebd3a04 100644
--- a/rasterio/tool.py
+++ b/rasterio/tool.py
@@ -1,12 +1,29 @@
+"""
+Implementations of various common operations, like `show()` for displaying an
+array or with matplotlib, and `stats()` for computing min/max/avg. Most can
+handle a numpy array or `rasterio.Band()`. Primarily supports `$ rio insp`.
+"""
+
+
+from __future__ import absolute_import
import code
import collections
import logging
+import warnings
try:
import matplotlib.pyplot as plt
except ImportError:
plt = None
+except RuntimeError as e:
+ # Certain environment configurations can trigger a RuntimeError like:
+
+ # Trying to import matplotlibRuntimeError: Python is not installed as a
+ # framework. The Mac OS X backend will not be able to function correctly
+ # if Python is not installed as a framework. See the Python ...
+ warnings.warn(str(e), RuntimeWarning, stacklevel=2)
+ plt = None
import numpy
@@ -22,19 +39,41 @@ Stats = collections.namedtuple('Stats', ['min', 'max', 'mean'])
funcs = locals()
-def show(source, cmap='gray'):
- """Show a raster using matplotlib.
+def show(source, cmap='gray', with_bounds=True):
+ """
+ Display a raster or raster band using matplotlib.
- The raster may be either an ndarray or a (dataset, bidx)
- tuple.
+ Parameters
+ ----------
+ source : array-like or (raster dataset, bidx)
+ If array-like, should be of format compatible with
+ matplotlib.pyplot.imshow. If the tuple (raster dataset, bidx),
+ selects band `bidx` from raster.
+ cmap : str (opt)
+ Specifies the colormap to use in plotting. See
+ matplotlib.Colors.Colormap. Default is 'gray'.
+ with_bounds : bool (opt)
+ Whether to change the image extent to the spatial bounds of the image,
+ rather than pixel coordinates. Only works when source is
+ (raster dataset, bidx).
"""
+
if isinstance(source, tuple):
arr = source[0].read(source[1])
+ xs = source[0].res[0] / 2.
+ ys = source[0].res[1] / 2.
+ if with_bounds:
+ extent = (source[0].bounds.left - xs, source[0].bounds.right - xs,
+ source[0].bounds.bottom - ys, source[0].bounds.top - ys)
+ else:
+ extent = None
else:
arr = source
+ extent = None
if plt is not None:
- plt.imshow(arr, cmap=cmap)
- plt.show()
+ imax = plt.imshow(arr, cmap=cmap, extent=extent)
+ fig = plt.gcf()
+ fig.show()
else:
raise ImportError("matplotlib could not be imported")
@@ -116,7 +155,8 @@ def show_hist(source, bins=10, masked=True, title='Histogram'):
plt.grid(True)
plt.xlabel('DN')
plt.ylabel('Frequency')
- plt.show()
+ fig = plt.gcf()
+ fig.show()
def main(banner, dataset, alt_interpreter=None):
diff --git a/rasterio/tools/mask.py b/rasterio/tools/mask.py
new file mode 100644
index 0000000..ec5423c
--- /dev/null
+++ b/rasterio/tools/mask.py
@@ -0,0 +1,92 @@
+from __future__ import absolute_import
+
+import warnings
+
+import rasterio
+from rasterio.features import geometry_mask
+
+
+def mask(raster, shapes, nodata=None, crop=False, all_touched=False,
+ invert=False):
+ """
+ For all regions in the input raster outside of the regions defined by
+ `shapes`, sets any data present to nodata.
+
+ Parameters
+ ----------
+ raster: rasterio RasterReader object
+ Raster to which the mask will be applied.
+ shapes: list of polygons
+ Polygons are GeoJSON-like dicts specifying the boundaries of features
+ in the raster to be kept. All data outside of specified polygons
+ will be set to nodata.
+ nodata: int or float (opt)
+ Value representing nodata within each raster band. If not set,
+ defaults to the nodata value for the input raster. If there is no
+ set nodata value for the raster, it defaults to 0.
+ crop: bool (opt)
+ Whether to crop the raster to the extent of the data. Defaults to
+ False.
+ all_touched: bool (opt)
+ Use all pixels touched by features. If False (default), use only
+ pixels whose center is within the polygon or that are selected by
+ Bresenhams line algorithm.
+ invert: bool (opt)
+ If True, mask will be True for pixels that overlap shapes.
+ False by default.
+
+ Returns
+ -------
+ masked: numpy ndarray
+ Data contained in raster after applying the mask.
+ out_transform: affine object
+ Information for mapping pixel coordinates in `masked` to another
+ coordinate system.
+ """
+
+ if crop and invert:
+ raise ValueError("crop and invert cannot both be True.")
+ if nodata is None:
+ if raster.nodata is not None:
+ nodata = raster.nodata
+ else:
+ nodata = 0
+
+ all_bounds = [rasterio.features.bounds(shape) for shape in shapes]
+ minxs, minys, maxxs, maxys = zip(*all_bounds)
+ mask_bounds = (min(minxs), min(minys), max(maxxs), max(maxys))
+
+ invert_y = raster.affine.e > 0
+ source_bounds = raster.bounds
+ if invert_y:
+ source_bounds = [source_bounds[0], source_bounds[3],
+ source_bounds[2], source_bounds[1]]
+ if rasterio.coords.disjoint_bounds(source_bounds, mask_bounds):
+ if crop:
+ raise ValueError("Input shapes do not overlap raster.")
+ else:
+ warnings.warn("GeoJSON outside bounds of existing output " +
+ "raster. Are they in different coordinate " +
+ "reference systems?")
+ if invert_y:
+ mask_bounds = [mask_bounds[0], mask_bounds[3],
+ mask_bounds[2], mask_bounds[1]]
+ if crop:
+ window = raster.window(*mask_bounds)
+ out_transform = raster.window_transform(window)
+ else:
+ window = None
+ out_transform = raster.affine
+
+ out_image = raster.read(window=window, masked=True)
+ out_shape = out_image.shape[1:]
+
+ shape_mask = geometry_mask(shapes, transform=out_transform, invert=invert,
+ out_shape=out_shape, all_touched=all_touched)
+ out_image.mask = out_image.mask | shape_mask
+ out_image.fill_value = nodata
+
+ for i in range(raster.count):
+ out_image[i] = out_image[i].filled(nodata)
+
+ return out_image, out_transform
diff --git a/rasterio/transform.py b/rasterio/transform.py
index 3bb073f..46db1b1 100644
--- a/rasterio/transform.py
+++ b/rasterio/transform.py
@@ -1,4 +1,5 @@
from __future__ import absolute_import
+from __future__ import division
import warnings
diff --git a/rasterio/warp.py b/rasterio/warp.py
index 34d384f..ed9bdc7 100644
--- a/rasterio/warp.py
+++ b/rasterio/warp.py
@@ -1,14 +1,26 @@
"""Raster warping and reprojection"""
-from affine import Affine
+from __future__ import absolute_import
+from __future__ import division
+
from math import ceil
+import warnings
+
+from affine import Affine
import numpy as np
+import rasterio
from rasterio._base import _transform
-from rasterio._warp import _transform_geom, _reproject, RESAMPLING
+from rasterio._warp import (_transform_geom, _reproject, Resampling,
+ _calculate_default_transform)
from rasterio.transform import guard_transform
+RESAMPLING = Resampling
+warnings.warn(
+ "RESAMPLING is deprecated, use Resampling instead.", DeprecationWarning)
+
+
def transform(src_crs, dst_crs, xs, ys, zs=None):
"""
Transform vectors of x, y and optionally z from source
@@ -150,7 +162,7 @@ def reproject(
dst_transform=None,
dst_crs=None,
dst_nodata=None,
- resampling=RESAMPLING.nearest,
+ resampling=Resampling.nearest,
**kwargs):
"""
Reproject a source raster to a destination raster.
@@ -195,13 +207,13 @@ def reproject(
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
+ Resampling.nearest,
+ Resampling.bilinear,
+ Resampling.cubic,
+ Resampling.cubic_spline,
+ Resampling.lanczos,
+ Resampling.average,
+ Resampling.mode
kwargs: dict, optional
Additional arguments passed to transformation function.
@@ -211,6 +223,9 @@ def reproject(
Output is written to destination.
"""
+ # Resampling guard.
+ _ = Resampling(resampling)
+
if src_transform:
src_transform = guard_transform(src_transform).to_gdal()
if dst_transform:
@@ -238,8 +253,7 @@ def calculate_default_transform(
bottom,
right,
top,
- resolution=None,
- densify_pts=21):
+ resolution=None):
"""
Transforms bounds to destination coordinate system, calculates resolution
if not provided, and returns destination transform and dimensions.
@@ -247,13 +261,8 @@ def calculate_default_transform(
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).
-
+ Destination width and height (and resolution if not provided), are
+ calculated using GDAL's method for suggest warp output.
Parameters
----------
@@ -270,36 +279,50 @@ def calculate_default_transform(
Bounding coordinates in src_crs, from the bounds property of a raster.
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(
- src_crs, dst_crs, left, bottom, right, top, densify_pts)
- x_dif = xmax - xmin
- y_dif = ymax - ymin
- size = float(width + height)
+ Note
+ ----
+ Must be called within a raster.drivers() context
- 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
+ Some behavior of this function is determined by the
+ CHECK_WITH_INVERT_PROJ environment variable
+ YES: constrain output raster to extents that can be inverted
+ avoids visual artifacts and coordinate discontinuties.
+ NO: reproject coordinates beyond valid bound limits
+ """
+ dst_affine, dst_width, dst_height = _calculate_default_transform(
+ src_crs, dst_crs,
+ width, height,
+ left, bottom, right, top)
+
+ # If resolution is specified, Keep upper-left anchored
+ # adjust the transform resolutions
+ # adjust the width/height by the ratio of estimated:specified res (ceil'd)
+ if resolution:
+
+ # resolutions argument into tuple
+ try:
+ res = (float(resolution), float(resolution))
+ except TypeError:
+ res = (resolution[0], resolution[0]) \
+ if len(resolution) == 1 else resolution[0:2]
+
+ # Assume yres is provided as positive,
+ # needs to be negative for north-up affine
+ xres = res[0]
+ yres = -res[1]
+
+ xratio = dst_affine.a / xres
+ yratio = dst_affine.e / yres
+
+ dst_affine = Affine(xres, dst_affine.b, dst_affine.c,
+ dst_affine.d, yres, dst_affine.f)
+
+ dst_width = ceil(dst_width * xratio)
+ dst_height = ceil(dst_height * yratio)
+
+ return dst_affine, dst_width, dst_height
diff --git a/tests/conftest.py b/tests/conftest.py
index 16f39e7..e7e7154 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -169,6 +169,42 @@ def diagonal_image():
@pytest.fixture()
+def basic_image_file(tmpdir, basic_image):
+ """
+ A basic raster file with a 10x10 array for testing sieve functions.
+ Contains data from pixelated_image.
+
+ Returns
+ -------
+
+ string
+ Filename of test raster file
+ """
+
+ from affine import Affine
+ import rasterio
+
+ image = basic_image
+
+ outfilename = str(tmpdir.join('basic_image.tif'))
+ kwargs = {
+ "crs": {'init': 'epsg:4326'},
+ "transform": Affine.identity(),
+ "count": 1,
+ "dtype": rasterio.uint8,
+ "driver": "GTiff",
+ "width": image.shape[1],
+ "height": image.shape[0],
+ "nodata": None
+ }
+ with rasterio.drivers():
+ with rasterio.open(outfilename, 'w', **kwargs) as out:
+ out.write_band(1, image)
+
+ return outfilename
+
+
+ at pytest.fixture()
def pixelated_image_file(tmpdir, pixelated_image):
"""
A basic raster file with a 10x10 array for testing sieve functions.
diff --git a/tests/data/no_band.h5 b/tests/data/no_band.h5
new file mode 100644
index 0000000..e1ce0ec
Binary files /dev/null and b/tests/data/no_band.h5 differ
diff --git a/tests/data/world.byte.tif b/tests/data/world.byte.tif
new file mode 100644
index 0000000..17f420b
Binary files /dev/null and b/tests/data/world.byte.tif differ
diff --git a/tests/test_crs.py b/tests/test_crs.py
index 07b0e0c..30548ad 100644
--- a/tests/test_crs.py
+++ b/tests/test_crs.py
@@ -41,24 +41,8 @@ def test_write_3857(tmpdir):
assert dst.crs == {'init': 'epsg:3857'}
info = subprocess.check_output([
'gdalinfo', dst_path])
- assert """PROJCS["WGS 84 / Pseudo-Mercator",
- GEOGCS["WGS 84",
- DATUM["WGS_1984",
- SPHEROID["WGS 84",6378137,298.257223563,
- AUTHORITY["EPSG","7030"]],
- AUTHORITY["EPSG","6326"]],
- PRIMEM["Greenwich",0],
- UNIT["degree",0.0174532925199433],
- AUTHORITY["EPSG","4326"]],
- PROJECTION["Mercator_1SP"],
- PARAMETER["central_meridian",0],
- PARAMETER["scale_factor",1],
- PARAMETER["false_easting",0],
- PARAMETER["false_northing",0],
- UNIT["metre",1,
- AUTHORITY["EPSG","9001"]],
- EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
- AUTHORITY["EPSG","3857"]]""" in info.decode('utf-8')
+ # WKT string may vary a bit w.r.t GDAL versions
+ assert 'PROJCS["WGS 84 / Pseudo-Mercator"' in info.decode('utf-8')
def test_from_proj4_json():
diff --git a/tests/test_read.py b/tests/test_read.py
index 572c01d..3bbc7d2 100644
--- a/tests/test_read.py
+++ b/tests/test_read.py
@@ -252,3 +252,10 @@ class ReaderContextTest(unittest.TestCase):
self.assertEqual(a.ndim, 3)
self.assertEqual(a.shape, (1, 2, 2))
self.assertTrue(hasattr(a, 'mask'))
+
+ def test_read_no_band(self):
+ with rasterio.open('tests/data/no_band.h5') as s:
+ self.assertEqual(s.count, 0)
+ self.assertEqual(s.meta['dtype'], 'float_')
+ self.assertIsNone(s.meta['nodata'])
+ self.assertRaises(ValueError, s.read)
diff --git a/tests/test_rio_features.py b/tests/test_rio_features.py
index b1f20a0..93c43b5 100644
--- a/tests/test_rio_features.py
+++ b/tests/test_rio_features.py
@@ -608,7 +608,7 @@ def test_rasterize_existing_output(tmpdir, runner, basic_feature):
result = runner.invoke(
features.rasterize,
- [output, '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1]],
+ ['-o', output, '--dimensions', DEFAULT_SHAPE[0], DEFAULT_SHAPE[1]],
input=json.dumps(basic_feature)
)
@@ -733,7 +733,7 @@ def test_rasterize_missing_parameters(tmpdir, runner, basic_feature):
output = str(tmpdir.join('test.tif'))
result = runner.invoke(
- features.rasterize, [output], input=json.dumps(basic_feature)
+ features.rasterize, ['-o', output], input=json.dumps(basic_feature)
)
assert result.exit_code == 2
diff --git a/tests/test_rio_helpers.py b/tests/test_rio_helpers.py
index 5e60740..89e0821 100644
--- a/tests/test_rio_helpers.py
+++ b/tests/test_rio_helpers.py
@@ -1,3 +1,6 @@
+import pytest
+
+from rasterio.errors import FileOverwriteError
from rasterio.rio import helpers
@@ -18,6 +21,30 @@ def test_resolve_files_inout__inout_files_output_o():
files=('a', 'b', 'c'), output='out') == ('out', ['a', 'b', 'c'])
+def test_fail_overwrite(tmpdir):
+ """Unforced overwrite of existing file fails."""
+ foo_tif = tmpdir.join('foo.tif')
+ foo_tif.write("content")
+ with pytest.raises(FileOverwriteError) as excinfo:
+ helpers.resolve_inout(files=[str(x) for x in tmpdir.listdir()])
+ assert "file exists and won't be overwritten without use of the " in str(excinfo.value)
+
+def test_force_overwrite(tmpdir):
+ """Forced overwrite of existing file succeeds."""
+ foo_tif = tmpdir.join('foo.tif')
+ foo_tif.write("content")
+ output, inputs = helpers.resolve_inout(
+ files=[str(x) for x in tmpdir.listdir()], force_overwrite=True)
+ assert output == str(foo_tif)
+
+
+def test_implicit_overwrite(tmpdir):
+ """Implicit overwrite of existing file succeeds."""
+ foo_tif = tmpdir.join('foo.tif')
+ foo_tif.write("content")
+ output, inputs = helpers.resolve_inout(output=str(foo_tif))
+ assert output == str(foo_tif)
+
+
def test_to_lower():
-
- assert helpers.to_lower(None, None, 'EPSG:3857') == 'epsg:3857'
\ No newline at end of file
+ assert helpers.to_lower(None, None, 'EPSG:3857') == 'epsg:3857'
diff --git a/tests/test_rio_info.py b/tests/test_rio_info.py
index eb4e1f3..f8efb44 100644
--- a/tests/test_rio_info.py
+++ b/tests/test_rio_info.py
@@ -127,12 +127,6 @@ class MockOption:
self.name = name
-def test_like_dataset_callback(data):
- ctx = MockContext()
- info.like_handler(ctx, 'like', str(data.join('RGB.byte.tif')))
- assert ctx.obj['like']['crs'] == {'init': 'epsg:32618'}
-
-
def test_all_callback_pass(data):
ctx = MockContext()
ctx.obj['like'] = {'transform': 'foo'}
@@ -170,31 +164,6 @@ def test_transform_callback(data):
assert info.transform_handler(ctx, MockOption('transform'), 'like') == 'foo'
-def test_nodata_callback_err(data):
- ctx = MockContext()
- ctx.obj['like'] = {'nodata': 'lolwut'}
- with pytest.raises(click.BadParameter):
- info.nodata_handler(ctx, MockOption('nodata'), 'lolwut')
-
-
-def test_nodata_callback_pass(data):
- """Always return None if the value is None"""
- ctx = MockContext()
- ctx.obj['like'] = {'nodata': -1}
- assert info.nodata_handler(ctx, MockOption('nodata'), None) is None
-
-
-def test_nodata_callback_0(data):
- ctx = MockContext()
- assert info.nodata_handler(ctx, MockOption('nodata'), '0') == 0.0
-
-
-def test_nodata_callback(data):
- ctx = MockContext()
- ctx.obj['like'] = {'nodata': -1}
- assert info.nodata_handler(ctx, MockOption('nodata'), 'like') == -1.0
-
-
def test_crs_callback_pass(data):
"""Always return None if the value is None"""
ctx = MockContext()
@@ -513,7 +482,7 @@ def test_bounds_obj_bbox():
'--precision', '2'
])
assert result.exit_code == 0
- assert result.output.strip() == '[-78.9, 23.56, -76.6, 25.55]'
+ assert result.output.strip() == '[-78.96, 23.56, -76.57, 25.55]'
def test_bounds_compact():
@@ -526,7 +495,7 @@ def test_bounds_compact():
'--compact'
])
assert result.exit_code == 0
- assert result.output.strip() == '[-78.9,23.56,-76.6,25.55]'
+ assert result.output.strip() == '[-78.96,23.56,-76.57,25.55]'
def test_bounds_indent():
@@ -553,7 +522,7 @@ def test_bounds_obj_bbox_mercator():
])
assert result.exit_code == 0
assert result.output.strip() == (
- '[-8782900.033, 2700489.278, -8527010.472, 2943560.235]')
+ '[-8789636.708, 2700489.278, -8524281.514, 2943560.235]')
def test_bounds_obj_bbox_projected():
@@ -605,7 +574,7 @@ def test_bounds_seq():
])
assert result.exit_code == 0
assert result.output == (
- '[-78.9, 23.56, -76.6, 25.55]\n[-78.9, 23.56, -76.6, 25.55]\n')
+ '[-78.96, 23.56, -76.57, 25.55]\n[-78.96, 23.56, -76.57, 25.55]\n')
assert '\x1e' not in result.output
@@ -622,8 +591,7 @@ def test_bounds_seq_rs():
])
assert result.exit_code == 0
assert result.output == (
- '\x1e[-78.9, 23.56, -76.6, 25.55]\n\x1e[-78.9, 23.56, -76.6, 25.55]\n')
-
+ '\x1e[-78.96, 23.56, -76.57, 25.55]\n\x1e[-78.96, 23.56, -76.57, 25.55]\n')
def test_insp():
runner = CliRunner()
diff --git a/tests/test_rio_merge.py b/tests/test_rio_merge.py
index db29e97..9dfc091 100644
--- a/tests/test_rio_merge.py
+++ b/tests/test_rio_merge.py
@@ -132,7 +132,7 @@ def test_merge_output_exists(tmpdir):
def test_merge_output_exists_without_nodata_fails(test_data_dir_2):
- """Fails without -f or --force-overwrite"""
+ """Fails without --force-overwrite"""
runner = CliRunner()
result = runner.invoke(
merge,
@@ -142,11 +142,11 @@ def test_merge_output_exists_without_nodata_fails(test_data_dir_2):
def test_merge_output_exists_without_nodata(test_data_dir_2):
- """Succeeds with -f"""
+ """Succeeds with --force-overwrite"""
runner = CliRunner()
result = runner.invoke(
merge,
- ['-f', str(test_data_dir_2.join('a.tif')),
+ ['--force-overwrite', str(test_data_dir_2.join('a.tif')),
str(test_data_dir_2.join('b.tif'))])
assert result.exit_code == 0
diff --git a/tests/test_rio_options.py b/tests/test_rio_options.py
index edb7245..43a52a4 100644
--- a/tests/test_rio_options.py
+++ b/tests/test_rio_options.py
@@ -4,7 +4,7 @@ import uuid
import click
import pytest
-from rasterio.rio.options import file_in_handler
+from rasterio.rio.options import file_in_handler, like_handler, nodata_handler
class MockContext:
@@ -61,3 +61,34 @@ def test_file_in_handler_with_vfs_file():
ctx = MockContext()
retval = file_in_handler(ctx, 'INPUT', 'file://tests/data/RGB.byte.tif')
assert retval.endswith('tests/data/RGB.byte.tif')
+
+
+def test_like_dataset_callback(data):
+ ctx = MockContext()
+ like_handler(ctx, 'like', str(data.join('RGB.byte.tif')))
+ assert ctx.obj['like']['crs'] == {'init': 'epsg:32618'}
+
+
+def test_nodata_callback_err(data):
+ ctx = MockContext()
+ ctx.obj['like'] = {'nodata': 'lolwut'}
+ with pytest.raises(click.BadParameter):
+ nodata_handler(ctx, MockOption('nodata'), 'lolwut')
+
+
+def test_nodata_callback_pass(data):
+ """Always return None if the value is None"""
+ ctx = MockContext()
+ ctx.obj['like'] = {'nodata': -1}
+ assert nodata_handler(ctx, MockOption('nodata'), None) is None
+
+
+def test_nodata_callback_0(data):
+ ctx = MockContext()
+ assert nodata_handler(ctx, MockOption('nodata'), '0') == 0.0
+
+
+def test_nodata_callback(data):
+ ctx = MockContext()
+ ctx.obj['like'] = {'nodata': -1}
+ assert nodata_handler(ctx, MockOption('nodata'), 'like') == -1.0
diff --git a/tests/test_rio_warp.py b/tests/test_rio_warp.py
index b5a57fc..8aa779b 100644
--- a/tests/test_rio_warp.py
+++ b/tests/test_rio_warp.py
@@ -2,7 +2,9 @@ import logging
import os
import re
import sys
+
import numpy
+import pytest
import rasterio
from rasterio.rio import warp
@@ -112,11 +114,12 @@ def test_warp_reproject_dst_crs(runner, tmpdir):
with rasterio.open(outputname) as output:
assert output.count == src.count
assert output.crs == {'init': 'epsg:4326'}
- assert output.width == 824
- assert output.height == 686
+ assert output.width == 835
+ assert output.height == 696
assert numpy.allclose(output.bounds,
- [-78.95864996545055, 23.564424693996177,
- -76.57259451863895, 25.550873767433984])
+ [-78.95864996545055, 23.564787976164418,
+ -76.5759177302349, 25.550873767433984])
+
def test_warp_reproject_dst_crs_error(runner, tmpdir):
srcname = 'tests/data/RGB.byte.tif'
@@ -184,14 +187,49 @@ def test_warp_reproject_bounds_no_res(runner, tmpdir):
assert result.exit_code == 2
-def test_warp_reproject_bounds_res(runner, tmpdir):
+def test_warp_reproject_multi_bounds_fail(runner, tmpdir):
+ """Mixing --bounds and --x-dst-bounds fails."""
srcname = 'tests/data/shade.tif'
outputname = str(tmpdir.join('test.tif'))
out_bounds = [-11850000, 4810000, -11849000, 4812000]
result = runner.invoke(warp.warp, [srcname, outputname,
'--dst-crs', 'EPSG:4326',
- '--res', 0.001, '--bounds', ]
+ '--x-dst-bounds'] + out_bounds +
+ ['--bounds'] + out_bounds)
+ assert result.exit_code == 2
+
+
+def test_warp_reproject_bounds_crossup_fail(runner, tmpdir):
+ """Crossed-up bounds raises click.BadParameter."""
+ srcname = 'tests/data/shade.tif'
+ outputname = str(tmpdir.join('test.tif'))
+ out_bounds = [-11850000, 4810000, -11849000, 4812000]
+ result = runner.invoke(warp.warp, [srcname, outputname,
+ '--dst-crs', 'EPSG:4326',
+ '--res', 0.001, '--x-dst-bounds', ]
+ out_bounds)
+ assert result.exit_code == 2
+
+
+def test_warp_reproject_bounds_res_future_warning(runner, tmpdir):
+ """Use of --bounds results in a warning from the 1.0 future."""
+ srcname = 'tests/data/shade.tif'
+ outputname = str(tmpdir.join('test.tif'))
+ out_bounds = [-11850000, 4810000, -11849000, 4812000]
+ result = runner.invoke(
+ warp.warp, [srcname, outputname, '--dst-crs', 'EPSG:4326',
+ '--res', 0.001, '--bounds'] + out_bounds)
+ assert "Future Warning" in result.output
+
+
+def test_warp_reproject_src_bounds_res(runner, tmpdir):
+ """--src-bounds option works."""
+ srcname = 'tests/data/shade.tif'
+ outputname = str(tmpdir.join('test.tif'))
+ out_bounds = [-11850000, 4810000, -11849000, 4812000]
+ result = runner.invoke(
+ warp.warp, [srcname, outputname, '--dst-crs', 'EPSG:4326',
+ '--res', 0.001, '--src-bounds'] + out_bounds)
assert result.exit_code == 0
assert os.path.exists(outputname)
@@ -206,6 +244,34 @@ def test_warp_reproject_bounds_res(runner, tmpdir):
assert output.height == 14
+def test_warp_reproject_dst_bounds(runner, tmpdir):
+ """--x-dst-bounds option works."""
+ srcname = 'tests/data/shade.tif'
+ outputname = str(tmpdir.join('test.tif'))
+ out_bounds = [-106.45036, 39.6138, -106.44136, 39.6278]
+ result = runner.invoke(
+ warp.warp, [srcname, outputname, '--dst-crs', 'EPSG:4326',
+ '--res', 0.001, '--x-dst-bounds'] + out_bounds)
+ assert result.exit_code == 0
+ assert os.path.exists(outputname)
+
+ with rasterio.open(srcname) as src:
+ with rasterio.open(outputname) as output:
+ assert output.crs == {'init': 'epsg:4326'}
+ assert numpy.allclose(output.bounds[0::3],
+ [-106.45036, 39.6278])
+ assert numpy.allclose([0.001, 0.001],
+ [output.affine.a, -output.affine.e])
+
+ # XXX: an extra row and column is produced in the dataset
+ # because we're using ceil instead of floor internally.
+ # Not necessarily a bug, but may change in the future.
+ assert numpy.allclose([output.bounds[2]-0.001, output.bounds[1]+0.001],
+ [-106.44136, 39.6138])
+ assert output.width == 10
+ assert output.height == 15
+
+
def test_warp_reproject_like(runner, tmpdir):
likename = str(tmpdir.join('like.tif'))
kwargs = {
@@ -236,3 +302,18 @@ def test_warp_reproject_like(runner, tmpdir):
assert numpy.allclose([0.001, 0.001], [output.affine.a, -output.affine.e])
assert output.width == 10
assert output.height == 10
+
+
+def test_warp_reproject_nolostdata(runner, tmpdir):
+ srcname = 'tests/data/world.byte.tif'
+ outputname = str(tmpdir.join('test.tif'))
+ result = runner.invoke(warp.warp, [srcname, outputname,
+ '--dst-crs', 'EPSG:3857'])
+ assert result.exit_code == 0
+ assert os.path.exists(outputname)
+
+ with rasterio.open(outputname) as output:
+ arr = output.read()
+ # 50 column swath on the right edge should have some ones (gdalwarped has 7223)
+ assert arr[0, :, -50:].sum() > 7000
+ assert output.crs == {'init': 'epsg:3857'}
diff --git a/tests/test_tool.py b/tests/test_tool.py
index b0e5285..9348686 100644
--- a/tests/test_tool.py
+++ b/tests/test_tool.py
@@ -8,7 +8,6 @@ except ImportError:
import rasterio
from rasterio.tool import show, show_hist, stats
-
def test_stats():
with rasterio.drivers():
with rasterio.open('tests/data/RGB.byte.tif') as src:
@@ -21,25 +20,50 @@ def test_stats():
assert np.allclose(np.array(results), np.array(results2))
-def test_show():
+def test_show_raster():
"""
This test only verifies that code up to the point of plotting with
matplotlib works correctly. Tests do not exercise matplotlib.
"""
- if plt:
- # Return because plotting causes the tests to block until the plot
- # window is closed.
- return
with rasterio.drivers():
with rasterio.open('tests/data/RGB.byte.tif') as src:
try:
show((src, 1))
+ fig = plt.gcf()
+ plt.close(fig)
+ except ImportError:
+ pass
+
+
+def test_show_raster_no_bounds():
+ """
+ This test only verifies that code up to the point of plotting with
+ matplotlib works correctly. Tests do not exercise matplotlib.
+ """
+
+ with rasterio.drivers():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ try:
+ show((src, 1), with_bounds=False)
+ fig = plt.gcf()
+ plt.close(fig)
except ImportError:
pass
+
+def test_show_array():
+ """
+ This test only verifies that code up to the point of plotting with
+ matplotlib works correctly. Tests do not exercise matplotlib.
+ """
+
+ with rasterio.drivers():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
try:
show(src.read(1))
+ fig = plt.gcf()
+ plt.close(fig)
except ImportError:
pass
@@ -49,17 +73,18 @@ def test_show_hist():
This test only verifies that code up to the point of plotting with
matplotlib works correctly. Tests do not exercise matplotlib.
"""
- if plt:
- # Return because plotting causes the tests to block until the plot
- # window is closed.
- return
+
with rasterio.open('tests/data/RGB.byte.tif') as src:
try:
show_hist((src, 1), bins=256)
+ fig = plt.gcf()
+ plt.close(fig)
except ImportError:
pass
try:
show_hist(src.read(), bins=256)
+ fig = plt.gcf()
+ plt.close(fig)
except ImportError:
pass
diff --git a/tests/test_tools_mask.py b/tests/test_tools_mask.py
new file mode 100644
index 0000000..8398f29
--- /dev/null
+++ b/tests/test_tools_mask.py
@@ -0,0 +1,51 @@
+import pytest
+
+import rasterio
+from rasterio.tools.mask import mask as mask_tool
+
+
+def test_nodata(basic_image_file, basic_geometry):
+ nodata_val = 0
+ geometries = [basic_geometry]
+ with rasterio.open(basic_image_file, "r") as src:
+ masked, transform = mask_tool(src, geometries, crop=False,
+ nodata=nodata_val, invert=True)
+ assert(masked.data.all() == nodata_val)
+
+
+def test_no_nodata(basic_image_file, basic_geometry):
+ default_nodata_val = 0
+ geometries = [basic_geometry]
+ with rasterio.open(basic_image_file, "r") as src:
+ masked, transform = mask_tool(src, geometries, crop=False, invert=True)
+ assert(masked.data.all() == default_nodata_val)
+
+
+def test_crop(basic_image, basic_image_file, basic_geometry):
+ geometries = [basic_geometry]
+ with rasterio.open(basic_image_file, "r") as src:
+ masked, transform = mask_tool(src, geometries, crop=True)
+
+ image = basic_image
+ image[4, :] = 0
+ image[:, 4] = 0
+ assert(masked.shape == (1, 4, 3))
+ assert((masked[0] == image[1:5, 2:5]).all())
+
+
+def test_crop_all_touched(basic_image, basic_image_file, basic_geometry):
+ geometries = [basic_geometry]
+ with rasterio.open(basic_image_file, "r") as src:
+ masked, transform = mask_tool(src, geometries, crop=True,
+ all_touched=True)
+
+ assert(masked.shape == (1, 4, 3))
+ assert((masked[0] == basic_image[1:5, 2:5]).all())
+
+
+def test_crop_and_invert(basic_image_file, basic_geometry):
+ geometries = [basic_geometry]
+ with rasterio.open(basic_image_file) as src:
+ with pytest.raises(ValueError):
+ masked, transform = mask_tool(src, geometries,
+ crop=True, invert=True)
diff --git a/tests/test_transform.py b/tests/test_transform.py
index 23a8efe..323b6a2 100644
--- a/tests/test_transform.py
+++ b/tests/test_transform.py
@@ -120,3 +120,24 @@ def test_affine_identity(tmpdir):
with rasterio.open(output) as out:
assert out.affine == Affine.identity()
+
+
+def test_from_bounds_two():
+ width = 80
+ height = 80
+ left = -120
+ top = 70
+ right = -80.5
+ bottom = 30.5
+ tr = transform.from_bounds(left, bottom, right, top, width, height)
+ # pixelwidth, rotation, ULX, rotation, pixelheight, ULY
+ expected = Affine(0.49375, 0.0, -120.0, 0.0, -0.49375, 70.0)
+ assert [round(v, 7) for v in tr] == [round(v, 7) for v in expected]
+
+ # Round right and bottom
+ right = -80
+ bottom = 30
+ tr = transform.from_bounds(left, bottom, right, top, width, height)
+ # pixelwidth, rotation, ULX, rotation, pixelheight, ULY
+ expected = Affine(0.5, 0.0, -120.0, 0.0, -0.5, 70.0)
+ assert [round(v, 7) for v in tr] == [round(v, 7) for v in expected]
diff --git a/tests/test_warp.py b/tests/test_warp.py
index d64dfae..3c1a4b8 100644
--- a/tests/test_warp.py
+++ b/tests/test_warp.py
@@ -7,7 +7,7 @@ import numpy
import rasterio
from rasterio.warp import (
- reproject, RESAMPLING, transform_geom, transform, transform_bounds,
+ reproject, Resampling, transform_geom, transform, transform_bounds,
calculate_default_transform)
@@ -131,9 +131,9 @@ def test_transform_bounds_densify_out_of_bounds():
def test_calculate_default_transform():
target_transform = Affine(
- 0.0028956983577810586, 0.0, -78.95864996545055,
- 0.0, -0.0028956983577810586, 25.550873767433984
- )
+ 0.0028535715391804096, 0.0, -78.95864996545055,
+ 0.0, -0.0028535715391804096, 25.550873767433984)
+
with rasterio.drivers():
with rasterio.open('tests/data/RGB.byte.tif') as src:
wgs84_crs = {'init': 'EPSG:4326'}
@@ -141,14 +141,13 @@ def test_calculate_default_transform():
src.crs, wgs84_crs, src.width, src.height, *src.bounds)
assert dst_transform.almost_equals(target_transform)
- assert width == 824
- assert height == 686
+ assert width == 835
+ assert height == 696
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,
@@ -209,7 +208,7 @@ def test_reproject_ndarray():
src_crs=src.crs,
dst_transform=DST_TRANSFORM,
dst_crs=dst_crs,
- resampling=RESAMPLING.nearest)
+ resampling=Resampling.nearest)
assert (out > 0).sum() == 438146
@@ -227,7 +226,7 @@ def test_reproject_epsg():
src_crs=src.crs,
dst_transform=DST_TRANSFORM,
dst_crs=dst_crs,
- resampling=RESAMPLING.nearest)
+ resampling=Resampling.nearest)
assert (out > 0).sum() == 438146
@@ -246,7 +245,7 @@ def test_reproject_out_of_bounds():
src_crs=src.crs,
dst_transform=DST_TRANSFORM,
dst_crs=dst_crs,
- resampling=RESAMPLING.nearest)
+ resampling=Resampling.nearest)
assert not out.any()
@@ -271,9 +270,36 @@ def test_reproject_nodata():
dst_nodata=nodata
)
- assert (out == 1).sum() == 4461
+ assert (out == 1).sum() == 6215
assert (out == nodata).sum() == (params.dst_width *
- params.dst_height - 4461)
+ params.dst_height - 6215)
+
+
+def test_reproject_nodata_nan():
+ params = default_reproject_params()
+ nodata = 215
+
+ with rasterio.drivers():
+ source = numpy.ones((params.width, params.height), dtype=numpy.float32)
+ 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=numpy.nan,
+ dst_transform=params.dst_transform,
+ dst_crs=params.dst_crs,
+ dst_nodata=numpy.nan
+ )
+
+ assert (out == 1).sum() == 6215
+ assert numpy.isnan(out).sum() == (params.dst_width *
+ params.dst_height - 6215)
+
def test_reproject_dst_nodata_default():
@@ -299,9 +325,9 @@ def test_reproject_dst_nodata_default():
dst_crs=params.dst_crs
)
- assert (out == 1).sum() == 4461
+ assert (out == 1).sum() == 6215
assert (out == 0).sum() == (params.dst_width *
- params.dst_height - 4461)
+ params.dst_height - 6215)
def test_reproject_invalid_dst_nodata():
@@ -392,7 +418,7 @@ def test_reproject_multi():
src_crs=src.crs,
dst_transform=DST_TRANSFORM,
dst_crs=dst_crs,
- resampling=RESAMPLING.nearest)
+ resampling=Resampling.nearest)
assert destin.any()
@@ -527,10 +553,30 @@ def test_transform_geom():
'EPSG:3373',
'EPSG:4326',
geom,
- antimeridian_cutting=True,
+ antimeridian_cutting=True,
antimeridian_offset=0)
assert result['type'] == 'MultiPolygon'
assert len(result['coordinates']) == 2
result = transform_geom('EPSG:3373', 'EPSG:4326', geom, precision=1)
assert int(result['coordinates'][0][0][0] * 10) == -1778
+
+
+def test_reproject_invalid_resampling():
+ # An invalid resampling number fails.
+ 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)
+ with pytest.raises(ValueError) as exc_info:
+ reproject(
+ source,
+ out,
+ src_transform=src.transform,
+ src_crs=src.crs,
+ dst_transform=DST_TRANSFORM,
+ dst_crs=dst_crs,
+ resampling=99)
+ assert "99 is not a valid Resampling" in exc_info.value.message
diff --git a/tests/test_warp_transform.py b/tests/test_warp_transform.py
new file mode 100644
index 0000000..8019df6
--- /dev/null
+++ b/tests/test_warp_transform.py
@@ -0,0 +1,63 @@
+import rasterio
+from rasterio._warp import _calculate_default_transform
+from rasterio.transform import Affine, from_bounds
+
+
+def test_identity():
+ """Get the same transform and dimensions back for same crs."""
+ # Tile: [53, 96, 8]
+ # [-11740727.544603072, 4852834.0517692715, -11584184.510675032, 5009377.085697309]
+ src_crs = dst_crs = 'EPSG:3857'
+ width = height = 1000
+ left, bottom, right, top = (
+ -11740727.544603072, 4852834.0517692715, -11584184.510675032,
+ 5009377.085697309)
+ transform = from_bounds(left, bottom, right, top, width, height)
+
+ with rasterio.drivers():
+ res_transform, res_width, res_height = _calculate_default_transform(
+ src_crs, dst_crs, width, height, left, bottom, right, top)
+
+ assert res_width == width
+ assert res_height == height
+ for res, exp in zip(res_transform, transform):
+ assert round(res, 7) == round(exp, 7)
+
+
+def test_gdal_transform_notnull():
+ with rasterio.drivers():
+ dt, dw, dh = _calculate_default_transform(
+ src_crs={'init': 'EPSG:4326'},
+ dst_crs={'init': 'EPSG:32610'},
+ width=80,
+ height=80,
+ left=-120,
+ bottom=30,
+ right=-80,
+ top=70)
+ assert True
+
+
+def test_gdal_transform_fail_dst_crs():
+ with rasterio.drivers():
+ dt, dw, dh = _calculate_default_transform(
+ {'init': 'EPSG:4326'},
+ '+proj=foobar',
+ width=80,
+ height=80,
+ left=-120,
+ bottom=30,
+ right=-80,
+ top=70)
+
+def test_gdal_transform_fail_src_crs():
+ with rasterio.drivers():
+ dt, dw, dh = _calculate_default_transform(
+ '+proj=foobar',
+ {'init': 'EPSG:32610'},
+ width=80,
+ height=80,
+ left=-120,
+ bottom=30,
+ right=-80,
+ top=70)
--
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