This is an automated email from the git hooks/post-receive script.
johanvdw-guest pushed a commit to branch master
in repository rasterio.
commit 498df05d2b5509e9d77d3ba48f99a19ff66843e4
Author: Johan Van de Wauw <johan.vandew...@gmail.com>
Date: Tue Feb 10 21:37:23 2015 +0100
Imported Upstream version 0.18
---
CHANGES.txt | 13 +
build-wheels.sh | 18 -
docs/cli.rst | 191 +++++++--
docs/reproject.rst | 50 ++-
rasterio/__init__.py | 2 +-
rasterio/_base.pxd | 1 +
rasterio/_base.pyx | 119 +++++-
rasterio/_drivers.pyx | 4 +-
rasterio/_err.pyx | 2 +-
rasterio/_example.pyx | 2 +
rasterio/_features.pyx | 64 ++-
rasterio/_fill.pyx | 84 ++++
rasterio/_gdal.pxd | 4 +-
rasterio/_io.pyx | 48 ++-
rasterio/_warp.pyx | 109 +-----
rasterio/features.py | 23 +-
rasterio/fill.py | 47 +++
rasterio/rasterfill.cpp | 886 ++++++++++++++++++++++++++++++++++++++++++
rasterio/rio/features.py | 236 ++++++++++-
rasterio/rio/info.py | 42 +-
rasterio/rio/main.py | 3 +-
rasterio/rio/sample.py | 94 +++++
rasterio/warp.py | 3 +-
requirements-dev.txt | 2 +-
setup.py | 22 +-
tests/conftest.py | 8 +
tests/test_features_bounds.py | 65 ++++
tests/test_fillnodata.py | 45 +++
tests/test_rio_features.py | 282 +++++++++++++-
tests/test_rio_info.py | 42 ++
tests/test_rio_sample.py | 81 ++++
tests/test_sampling.py | 17 +
tests/test_tags.py | 8 +
tests/test_transform.py | 14 +-
34 files changed, 2398 insertions(+), 233 deletions(-)
diff --git a/CHANGES.txt b/CHANGES.txt
index fc51402..9dce1f6 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,6 +1,19 @@
Changes
=======
+0.18.0 (2015-02-10)
+-------------------
+- New rio-rasterize command (#187).
+- New window_transform method (#215).
+- New sample method and rio-sample command (#251, #275).
+- New fillnodata function based on GDAL's rasterfill.cpp (#253).
+- Speedups for _features and _warp modules (#259).
+- Enhancements for rio-info: 'res', 'lnglat', and 'stats' (#269, #270).
+
+0.17.1 (2015-01-20)
+-------------------
+- Properly handle metadata tags with values that contain "=" (#254).
+
0.17.0 (2015-01-15)
-------------------
- Enhancements to rio-merge: relaxation of same-extent and same-resolution
diff --git a/build-wheels.sh b/build-wheels.sh
deleted file mode 100644
index 36e68dc..0000000
--- a/build-wheels.sh
+++ /dev/null
@@ -1,18 +0,0 @@
-#!/bin/bash
-
-# Automation of this is a TODO. For now, it depends on manually built libraries
-# as detailed in https://gist.github.com/sgillies/a8a2fb910a98a8566d0a.
-
-export MACOSX_DEPLOYMENT_TARGET=10.6
-export GDAL_CONFIG="/usr/local/bin/gdal-config"
-export PACKAGE_DATA=1
-
-VERSION=$1
-
-source $HOME/envs/riowhl27/bin/activate
-CFLAGS="`$GDAL_CONFIG --cflags`" LDFLAGS="`$GDAL_CONFIG --libs` `$GDAL_CONFIG
--dep-libs`" python setup.py bdist_wheel -d wheels/$VERSION
-source $HOME/envs/riowhl34/bin/activate
-CFLAGS="`$GDAL_CONFIG --cflags`" LDFLAGS="`$GDAL_CONFIG --libs` `$GDAL_CONFIG
--dep-libs`" python setup.py bdist_wheel -d wheels/$VERSION
-
-parallel delocate-wheel -w fixed_wheels/$VERSION --require-archs=intel -v {}
::: wheels/$VERSION/*.whl
-parallel cp {}
{.}.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl
::: fixed_wheels/$VERSION/*.whl
diff --git a/docs/cli.rst b/docs/cli.rst
index a744c4a..4e6adf2 100644
--- a/docs/cli.rst
+++ b/docs/cli.rst
@@ -13,19 +13,21 @@ Rasterio's new command line interface is a program named
"rio".
Options:
-v, --verbose Increase verbosity.
-q, --quiet Decrease verbosity.
+ --version Show the version and exit.
--help Show this message and exit.
Commands:
bounds Write bounding boxes to stdout as GeoJSON.
+ env Print information about the rio environment.
info Print information about a data file.
insp Open a data file and start an interpreter.
merge Merge a stack of raster datasets.
+ rasterize Rasterize features.
+ sample Sample a dataset.
shapes Write the shapes of features.
stack Stack a number of bands into a multiband dataset.
- transform Transform coordinates.
-
-It is developed using the ``click`` package.
+It is developed using `Click <http://click.pocoo.org/3/>`__.
bounds
------
@@ -84,26 +86,89 @@ Shoot the GeoJSON into a Leaflet map using geojsonio-cli by
typing
info
----
-Rio's info command intends to serve some of the same uses as gdalinfo.
+Rio's info command prints structured information about a dataset.
.. code-block:: console
- $ rio info tests/data/RGB.byte.tif
- { 'affine': Affine(300.0379266750948, 0.0, 101985.0,
- 0.0, -300.041782729805, 2826915.0),
- 'count': 3,
- 'crs': { 'init': u'epsg:32618'},
- 'driver': u'GTiff',
- 'dtype': <type 'numpy.uint8'>,
- 'height': 718,
- 'nodata': 0.0,
- 'transform': ( 101985.0,
- 300.0379266750948,
- 0.0,
- 2826915.0,
- 0.0,
- -300.041782729805),
- 'width': 791}
+ $ rio info tests/data/RGB.byte.tif --indent 2
+ {
+ "count": 3,
+ "crs": "EPSG:32618",
+ "dtype": "uint8",
+ "driver": "GTiff",
+ "bounds": [
+ 101985.0,
+ 2611485.0,
+ 339315.0,
+ 2826915.0
+ ],
+ "lnglat": [
+ -77.75790625255473,
+ 24.561583285327067
+ ],
+ "height": 718,
+ "width": 791,
+ "shape": [
+ 718,
+ 791
+ ],
+ "res": [
+ 300.0379266750948,
+ 300.041782729805
+ ],
+ "nodata": 0.0
+ }
+
+More information, such as band statistics, can be had using the `--verbose`
+option.
+
+.. code-block:: console
+
+ $ rio info tests/data/RGB.byte.tif --indent 2
+ {
+ "count": 3,
+ "crs": "EPSG:32618",
+ "stats": [
+ {
+ "max": 255.0,
+ "mean": 44.434478650699106,
+ "min": 1.0
+ },
+ {
+ "max": 255.0,
+ "mean": 66.02203484105824,
+ "min": 1.0
+ },
+ {
+ "max": 255.0,
+ "mean": 71.39316199120559,
+ "min": 1.0
+ }
+ ],
+ "dtype": "uint8",
+ "driver": "GTiff",
+ "bounds": [
+ 101985.0,
+ 2611485.0,
+ 339315.0,
+ 2826915.0
+ ],
+ "lnglat": [
+ -77.75790625255473,
+ 24.561583285327067
+ ],
+ "height": 718,
+ "width": 791,
+ "shape": [
+ 718,
+ 791
+ ],
+ "res": [
+ 300.0379266750948,
+ 300.041782729805
+ ],
+ "nodata": 0.0
+ }
insp
----
@@ -113,25 +178,12 @@ The insp command opens a dataset and an interpreter.
.. code-block:: console
$ rio insp tests/data/RGB.byte.tif
- Rasterio 0.9 Interactive Inspector (Python 2.7.5)
+ Rasterio 0.18 Interactive Inspector (Python 2.7.9)
Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
- >>> import pprint
- >>> pprint.pprint(src.meta)
- {'affine': Affine(300.0379266750948, 0.0, 101985.0,
- 0.0, -300.041782729805, 2826915.0),
- 'count': 3,
- 'crs': {'init': u'epsg:32618'},
- 'driver': u'GTiff',
- 'dtype': <type 'numpy.uint8'>,
- 'height': 718,
- 'nodata': 0.0,
- 'transform': (101985.0,
- 300.0379266750948,
- 0.0,
- 2826915.0,
- 0.0,
- -300.041782729805),
- 'width': 791}
+ >>> print src.name
+ tests/data/RGB.byte.tif
+ >>> print src.bounds
+ BoundingBox(left=101985.0, bottom=2611485.0, right=339315.0, top=2826915.0)
merge
-----
@@ -143,6 +195,69 @@ datasets.
$ rio merge rasterio/tests/data/R*.tif merged.tif
+rasterize
+---------
+
+New in 0.18.
+
+The rasterize command rasterizes GeoJSON features into a new or existing
+raster.
+
+.. code-block:: console
+
+ $ rio rasterize test.tif --res 0.0167 < input.geojson
+
+The resulting file will have an upper left coordinate determined by the bounds
+of the GeoJSON (in EPSG:4326, which is the default), with a
+pixel size of approximately 30 arc seconds. Pixels whose center is within the
+polygon or that are selected by brezenhams line algorithm will be burned in
+with a default value of 1.
+
+It is possible to rasterize into an existing raster and use an alternative
+default value:
+
+.. code-block:: console
+
+ $ rio rasterize existing.tif --default_value 10 < input.geojson
+
+It is also possible to rasterize using a template raster, which will be used
+to determine the transform, dimensions, and coordinate reference system of the
+output raster:
+
+.. code-block:: console
+
+ $ rio rasterize test.tif --like tests/data/shade.tif < input.geojson
+
+GeoJSON features may be provided using stdin or specified directly as first
+argument, and dimensions may be provided in place of pixel resolution:
+
+.. code-block:: console
+
+ $ rio rasterize input.geojson test.tif --dimensions 1024 1024
+
+Other options are available, see:
+
+.. code-block:: console
+
+ $ rio rasterize --help
+
+sample
+------
+
+New in 0.18.
+
+The sample command reads ``x, y`` positions from stdin and writes the dataset
+values at that position to stdout.
+
+.. code-block:: console
+
+ $ cat << EOF | rio sample tests/data/RGB.byte.tif
+ > [220649.99999832606, 2719199.999999095]
+ > EOF
+ [18, 25, 14]
+
+The output of the transform command (see below) makes good input for sample.
+
shapes
------
diff --git a/docs/reproject.rst b/docs/reproject.rst
index ad290b4..0b467f4 100644
--- a/docs/reproject.rst
+++ b/docs/reproject.rst
@@ -53,9 +53,57 @@ transform.
assert destination.any()
assert not destination.all()
+
See `examples/reproject.py
<https://github.com/mapbox/rasterio/blob/master/examples/reproject.py>`__ for
code that writes the destination array to a GeoTIFF file. I've
uploaded the resulting file to a Mapbox map to demonstrate that the
reprojection is
-correct:
https://a.tiles.mapbox.com/v3/sgillies.hfek2oko/page.html?secure=1#6/0.000/0.033
+correct:
https://a.tiles.mapbox.com/v3/sgillies.hfek2oko/page.html?secure=1#6/0.000/0.033.
+
+Reprojecting a GeoTIFF dataset
+------------------------------
+
+Here's a more real-world example of using ``reproject()`` to make an output
dataset zoomed out by a factor of 2.
+Methods of the ``rasterio.Affine`` class help us generate the output dataset's
transform matrix and, thereby, its
+spatial extent.
+
+.. code-block:: python
+
+ import numpy
+ import rasterio
+ from rasterio import Affine as A
+ from rasterio.warp import reproject, RESAMPLING
+
+ with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
+ src_transform = src.affine
+
+ # Zoom out by a factor of 2 from the center of the source
+ # dataset. The destination transform is the product of the
+ # source transform, a translation down and to the right, and
+ # a scaling.
+ dst_transform = src_transform*A.translation(
+ -src.width/2.0, -src.height/2.0)*A.scale(2.0)
+
+ data = src.read()
+
+ kwargs = src.meta
+ kwargs['transform'] = dst_transform
+
+ with rasterio.open('/tmp/zoomed-out.tif', 'w', **kwargs) as dst:
+
+ for i, band in enumerate(data, 1):
+ dest = numpy.zeros_like(band)
+
+ reproject(
+ band,
+ dest,
+ src_transform=src_transform,
+ src_crs=src.crs,
+ dst_transform=dst_transform,
+ dst_crs=src.crs,
+ resampling=RESAMPLING.nearest)
+
+ dst.write_band(i, dest)
+
+.. image:: https://farm8.staticflickr.com/7399/16390100651_54f01b8601_b_d.jpg)
References
----------
diff --git a/rasterio/__init__.py b/rasterio/__init__.py
index 531dae2..0ea3dd2 100644
--- a/rasterio/__init__.py
+++ b/rasterio/__init__.py
@@ -18,7 +18,7 @@ from rasterio.transform import Affine, guard_transform
__all__ = [
'band', 'open', 'drivers', 'copy', 'pad']
-__version__ = "0.17"
+__version__ = "0.18"
log = logging.getLogger('rasterio')
class NullHandler(logging.Handler):
diff --git a/rasterio/_base.pxd b/rasterio/_base.pxd
index ad5440b..ca070b8 100644
--- a/rasterio/_base.pxd
+++ b/rasterio/_base.pxd
@@ -23,3 +23,4 @@ cdef class DatasetReader:
cdef void *band(self, int bidx)
+cdef void *_osr_from_crs(object crs)
diff --git a/rasterio/_base.pyx b/rasterio/_base.pyx
index 660d112..96c7540 100644
--- a/rasterio/_base.pyx
+++ b/rasterio/_base.pyx
@@ -99,13 +99,13 @@ cdef class DatasetReader(object):
def read_crs(self):
cdef char *proj_c = NULL
- cdef char *auth_key = NULL
- cdef char *auth_val = NULL
+ cdef const char * auth_key = NULL
+ cdef const char * auth_val = NULL
cdef void *osr = NULL
if self._hds == NULL:
raise ValueError("Null dataset")
crs = {}
- cdef char *wkt = _gdal.GDALGetProjectionRef(self._hds)
+ cdef const char * wkt = _gdal.GDALGetProjectionRef(self._hds)
if wkt is NULL:
raise ValueError("Unexpected NULL spatial reference")
wkt_b = wkt
@@ -166,7 +166,7 @@ cdef class DatasetReader(object):
cdef char *proj_c = NULL
cdef char *key_c = NULL
cdef void *osr = NULL
- cdef char *wkt = NULL
+ cdef const char * wkt = NULL
if self._hds == NULL:
raise ValueError("Null dataset")
wkt = _gdal.GDALGetProjectionRef(self._hds)
@@ -379,6 +379,11 @@ cdef class DatasetReader(object):
lr = self.index(right, bottom)
return tuple(zip(ul, lr))
+ def window_transform(self, window):
+ """Returns the affine transform for a dataset window."""
+ (r, _), (c, _) = window
+ return self.affine * Affine.translation(c or 0, r or 0)
+
@property
def meta(self):
m = {
@@ -394,7 +399,14 @@ cdef class DatasetReader(object):
self._read = True
return m
-
+ def lnglat(self):
+ w, s, e, n = self.bounds
+ cx = (w + e)/2.0
+ cy = (s + n)/2.0
+ lng, lat = _transform(
+ self.crs, {'init': 'epsg:4326'}, [cx], [cy], None)
+ return lng.pop(), lat.pop()
+
def get_crs(self):
# _read tells us that the CRS was read before and really is
# None.
@@ -498,7 +510,7 @@ cdef class DatasetReader(object):
item_c = papszStrList[i]
item_b = item_c
item = item_b.decode('utf-8')
- key, value = item.split('=')
+ key, value = item.split('=', 1)
retval[key] = value
return retval
@@ -521,7 +533,7 @@ cdef class DatasetReader(object):
cdef void *hBand
cdef void *hTable
cdef int i
- cdef _gdal.GDALColorEntry *color
+ cdef const _gdal.GDALColorEntry * color
if self._hds == NULL:
raise ValueError("can't read closed raster file")
if bidx not in self.indexes:
@@ -590,6 +602,7 @@ cpdef eval_window(object window, int height, int width):
"invalid window: col range (%d, %d)" % (c_start, c_stop))
return (r_start, r_stop), (c_start, c_stop)
+
def window_shape(window, height=-1, width=-1):
"""Returns shape of a window.
@@ -599,8 +612,100 @@ def window_shape(window, height=-1, width=-1):
(a, b), (c, d) = eval_window(window, height, width)
return b-a, d-c
+
def window_index(window):
return tuple(slice(*w) for w in window)
+
def tastes_like_gdal(t):
return t[2] == t[4] == 0.0 and t[1] > 0 and t[5] < 0
+
+
+cdef void *_osr_from_crs(object crs):
+ cdef char *proj_c = NULL
+ cdef void *osr = _gdal.OSRNewSpatialReference(NULL)
+ params = []
+ # Normally, we expect a CRS dict.
+ if isinstance(crs, dict):
+ # EPSG is a special case.
+ init = crs.get('init')
+ if init:
+ auth, val = init.split(':')
+ if auth.upper() == 'EPSG':
+ _gdal.OSRImportFromEPSG(osr, int(val))
+ else:
+ crs['wktext'] = True
+ for k, v in crs.items():
+ if v is True or (k in ('no_defs', 'wktext') and v):
+ params.append("+%s" % k)
+ else:
+ params.append("+%s=%s" % (k, v))
+ proj = " ".join(params)
+ log.debug("PROJ.4 to be imported: %r", proj)
+ proj_b = proj.encode('utf-8')
+ proj_c = proj_b
+ _gdal.OSRImportFromProj4(osr, proj_c)
+ # Fall back for CRS strings like "EPSG:3857."
+ else:
+ proj_b = crs.encode('utf-8')
+ proj_c = proj_b
+ _gdal.OSRSetFromUserInput(osr, proj_c)
+ return osr
+
+
+def _transform(src_crs, dst_crs, xs, ys, zs):
+ cdef double *x = NULL
+ cdef double *y = NULL
+ cdef double *z = NULL
+ cdef char *proj_c = NULL
+ cdef void *src = NULL
+ cdef void *dst = NULL
+ cdef void *transform = NULL
+ cdef int i
+
+ assert len(xs) == len(ys)
+ assert zs is None or len(xs) == len(zs)
+
+ src = _osr_from_crs(src_crs)
+ dst = _osr_from_crs(dst_crs)
+
+ n = len(xs)
+ x = <double *>_gdal.CPLMalloc(n*sizeof(double))
+ y = <double *>_gdal.CPLMalloc(n*sizeof(double))
+ for i in range(n):
+ x[i] = xs[i]
+ y[i] = ys[i]
+
+ if zs is not None:
+ z = <double *>_gdal.CPLMalloc(n*sizeof(double))
+ for i in range(n):
+ z[i] = zs[i]
+
+ transform = _gdal.OCTNewCoordinateTransformation(src, dst)
+ res = _gdal.OCTTransform(transform, n, x, y, z)
+ #if res:
+ # raise ValueError("Failed coordinate transformation")
+
+ res_xs = [0]*n
+ res_ys = [0]*n
+
+ for i in range(n):
+ res_xs[i] = x[i]
+ res_ys[i] = y[i]
+
+ if zs is not None:
+ res_zs = [0]*n
+ for i in range(n):
+ res_zs[i] = z[i]
+ _gdal.CPLFree(z)
+
+ retval = (res_xs, res_ys, res_zs)
+ else:
+ retval = (res_xs, res_ys)
+
+ _gdal.CPLFree(x)
+ _gdal.CPLFree(y)
+ _gdal.OCTDestroyCoordinateTransformation(transform)
+ _gdal.OSRDestroySpatialReference(src)
+ _gdal.OSRDestroySpatialReference(dst)
+ return retval
diff --git a/rasterio/_drivers.pyx b/rasterio/_drivers.pyx
index b8080da..40b8271 100644
--- a/rasterio/_drivers.pyx
+++ b/rasterio/_drivers.pyx
@@ -129,8 +129,8 @@ cdef class GDALEnv(object):
def drivers(self):
cdef void *drv = NULL
- cdef char *key = NULL
- cdef char *val = NULL
+ cdef const char *key = NULL
+ cdef const char *val = NULL
cdef int i
result = {}
for i in range(GDALGetDriverCount()):
diff --git a/rasterio/_err.pyx b/rasterio/_err.pyx
index 44bb201..b59ce44 100644
--- a/rasterio/_err.pyx
+++ b/rasterio/_err.pyx
@@ -61,7 +61,7 @@ cdef class GDALErrCtxManager:
def __exit__(self, exc_type=None, exc_val=None, exc_tb=None):
cdef int err_type = CPLGetLastErrorType()
cdef int err_no = CPLGetLastErrorNo()
- cdef char *msg = CPLGetLastErrorMsg()
+ cdef const char *msg = CPLGetLastErrorMsg()
# TODO: warn for err_type 2?
if err_type >= 3:
raise exception_map[err_no](msg)
diff --git a/rasterio/_example.pyx b/rasterio/_example.pyx
index c203847..e1b4396 100644
--- a/rasterio/_example.pyx
+++ b/rasterio/_example.pyx
@@ -1,3 +1,5 @@
+# cython: boundscheck=False
+
import numpy
cimport numpy
diff --git a/rasterio/_features.pyx b/rasterio/_features.pyx
index 3858af4..c7c63c6 100644
--- a/rasterio/_features.pyx
+++ b/rasterio/_features.pyx
@@ -1,4 +1,3 @@
-# cython: profile=True
import logging
import json
@@ -50,12 +49,12 @@ def _shapes(image, mask, connectivity, transform):
"""
cdef int retval, rows, cols
- cdef void *hband
- cdef void *hmaskband
- cdef void *hfdriver
- cdef void *hfs
- cdef void *hlayer
- cdef void *fielddefn
+ cdef void *hband = NULL
+ cdef void *hmaskband = NULL
+ cdef void *hfdriver = NULL
+ cdef void *hfs = NULL
+ cdef void *hlayer = NULL
+ cdef void *fielddefn = NULL
cdef _io.RasterReader rdr
cdef _io.RasterReader mrdr
cdef char **options = NULL
@@ -162,9 +161,9 @@ def _sieve(image, size, output, mask, connectivity):
cdef InMemoryRaster in_mem_ds = None
cdef InMemoryRaster out_mem_ds = None
cdef InMemoryRaster mask_mem_ds = None
- cdef void *in_band
- cdef void *out_band
- cdef void *mask_band
+ cdef void *in_band = NULL
+ cdef void *out_band = NULL
+ cdef void *mask_band = NULL
cdef _io.RasterReader rdr
cdef _io.RasterUpdater udr
cdef _io.RasterReader mask_reader
@@ -291,6 +290,45 @@ def _rasterize(shapes, image, transform, all_touched):
_gdal.CSLDestroy(options)
+def _explode(coords):
+ """Explode a GeoJSON geometry's coordinates object and yield
+ coordinate tuples. As long as the input is conforming, the type of
+ the geometry doesn't matter. From Fiona 1.4.8"""
+ for e in coords:
+ if isinstance(e, (float, int)):
+ yield coords
+ break
+ else:
+ for f in _explode(e):
+ yield f
+
+
+def _bounds(geometry):
+ """Bounding box of a GeoJSON geometry.
+ From Fiona 1.4.8 with updates here to handle feature collections.
+ TODO: add to Fiona.
+ """
+
+ try:
+ if 'features' in geometry:
+ xmins = []
+ ymins = []
+ xmaxs = []
+ ymaxs = []
+ for feature in geometry['features']:
+ xmin, ymin, xmax, ymax = _bounds(feature['geometry'])
+ xmins.append(xmin)
+ ymins.append(ymin)
+ xmaxs.append(xmax)
+ ymaxs.append(ymax)
+ return min(xmins), min(ymins), max(xmaxs), max(ymaxs)
+ else:
+ xyz = tuple(zip(*list(_explode(geometry['coordinates']))))
+ return min(xyz[0]), min(xyz[1]), max(xyz[0]), max(xyz[1])
+ except (KeyError, TypeError):
+ return None
+
+
# Mapping of OGR integer geometry types to GeoJSON type names.
GEOMETRY_TYPES = {
0: 'Unknown',
@@ -416,12 +454,6 @@ cdef class GeomBuilder:
return result
-cdef geometry(void *geom):
- """Returns a GeoJSON object from an OGR geometry object."""
-
- return GeomBuilder().build(geom)
-
-
cdef class OGRGeomBuilder:
"""
Builds an OGR geometry from GeoJSON geometry.
diff --git a/rasterio/_fill.pyx b/rasterio/_fill.pyx
new file mode 100644
index 0000000..8d7c998
--- /dev/null
+++ b/rasterio/_fill.pyx
@@ -0,0 +1,84 @@
+# distutils: language = c++
+# cython: profile=True
+#
+
+import numpy as np
+cimport numpy as np
+
+from rasterio import dtypes
+from rasterio._err import cpl_errs
+from rasterio cimport _gdal, _io
+
+from rasterio._io cimport InMemoryRaster
+
+
+def _fillnodata(image, mask, double max_search_distance=100.0,
+ int smoothing_iterations=0):
+ cdef void *memdriver = _gdal.GDALGetDriverByName("MEM")
+ cdef void *image_dataset
+ cdef void *image_band
+ cdef void *mask_dataset
+ cdef void *mask_band
+ cdef char **alg_options = NULL
+
+ if isinstance(image, np.ndarray):
+ # copy numpy ndarray into an in-memory dataset
+ image_dataset = _gdal.GDALCreate(
+ memdriver,
+ "image",
+ image.shape[1],
+ image.shape[0],
+ 1,
+ <_gdal.GDALDataType>dtypes.dtype_rev[image.dtype.name],
+ NULL)
+ image_band = _gdal.GDALGetRasterBand(image_dataset, 1)
+ _io.io_auto(image, image_band, True)
+ elif isinstance(image, tuple):
+ # TODO
+ raise NotImplementedError()
+ else:
+ raise ValueError("Invalid source image")
+
+ if isinstance(mask, np.ndarray):
+ mask_cast = mask.astype('uint8')
+ mask_dataset = _gdal.GDALCreate(
+ memdriver,
+ "mask",
+ mask.shape[1],
+ mask.shape[0],
+ 1,
+ <_gdal.GDALDataType>dtypes.dtype_rev['uint8'],
+ NULL)
+ mask_band = _gdal.GDALGetRasterBand(mask_dataset, 1)
+ _io.io_auto(mask_cast, mask_band, True)
+ elif isinstance(mask, tuple):
+ # TODO
+ raise NotImplementedError()
+ elif mask is None:
+ mask_band = NULL
+ else:
+ raise ValueError("Invalid source image mask")
+
+ with cpl_errs:
+ alg_options = _gdal.CSLSetNameValue(
+ alg_options, "TEMP_FILE_DRIVER", "MEM")
+
+ _gdal.GDALFillNodata(
+ image_band,
+ mask_band,
+ max_search_distance,
+ 0,
+ smoothing_iterations,
+ alg_options,
+ NULL,
+ NULL)
+
+ # read the result into a numpy ndarray
+ result = np.empty(image.shape, dtype=image.dtype)
+ _io.io_auto(result, image_band, False)
+
+ _gdal.GDALClose(image_dataset)
+ _gdal.GDALClose(mask_dataset)
+ _gdal.CSLDestroy(alg_options)
+
+ return result
diff --git a/rasterio/_gdal.pxd b/rasterio/_gdal.pxd
index 1214678..0ae46f6 100644
--- a/rasterio/_gdal.pxd
+++ b/rasterio/_gdal.pxd
@@ -101,7 +101,7 @@ cdef extern from "gdal.h" nogil:
short c3
short c4
- const GDALColorEntry *GDALGetColorEntry (void *hTable, int)
+ const GDALColorEntry * GDALGetColorEntry (void *hTable, int)
void GDALSetColorEntry (void *hTable, int i, const GDALColorEntry *poEntry)
int GDALSetRasterColorTable (void *hBand, void *hTable)
void *GDALGetRasterColorTable (void *hBand)
@@ -211,5 +211,5 @@ cdef extern from "gdal_alg.h":
int GDALApproxTransform(void *pTransformArg, int bDstToSrc, int
nPointCount, double *x, double *y, double *z, int *panSuccess )
void GDALDestroyApproxTransformer( void * )
-
+ 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)
diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx
index 98e18e1..e9339a5 100644
--- a/rasterio/_io.pyx
+++ b/rasterio/_io.pyx
@@ -559,8 +559,13 @@ cdef class RasterReader(_base.DatasetReader):
a 2D array if it is a band index number.
out: numpy ndarray, optional
- An optional reference to an output array with the same
- dimensions and shape.
+ As with Numpy ufuncs, this is an optional reference to an
+ output array with the same dimensions and shape into which
+ data will be placed.
+
+ *Note*: the method's return value may be a view on this
+ array. In other words, `out` is likely to be an
+ incomplete representation of the method's results.
window : a pair (tuple) of pairs of ints, optional
The optional `window` argument is a 2 item tuple. The first
@@ -585,8 +590,11 @@ cdef class RasterReader(_base.DatasetReader):
Returns
-------
- Numpy ndarray
+ Numpy ndarray or a view on a Numpy ndarray
+ Note: as with Numpy ufuncs, an object is returned even if you
+ use the optional `out` argument and the return value shall be
+ preferentially used by callers.
"""
return2d = False
@@ -855,6 +863,29 @@ cdef class RasterReader(_base.DatasetReader):
hmask, 0, xoff, yoff, width, height, out)
return out
+ def sample(self, xy, indexes=None):
+ """Get the values of a dataset at certain positions
+
+ Parameters
+ ----------
+ xy : iterable, pairs of floats
+ A sequence or generator of (x, y) pairs.
+
+ indexes : list of ints or a single int, optional
+ If `indexes` is a list, the result is a 3D array, but is
+ a 2D array if it is a band index number.
+
+ Returns
+ -------
+ Iterable, yielding dataset values for the specified `indexes`
+ as an ndarray.
+ """
+ for x, y in xy:
+ r, c = self.index(x, y)
+ window = ((r, r+1), (c, c+1))
+ data = self.read(
+ indexes, window=window, masked=False, boundless=True)
+ yield data[:,0,0]
cdef class RasterUpdater(RasterReader):
# Read-write access to raster data and metadata.
@@ -911,7 +942,8 @@ cdef class RasterUpdater(RasterReader):
def start(self):
cdef const char *drv_name = NULL
cdef char **options = NULL
- cdef char *key_c, *val_c = NULL
+ cdef char *key_c = NULL
+ cdef char *val_c = NULL
cdef void *drv = NULL
cdef void *hband = NULL
cdef int success
@@ -1308,7 +1340,7 @@ cdef class RasterUpdater(RasterReader):
def write_colormap(self, bidx, colormap):
"""Write a colormap for a band to the dataset."""
- cdef void *hBand
+ cdef void *hBand = NULL
cdef void *hTable
cdef _gdal.GDALColorEntry color
if self._hds == NULL:
@@ -1412,6 +1444,7 @@ cdef class InMemoryRaster:
self.dataset = NULL
cdef void *memdriver = _gdal.GDALGetDriverByName("MEM")
+ cdef int i = 0 # avoids Cython warning in for loop below
# Several GDAL operations require the array of band IDs as input
self.band_ids[0] = 1
@@ -1548,7 +1581,8 @@ cdef class IndirectRasterUpdater(RasterUpdater):
def close(self):
cdef const char *drv_name = NULL
cdef char **options = NULL
- cdef char *key_c, *val_c = NULL
+ cdef char *key_c = NULL
+ cdef char *val_c = NULL
cdef void *drv = NULL
cdef void *temp = NULL
cdef int success
@@ -1598,7 +1632,7 @@ def writer(path, mode, **kwargs):
# format driver's capabilities.
cdef void *hds = NULL
cdef void *drv = NULL
- cdef char *drv_name = NULL
+ cdef const char *drv_name = NULL
cdef const char *fname = NULL
if mode == 'w' and 'driver' in kwargs:
diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx
index d9faaed..b84470a 100644
--- a/rasterio/_warp.pyx
+++ b/rasterio/_warp.pyx
@@ -1,6 +1,4 @@
# distutils: language = c++
-# cython: profile=True
-#
from collections import namedtuple
import logging
@@ -8,7 +6,7 @@ import logging
import numpy as np
cimport numpy as np
-from rasterio cimport _gdal, _ogr, _io, _features
+from rasterio cimport _base, _gdal, _ogr, _io, _features
from rasterio import dtypes
@@ -78,94 +76,6 @@ def tastes_like_gdal(t):
return t[2] == t[4] == 0.0 and t[1] > 0 and t[5] < 0
-cdef void *_osr_from_crs(object crs):
- cdef char *proj_c = NULL
- cdef void *osr
- osr = _gdal.OSRNewSpatialReference(NULL)
- params = []
- # Normally, we expect a CRS dict.
- if isinstance(crs, dict):
- # EPSG is a special case.
- init = crs.get('init')
- if init:
- auth, val = init.split(':')
- if auth.upper() == 'EPSG':
- _gdal.OSRImportFromEPSG(osr, int(val))
- else:
- crs['wktext'] = True
- for k, v in crs.items():
- if v is True or (k in ('no_defs', 'wktext') and v):
- params.append("+%s" % k)
- else:
- params.append("+%s=%s" % (k, v))
- proj = " ".join(params)
- log.debug("PROJ.4 to be imported: %r", proj)
- proj_b = proj.encode('utf-8')
- proj_c = proj_b
- _gdal.OSRImportFromProj4(osr, proj_c)
- # Fall back for CRS strings like "EPSG:3857."
- else:
- proj_b = crs.encode('utf-8')
- proj_c = proj_b
- _gdal.OSRSetFromUserInput(osr, proj_c)
- return osr
-
-
-def _transform(src_crs, dst_crs, xs, ys, zs):
- cdef double *x, *y, *z = NULL
- cdef char *proj_c = NULL
- cdef void *src, *dst
- cdef void *transform
- cdef int i
-
- assert len(xs) == len(ys)
- assert zs is None or len(xs) == len(zs)
-
- src = _osr_from_crs(src_crs)
- dst = _osr_from_crs(dst_crs)
-
- n = len(xs)
- x = <double *>_gdal.CPLMalloc(n*sizeof(double))
- y = <double *>_gdal.CPLMalloc(n*sizeof(double))
- for i in range(n):
- x[i] = xs[i]
- y[i] = ys[i]
-
- if zs is not None:
- z = <double *>_gdal.CPLMalloc(n*sizeof(double))
- for i in range(n):
- z[i] = zs[i]
-
- transform = _gdal.OCTNewCoordinateTransformation(src, dst)
- res = _gdal.OCTTransform(transform, n, x, y, z)
- #if res:
- # raise ValueError("Failed coordinate transformation")
-
- res_xs = [0]*n
- res_ys = [0]*n
-
- for i in range(n):
- res_xs[i] = x[i]
- res_ys[i] = y[i]
-
- if zs is not None:
- res_zs = [0]*n
- for i in range(n):
- res_zs[i] = z[i]
- _gdal.CPLFree(z)
-
- retval = (res_xs, res_ys, res_zs)
- else:
- retval = (res_xs, res_ys)
-
- _gdal.CPLFree(x)
- _gdal.CPLFree(y)
- _gdal.OCTDestroyCoordinateTransformation(transform)
- _gdal.OSRDestroySpatialReference(src)
- _gdal.OSRDestroySpatialReference(dst)
- return retval
-
-
def _transform_geom(
src_crs, dst_crs, geom, antimeridian_cutting, antimeridian_offset,
precision):
@@ -174,15 +84,16 @@ def _transform_geom(
cdef char *key_c = NULL
cdef char *val_c = NULL
cdef char **options = NULL
- cdef void *src, *dst
- cdef void *transform
- cdef OGRGeometryFactory *factory
- cdef void *src_ogr_geom
- cdef void *dst_ogr_geom
+ cdef void *src = NULL
+ cdef void *dst = NULL
+ cdef void *transform = NULL
+ cdef OGRGeometryFactory *factory = NULL
+ cdef void *src_ogr_geom = NULL
+ cdef void *dst_ogr_geom = NULL
cdef int i
- src = _osr_from_crs(src_crs)
- dst = _osr_from_crs(dst_crs)
+ src = _base._osr_from_crs(src_crs)
+ dst = _base._osr_from_crs(dst_crs)
transform = _gdal.OCTNewCoordinateTransformation(src, dst)
# Transform options.
@@ -259,7 +170,7 @@ def _reproject(
bands of datasets on disk, the coordinate reference systems and
transforms will be read from the appropriate datasets.
"""
- cdef int retval, rows, cols
+ cdef int retval=0, rows, cols
cdef void *hrdriver
cdef void *hdsin
cdef void *hdsout
diff --git a/rasterio/features.py b/rasterio/features.py
index ea5884a..3673c3a 100644
--- a/rasterio/features.py
+++ b/rasterio/features.py
@@ -8,7 +8,7 @@ import warnings
import numpy as np
import rasterio
-from rasterio._features import _shapes, _sieve, _rasterize
+from rasterio._features import _shapes, _sieve, _rasterize, _bounds
from rasterio.transform import IDENTITY, guard_transform
from rasterio.dtypes import get_minimum_int_dtype
@@ -321,3 +321,24 @@ def rasterize(
_rasterize(valid_shapes, out, transform.to_gdal(), all_touched)
return out
+
+
+def bounds(geometry):
+ """Returns a (minx, miny, maxx, maxy) bounding box. From Fiona 1.4.8.
+ Modified to return bbox from geometry if available.
+
+ Parameters
+ ----------
+ geometry: GeoJSON-like feature, feature collection, or geometry.
+
+ Returns
+ -------
+ tuple
+ Bounding box: (minx, miny, maxx, maxy)
+ """
+
+ if 'bbox' in geometry:
+ return tuple(geometry['bbox'])
+
+ geom = geometry.get('geometry') or geometry
+ return _bounds(geom)
diff --git a/rasterio/fill.py b/rasterio/fill.py
new file mode 100644
index 0000000..37ec58e
--- /dev/null
+++ b/rasterio/fill.py
@@ -0,0 +1,47 @@
+import rasterio
+from rasterio._fill import _fillnodata
+
+def fillnodata(
+ image,
+ mask=None,
+ max_search_distance=100.0,
+ smoothing_iterations=0):
+ """
+ Fill selected raster regions by interpolation from the edges.
+
+ This algorithm will interpolate values for all designated nodata pixels
+ (marked by zeros in `mask`). For each pixel a four direction conic search
+ is done to find values to interpolate from (using inverse distance
+ weighting). Once all values are interpolated, zero or more smoothing
+ iterations (3x3 average filters on interpolated pixels) are applied to
+ smooth out artifacts.
+
+ This algorithm is generally suitable for interpolating missing regions of
+ fairly continuously varying rasters (such as elevation models for
+ instance). It is also suitable for filling small holes and cracks in more
+ irregularly varying images (like aerial photos). It is generally not so
+ great for interpolating a raster from sparse point data.
+
+ Parameters
+ ----------
+ image : numpy ndarray
+ The source containing nodata holes.
+ mask : numpy ndarray
+ A mask band indicating which pixels to interpolate. Pixels to
+ interpolate into are indicated by the value 0. Values of 1 indicate
+ areas to use during interpolation. Must be same shape as image.
+ max_search_distance : float, optional
+ The maxmimum number of pixels to search in all directions to find
+ values to interpolate from. The default is 100.
+ smoothing_iterations : integer, optional
+ The number of 3x3 smoothing filter passes to run. The default is 0.
+
+ Returns
+ -------
+ out : numpy ndarray
+ The filled raster array.
+ """
+ max_search_distance = float(max_search_distance)
+ smoothing_iterations = int(smoothing_iterations)
+ with rasterio.drivers():
+ return _fillnodata(image, mask, max_search_distance,
smoothing_iterations)
diff --git a/rasterio/rasterfill.cpp b/rasterio/rasterfill.cpp
new file mode 100644
index 0000000..df2297b
--- /dev/null
+++ b/rasterio/rasterfill.cpp
@@ -0,0 +1,886 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project: GDAL
+ * Purpose: Interpolate in nodata areas.
+ * Author: Frank Warmerdam, warmer...@pobox.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2008, Frank Warmerdam
+ * Copyright (c) 2015, Sean Gillies <s...@mapbox.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ***************************************************************************/
+
+#include "gdal_alg.h"
+#include "cpl_conv.h"
+#include "cpl_string.h"
+
+CPL_CVSID("$Id$");
+
+/************************************************************************/
+/* GDALFilterLine() */
+/* */
+/* Apply 3x3 filtering one one scanline with masking for which */
+/* pixels are to be interpolated (ThisFMask) and which window */
+/* pixels are valid to include in the interpolation (TMask). */
+/************************************************************************/
+
+static void
+GDALFilterLine( float *pafLastLine, float *pafThisLine, float *pafNextLine,
+ float *pafOutLine,
+ GByte *pabyLastTMask, GByte *pabyThisTMask,
GByte*pabyNextTMask,
+ GByte *pabyThisFMask, int nXSize )
+
+{
+ int iX;
+
+ for( iX = 0; iX < nXSize; iX++ )
+ {
+ if( !pabyThisFMask[iX] )
+ {
+ pafOutLine[iX] = pafThisLine[iX];
+ continue;
+ }
+
+ CPLAssert( pabyThisTMask[iX] );
+
+ double dfValSum = 0.0;
+ double dfWeightSum = 0.0;
+
+ // Previous line
+ if( pafLastLine != NULL )
+ {
+ if( iX > 0 && pabyLastTMask[iX-1] )
+ {
+ dfValSum += pafLastLine[iX-1];
+ dfWeightSum += 1.0;
+ }
+ if( pabyLastTMask[iX] )
+ {
+ dfValSum += pafLastLine[iX];
+ dfWeightSum += 1.0;
+ }
+ if( iX < nXSize-1 && pabyLastTMask[iX+1] )
+ {
+ dfValSum += pafLastLine[iX+1];
+ dfWeightSum += 1.0;
+ }
+ }
+
+ // Current Line
+ if( iX > 0 && pabyThisTMask[iX-1] )
+ {
+ dfValSum += pafThisLine[iX-1];
+ dfWeightSum += 1.0;
+ }
+ if( pabyThisTMask[iX] )
+ {
+ dfValSum += pafThisLine[iX];
+ dfWeightSum += 1.0;
+ }
+ if( iX < nXSize-1 && pabyThisTMask[iX+1] )
+ {
+ dfValSum += pafThisLine[iX+1];
+ dfWeightSum += 1.0;
+ }
+
+ // Next line
+ if( pafNextLine != NULL )
+ {
+ if( iX > 0 && pabyNextTMask[iX-1] )
+ {
+ dfValSum += pafNextLine[iX-1];
+ dfWeightSum += 1.0;
+ }
+ if( pabyNextTMask[iX] )
+ {
+ dfValSum += pafNextLine[iX];
+ dfWeightSum += 1.0;
+ }
+ if( iX < nXSize-1 && pabyNextTMask[iX+1] )
+ {
+ dfValSum += pafNextLine[iX+1];
+ dfWeightSum += 1.0;
+ }
+ }
+
+ pafOutLine[iX] = (float) (dfValSum / dfWeightSum);
+ }
+}
+
+/************************************************************************/
+/* GDALMultiFilter() */
+/* */
+/* Apply multiple iterations of a 3x3 smoothing filter over a */
+/* band with masking controlling what pixels should be */
+/* filtered (FiltMaskBand non zero) and which pixels can be */
+/* considered valid contributors to the filter */
+/* (TargetMaskBand non zero). */
+/* */
+/* This implementation attempts to apply many iterations in */
+/* one IO pass by managing the filtering over a rolling buffer */
+/* of nIternations+2 scanlines. While possibly clever this */
+/* makes the algorithm implementation largely */
+/* incomprehensible. */
+/************************************************************************/
+
+static CPLErr
+GDALMultiFilter( GDALRasterBandH hTargetBand,
+ GDALRasterBandH hTargetMaskBand,
+ GDALRasterBandH hFiltMaskBand,
+ int nIterations,
+ GDALProgressFunc pfnProgress,
+ void * pProgressArg )
+
+{
+ float *paf3PassLineBuf;
+ GByte *pabyTMaskBuf;
+ GByte *pabyFMaskBuf;
+ float *pafThisPass, *pafLastPass, *pafSLastPass;
+
+ int nBufLines = nIterations + 2;
+ int iPassCounter = 0;
+ int nNewLine; // the line being loaded this time (zero based scanline)
+ int nXSize = GDALGetRasterBandXSize( hTargetBand );
+ int nYSize = GDALGetRasterBandYSize( hTargetBand );
+ CPLErr eErr = CE_None;
+
+/* -------------------------------------------------------------------- */
+/* Report starting progress value. */
+/* -------------------------------------------------------------------- */
+ if( !pfnProgress( 0.0, "Smoothing Filter...", pProgressArg ) )
+ {
+ CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+ return CE_Failure;
+ }
+
+/* -------------------------------------------------------------------- */
+/* Allocate rotating buffers. */
+/* -------------------------------------------------------------------- */
+ pabyTMaskBuf = (GByte *) VSIMalloc2(nXSize, nBufLines);
+ pabyFMaskBuf = (GByte *) VSIMalloc2(nXSize, nBufLines);
+
+ paf3PassLineBuf = (float *) VSIMalloc3(nXSize, nBufLines, 3 *
sizeof(float));
+ if (pabyTMaskBuf == NULL || pabyFMaskBuf == NULL || paf3PassLineBuf ==
NULL)
+ {
+ CPLError(CE_Failure, CPLE_OutOfMemory,
+ "Could not allocate enough memory for temporary buffers");
+ eErr = CE_Failure;
+ goto end;
+ }
+
+/* -------------------------------------------------------------------- */
+/* Process rotating buffers. */
+/* -------------------------------------------------------------------- */
+ for( nNewLine = 0;
+ eErr == CE_None && nNewLine < nYSize+nIterations;
+ nNewLine++ )
+ {
+/* -------------------------------------------------------------------- */
+/* Rotate pass buffers. */
+/* -------------------------------------------------------------------- */
+ iPassCounter = (iPassCounter + 1) % 3;
+
+ pafSLastPass = paf3PassLineBuf
+ + ((iPassCounter+0)%3) * nXSize*nBufLines;
+ pafLastPass = paf3PassLineBuf
+ + ((iPassCounter+1)%3) * nXSize*nBufLines;
+ pafThisPass = paf3PassLineBuf
+ + ((iPassCounter+2)%3) * nXSize*nBufLines;
+
+/* -------------------------------------------------------------------- */
+/* Where does the new line go in the rotating buffer? */
+/* -------------------------------------------------------------------- */
+ int iBufOffset = nNewLine % nBufLines;
+
+/* -------------------------------------------------------------------- */
+/* Read the new data line if it is't off the bottom of the */
+/* image. */
+/* -------------------------------------------------------------------- */
+ if( nNewLine < nYSize )
+ {
+ eErr =
+ GDALRasterIO( hTargetMaskBand, GF_Read,
+ 0, nNewLine, nXSize, 1,
+ pabyTMaskBuf + nXSize * iBufOffset, nXSize, 1,
+ GDT_Byte, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+ eErr =
+ GDALRasterIO( hFiltMaskBand, GF_Read,
+ 0, nNewLine, nXSize, 1,
+ pabyFMaskBuf + nXSize * iBufOffset, nXSize, 1,
+ GDT_Byte, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+ eErr =
+ GDALRasterIO( hTargetBand, GF_Read,
+ 0, nNewLine, nXSize, 1,
+ pafThisPass + nXSize * iBufOffset, nXSize, 1,
+ GDT_Float32, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+ }
+
+/* -------------------------------------------------------------------- */
+/* Loop over the loaded data, applying the filter to all loaded */
+/* lines with neighbours. */
+/* -------------------------------------------------------------------- */
+ int iFLine;
+
+ for( iFLine = nNewLine-1;
+ eErr == CE_None && iFLine >= nNewLine-nIterations;
+ iFLine-- )
+ {
+ int iLastOffset, iThisOffset, iNextOffset;
+
+ iLastOffset = (iFLine-1) % nBufLines;
+ iThisOffset = (iFLine ) % nBufLines;
+ iNextOffset = (iFLine+1) % nBufLines;
+
+ // default to preserving the old value.
+ if( iFLine >= 0 )
+ memcpy( pafThisPass + iThisOffset * nXSize,
+ pafLastPass + iThisOffset * nXSize,
+ sizeof(float) * nXSize );
+
+ // currently this skips the first and last line. Eventually
+ // we will enable these too. TODO
+ if( iFLine < 1 || iFLine >= nYSize-1 )
+ {
+ continue;
+ }
+
+ GDALFilterLine(
+ pafSLastPass + iLastOffset * nXSize,
+ pafLastPass + iThisOffset * nXSize,
+ pafThisPass + iNextOffset * nXSize,
+ pafThisPass + iThisOffset * nXSize,
+ pabyTMaskBuf + iLastOffset * nXSize,
+ pabyTMaskBuf + iThisOffset * nXSize,
+ pabyTMaskBuf + iNextOffset * nXSize,
+ pabyFMaskBuf + iThisOffset * nXSize,
+ nXSize );
+ }
+
+/* -------------------------------------------------------------------- */
+/* Write out the top data line that will be rolling out of our */
+/* buffer. */
+/* -------------------------------------------------------------------- */
+ int iLineToSave = nNewLine - nIterations;
+
+ if( iLineToSave >= 0 && eErr == CE_None )
+ {
+ iBufOffset = iLineToSave % nBufLines;
+
+ eErr =
+ GDALRasterIO( hTargetBand, GF_Write,
+ 0, iLineToSave, nXSize, 1,
+ pafThisPass + nXSize * iBufOffset, nXSize, 1,
+ GDT_Float32, 0, 0 );
+ }
+
+/* -------------------------------------------------------------------- */
+/* Report progress. */
+/* -------------------------------------------------------------------- */
+ if( eErr == CE_None
+ && !pfnProgress( (nNewLine+1) / (double) (nYSize+nIterations),
+ "Smoothing Filter...", pProgressArg ) )
+ {
+ CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+ eErr = CE_Failure;
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Cleanup */
+/* -------------------------------------------------------------------- */
+end:
+ CPLFree( pabyTMaskBuf );
+ CPLFree( pabyFMaskBuf );
+ CPLFree( paf3PassLineBuf );
+
+ return eErr;
+}
+
+/************************************************************************/
+/* QUAD_CHECK() */
+/* */
+/* macro for checking whether a point is nearer than the */
+/* existing closest point. */
+/************************************************************************/
+#define QUAD_CHECK(quad_dist, quad_value, \
+target_x, target_y, origin_x, origin_y, target_value ) \
+ \
+if( quad_value != nNoDataVal )
\
+{ \
+ double dfDx = (double)target_x - (double)origin_x; \
+ double dfDy = (double)target_y - (double)origin_y; \
+ double dfDistSq = dfDx * dfDx + dfDy * dfDy; \
+ \
+ if( dfDistSq < quad_dist*quad_dist ) \
+ { \
+ CPLAssert( dfDistSq > 0.0 ); \
+ quad_dist = sqrt(dfDistSq); \
+ quad_value = target_value; \
+ } \
+}
+
+/************************************************************************/
+/* GDALFillNodata() */
+/************************************************************************/
+
+/**
+ * Fill selected raster regions by interpolation from the edges.
+ *
+ * This algorithm will interpolate values for all designated
+ * nodata pixels (marked by zeros in hMaskBand). For each pixel
+ * a four direction conic search is done to find values to interpolate
+ * from (using inverse distance weighting). Once all values are
+ * interpolated, zero or more smoothing iterations (3x3 average
+ * filters on interpolated pixels) are applied to smooth out
+ * artifacts.
+ *
+ * This algorithm is generally suitable for interpolating missing
+ * regions of fairly continuously varying rasters (such as elevation
+ * models for instance). It is also suitable for filling small holes
+ * and cracks in more irregularly varying images (like airphotos). It
+ * is generally not so great for interpolating a raster from sparse
+ * point data - see the algorithms defined in gdal_grid.h for that case.
+ *
+ * @param hTargetBand the raster band to be modified in place.
+ * @param hMaskBand a mask band indicating pixels to be interpolated (zero
valued
+ * @param dfMaxSearchDist the maximum number of pixels to search in all
+ * directions to find values to interpolate from.
+ * @param bDeprecatedOption unused argument, should be zero.
+ * @param nSmoothingIterations the number of 3x3 smoothing filter passes to
+ * run (0 or more).
+ * @param papszOptions additional name=value options in a string list (the
+ * temporary file driver can be specified like TEMP_FILE_DRIVER=MEM).
+ * @param pfnProgress the progress function to report completion.
+ * @param pProgressArg callback data for progress function.
+ *
+ * @return CE_None on success or CE_Failure if something goes wrong.
+ */
+
+CPLErr CPL_STDCALL
+GDALFillNodata( GDALRasterBandH hTargetBand,
+ GDALRasterBandH hMaskBand,
+ double dfMaxSearchDist,
+ int bDeprecatedOption,
+ int nSmoothingIterations,
+ char **papszOptions,
+ GDALProgressFunc pfnProgress,
+ void * pProgressArg )
+
+{
+ VALIDATE_POINTER1( hTargetBand, "GDALFillNodata", CE_Failure );
+
+ int nXSize = GDALGetRasterBandXSize( hTargetBand );
+ int nYSize = GDALGetRasterBandYSize( hTargetBand );
+ CPLErr eErr = CE_None;
+
+ // Special "x" pixel values identifying pixels as special.
+ GUInt32 nNoDataVal;
+ GDALDataType eType;
+
+ if( dfMaxSearchDist == 0.0 )
+ dfMaxSearchDist = MAX(nXSize,nYSize) + 1;
+
+ int nMaxSearchDist = (int) floor(dfMaxSearchDist);
+
+ if( nXSize > 65533 || nYSize > 65533 )
+ {
+ eType = GDT_UInt32;
+ nNoDataVal = 4000002;
+ }
+ else
+ {
+ eType = GDT_UInt16;
+ nNoDataVal = 65535;
+ }
+
+ if( hMaskBand == NULL )
+ hMaskBand = GDALGetMaskBand( hTargetBand );
+
+ /* If there are smoothing iterations, reserve 10% of the progress for them
*/
+ double dfProgressRatio = (nSmoothingIterations > 0) ? 0.9 : 1.0;
+
+/* -------------------------------------------------------------------- */
+/* Initialize progress counter. */
+/* -------------------------------------------------------------------- */
+ if( pfnProgress == NULL )
+ pfnProgress = GDALDummyProgress;
+
+ if( !pfnProgress( 0.0, "Filling...", pProgressArg ) )
+ {
+ CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+ return CE_Failure;
+ }
+
+/* -------------------------------------------------------------------- */
+/* Determine format driver for temp work files. */
+/* -------------------------------------------------------------------- */
+ CPLString osTmpFileDriver = CSLFetchNameValueDef(
+ papszOptions, "TEMP_FILE_DRIVER", "MEM");
+ GDALDriverH hDriver = GDALGetDriverByName((const char *) osTmpFileDriver);
+
+ if (hDriver == NULL)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "Given driver is not registered");
+ return CE_Failure;
+ }
+
+ if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, NULL) == NULL)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "Given driver is incapable of creating temp work files");
+ return CE_Failure;
+ }
+
+ char **papszWorkFileOptions = NULL;
+ if (osTmpFileDriver == "GTiff") {
+ papszWorkFileOptions = CSLSetNameValue(
+ papszWorkFileOptions, "COMPRESS", "LZW");
+ papszWorkFileOptions = CSLSetNameValue(
+ papszWorkFileOptions, "BIGTIFF", "IF_SAFER");
+ }
+
+/* -------------------------------------------------------------------- */
+/* Create a work file to hold the Y "last value" indices. */
+/* -------------------------------------------------------------------- */
+ GDALDatasetH hYDS;
+ GDALRasterBandH hYBand;
+
+ CPLString osTmpFile = CPLGenerateTempFilename("");
+ CPLString osYTmpFile = osTmpFile + "fill_y_work.tif";
+
+ hYDS = GDALCreate( hDriver, osYTmpFile, nXSize, nYSize, 1,
+ eType, (char **) papszWorkFileOptions );
+
+ if ( hYDS == NULL )
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "Could not create Y index work file. Check driver capabilities.");
+ return CE_Failure;
+ }
+
+ hYBand = GDALGetRasterBand( hYDS, 1 );
+
+/* -------------------------------------------------------------------- */
+/* Create a work file to hold the pixel value associated with */
+/* the "last xy value" pixel. */
+/* -------------------------------------------------------------------- */
+ GDALDatasetH hValDS;
+ GDALRasterBandH hValBand;
+ CPLString osValTmpFile = osTmpFile + "fill_val_work.tif";
+
+ hValDS = GDALCreate( hDriver, osValTmpFile, nXSize, nYSize, 1,
+ GDALGetRasterDataType( hTargetBand ),
+ (char **) papszWorkFileOptions );
+
+ if ( hValDS == NULL )
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "Could not create XY value work file. Check driver capabilities.");
+ return CE_Failure;
+ }
+
+ hValBand = GDALGetRasterBand( hValDS, 1 );
+
+/* -------------------------------------------------------------------- */
+/* Create a mask file to make it clear what pixels can be filtered */
+/* on the filtering pass. */
+/* -------------------------------------------------------------------- */
+ GDALDatasetH hFiltMaskDS;
+ GDALRasterBandH hFiltMaskBand;
+ CPLString osFiltMaskTmpFile = osTmpFile + "fill_filtmask_work.tif";
+
+ hFiltMaskDS =
+ GDALCreate( hDriver, osFiltMaskTmpFile, nXSize, nYSize, 1,
+ GDT_Byte, (char **) papszWorkFileOptions );
+
+ if ( hFiltMaskDS == NULL )
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "Could not create mask work file. Check driver capabilities.");
+ return CE_Failure;
+ }
+
+ hFiltMaskBand = GDALGetRasterBand( hFiltMaskDS, 1 );
+
+/* -------------------------------------------------------------------- */
+/* Allocate buffers for last scanline and this scanline. */
+/* -------------------------------------------------------------------- */
+ GUInt32 *panLastY, *panThisY, *panTopDownY;
+ float *pafLastValue, *pafThisValue, *pafScanline, *pafTopDownValue;
+ GByte *pabyMask, *pabyFiltMask;
+ int iX;
+ int iY;
+
+ panLastY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
+ panThisY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
+ panTopDownY = (GUInt32 *) VSICalloc(nXSize,sizeof(GUInt32));
+ pafLastValue = (float *) VSICalloc(nXSize,sizeof(float));
+ pafThisValue = (float *) VSICalloc(nXSize,sizeof(float));
+ pafTopDownValue = (float *) VSICalloc(nXSize,sizeof(float));
+ pafScanline = (float *) VSICalloc(nXSize,sizeof(float));
+ pabyMask = (GByte *) VSICalloc(nXSize,1);
+ pabyFiltMask = (GByte *) VSICalloc(nXSize,1);
+ if (panLastY == NULL || panThisY == NULL || panTopDownY == NULL ||
+ pafLastValue == NULL || pafThisValue == NULL || pafTopDownValue ==
NULL ||
+ pafScanline == NULL || pabyMask == NULL || pabyFiltMask == NULL)
+ {
+ CPLError(CE_Failure, CPLE_OutOfMemory,
+ "Could not allocate enough memory for temporary buffers");
+
+ eErr = CE_Failure;
+ goto end;
+ }
+
+ for( iX = 0; iX < nXSize; iX++ )
+ {
+ panLastY[iX] = nNoDataVal;
+ }
+
+/* ==================================================================== */
+/* Make first pass from top to bottom collecting the "last */
+/* known value" for each column and writing it out to the work */
+/* files. */
+/* ==================================================================== */
+
+ for( iY = 0; iY < nYSize && eErr == CE_None; iY++ )
+ {
+/* -------------------------------------------------------------------- */
+/* Read data and mask for this line. */
+/* -------------------------------------------------------------------- */
+ eErr =
+ GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1,
+ pabyMask, nXSize, 1, GDT_Byte, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+ eErr =
+ GDALRasterIO( hTargetBand, GF_Read, 0, iY, nXSize, 1,
+ pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+/* -------------------------------------------------------------------- */
+/* Figure out the most recent pixel for each column. */
+/* -------------------------------------------------------------------- */
+
+ for( iX = 0; iX < nXSize; iX++ )
+ {
+ if( pabyMask[iX] )
+ {
+ pafThisValue[iX] = pafScanline[iX];
+ panThisY[iX] = iY;
+ }
+ else if( iY <= dfMaxSearchDist + panLastY[iX] )
+ {
+ pafThisValue[iX] = pafLastValue[iX];
+ panThisY[iX] = panLastY[iX];
+ }
+ else
+ {
+ panThisY[iX] = nNoDataVal;
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Write out best index/value to working files. */
+/* -------------------------------------------------------------------- */
+ eErr = GDALRasterIO( hYBand, GF_Write, 0, iY, nXSize, 1,
+ panThisY, nXSize, 1, GDT_UInt32, 0, 0 );
+ if( eErr != CE_None )
+ break;
+
+ eErr = GDALRasterIO( hValBand, GF_Write, 0, iY, nXSize, 1,
+ pafThisValue, nXSize, 1, GDT_Float32, 0, 0 );
+ if( eErr != CE_None )
+ break;
+
+/* -------------------------------------------------------------------- */
+/* Flip this/last buffers. */
+/* -------------------------------------------------------------------- */
+ {
+ float *pafTmp = pafThisValue;
+ pafThisValue = pafLastValue;
+ pafLastValue = pafTmp;
+
+ GUInt32 *panTmp = panThisY;
+ panThisY = panLastY;
+ panLastY = panTmp;
+ }
+
+/* -------------------------------------------------------------------- */
+/* report progress. */
+/* -------------------------------------------------------------------- */
+ if( eErr == CE_None
+ && !pfnProgress( dfProgressRatio * (0.5*(iY+1) / (double)nYSize),
+ "Filling...", pProgressArg ) )
+ {
+ CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+ eErr = CE_Failure;
+ }
+ }
+
+/* ==================================================================== */
+/* Now we will do collect similar this/last information from */
+/* bottom to top and use it in combination with the top to */
+/* bottom search info to interpolate. */
+/* ==================================================================== */
+ for( iY = nYSize-1; iY >= 0 && eErr == CE_None; iY-- )
+ {
+ eErr =
+ GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1,
+ pabyMask, nXSize, 1, GDT_Byte, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+ eErr =
+ GDALRasterIO( hTargetBand, GF_Read, 0, iY, nXSize, 1,
+ pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+/* -------------------------------------------------------------------- */
+/* Figure out the most recent pixel for each column. */
+/* -------------------------------------------------------------------- */
+
+ for( iX = 0; iX < nXSize; iX++ )
+ {
+ if( pabyMask[iX] )
+ {
+ pafThisValue[iX] = pafScanline[iX];
+ panThisY[iX] = iY;
+ }
+ else if( panLastY[iX] - iY <= dfMaxSearchDist )
+ {
+ pafThisValue[iX] = pafLastValue[iX];
+ panThisY[iX] = panLastY[iX];
+ }
+ else
+ {
+ panThisY[iX] = nNoDataVal;
+ }
+ }
+
+/* -------------------------------------------------------------------- */
+/* Load the last y and corresponding value from the top down pass. */
+/* -------------------------------------------------------------------- */
+ eErr =
+ GDALRasterIO( hYBand, GF_Read, 0, iY, nXSize, 1,
+ panTopDownY, nXSize, 1, GDT_UInt32, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+ eErr =
+ GDALRasterIO( hValBand, GF_Read, 0, iY, nXSize, 1,
+ pafTopDownValue, nXSize, 1, GDT_Float32, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+/* -------------------------------------------------------------------- */
+/* Attempt to interpolate any pixels that are nodata. */
+/* -------------------------------------------------------------------- */
+ memset( pabyFiltMask, 0, nXSize );
+ for( iX = 0; iX < nXSize; iX++ )
+ {
+ int iStep, iQuad;
+ int nThisMaxSearchDist = nMaxSearchDist;
+
+ // If this was a valid target - no change.
+ if( pabyMask[iX] )
+ continue;
+
+ // Quadrants 0:topleft, 1:bottomleft, 2:topright, 3:bottomright
+ double adfQuadDist[4];
+ double adfQuadValue[4];
+
+ for( iQuad = 0; iQuad < 4; iQuad++ )
+ {
+ adfQuadDist[iQuad] = dfMaxSearchDist + 1.0;
+ adfQuadValue[iQuad] = 0.0;
+ }
+
+ // Step left and right by one pixel searching for the closest
+ // target value for each quadrant.
+ for( iStep = 0; iStep < nThisMaxSearchDist; iStep++ )
+ {
+ int iLeftX = MAX(0,iX - iStep);
+ int iRightX = MIN(nXSize-1,iX + iStep);
+
+ // top left includes current line
+ QUAD_CHECK(adfQuadDist[0],adfQuadValue[0],
+ iLeftX, panTopDownY[iLeftX], iX, iY,
+ pafTopDownValue[iLeftX] );
+
+ // bottom left
+ QUAD_CHECK(adfQuadDist[1],adfQuadValue[1],
+ iLeftX, panLastY[iLeftX], iX, iY,
+ pafLastValue[iLeftX] );
+
+ // top right and bottom right do no include center pixel.
+ if( iStep == 0 )
+ continue;
+
+ // top right includes current line
+ QUAD_CHECK(adfQuadDist[2],adfQuadValue[2],
+ iRightX, panTopDownY[iRightX], iX, iY,
+ pafTopDownValue[iRightX] );
+
+ // bottom right
+ QUAD_CHECK(adfQuadDist[3],adfQuadValue[3],
+ iRightX, panLastY[iRightX], iX, iY,
+ pafLastValue[iRightX] );
+
+ // every four steps, recompute maximum distance.
+ if( (iStep & 0x3) == 0 )
+ nThisMaxSearchDist = (int) floor(
+ MAX(MAX(adfQuadDist[0],adfQuadDist[1]),
+ MAX(adfQuadDist[2],adfQuadDist[3])) );
+ }
+
+ double dfWeightSum = 0.0;
+ double dfValueSum = 0.0;
+
+ for( iQuad = 0; iQuad < 4; iQuad++ )
+ {
+ if( adfQuadDist[iQuad] <= dfMaxSearchDist )
+ {
+ double dfWeight = 1.0 / adfQuadDist[iQuad];
+
+ dfWeightSum += dfWeight;
+ dfValueSum += adfQuadValue[iQuad] * dfWeight;
+ }
+ }
+
+ if( dfWeightSum > 0.0 )
+ {
+ pabyMask[iX] = 255;
+ pabyFiltMask[iX] = 255;
+ pafScanline[iX] = (float) (dfValueSum / dfWeightSum);
+ }
+
+ }
+
+/* -------------------------------------------------------------------- */
+/* Write out the updated data and mask information. */
+/* -------------------------------------------------------------------- */
+ eErr =
+ GDALRasterIO( hTargetBand, GF_Write, 0, iY, nXSize, 1,
+ pafScanline, nXSize, 1, GDT_Float32, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+ eErr =
+ GDALRasterIO( hFiltMaskBand, GF_Write, 0, iY, nXSize, 1,
+ pabyFiltMask, nXSize, 1, GDT_Byte, 0, 0 );
+
+ if( eErr != CE_None )
+ break;
+
+/* -------------------------------------------------------------------- */
+/* Flip this/last buffers. */
+/* -------------------------------------------------------------------- */
+ {
+ float *pafTmp = pafThisValue;
+ pafThisValue = pafLastValue;
+ pafLastValue = pafTmp;
+
+ GUInt32 *panTmp = panThisY;
+ panThisY = panLastY;
+ panLastY = panTmp;
+ }
+
+/* -------------------------------------------------------------------- */
+/* report progress. */
+/* -------------------------------------------------------------------- */
+ if( eErr == CE_None
+ && !pfnProgress( dfProgressRatio*(0.5+0.5*(nYSize-iY) /
(double)nYSize),
+ "Filling...", pProgressArg ) )
+ {
+ CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
+ eErr = CE_Failure;
+ }
+ }
+
+/* ==================================================================== */
+/* Now we will do iterative average filters over the */
+/* interpolated values to smooth things out and make linear */
+/* artifacts less obvious. */
+/* ==================================================================== */
+ if( eErr == CE_None && nSmoothingIterations > 0 )
+ {
+ // force masks to be to flushed and recomputed.
+ GDALFlushRasterCache( hMaskBand );
+
+ void *pScaledProgress;
+ pScaledProgress =
+ GDALCreateScaledProgress( dfProgressRatio, 1.0, pfnProgress, NULL
);
+
+ eErr = GDALMultiFilter( hTargetBand, hMaskBand, hFiltMaskBand,
+ nSmoothingIterations,
+ GDALScaledProgress, pScaledProgress );
+
+ GDALDestroyScaledProgress( pScaledProgress );
+ }
+
+/* -------------------------------------------------------------------- */
+/* Close and clean up temporary files. Free working buffers */
+/* -------------------------------------------------------------------- */
+end:
+ CPLFree(panLastY);
+ CPLFree(panThisY);
+ CPLFree(panTopDownY);
+ CPLFree(pafLastValue);
+ CPLFree(pafThisValue);
+ CPLFree(pafTopDownValue);
+ CPLFree(pafScanline);
+ CPLFree(pabyMask);
+ CPLFree(pabyFiltMask);
+
+ GDALClose( hYDS );
+ GDALClose( hValDS );
+ GDALClose( hFiltMaskDS );
+
+ CSLDestroy(papszWorkFileOptions);
+
+ GDALDeleteDataset( hDriver, osYTmpFile );
+ GDALDeleteDataset( hDriver, osValTmpFile );
+ GDALDeleteDataset( hDriver, osFiltMaskTmpFile );
+
+ return eErr;
+}
diff --git a/rasterio/rio/features.py b/rasterio/rio/features.py
index f6d4c24..ade9b31 100644
--- a/rasterio/rio/features.py
+++ b/rasterio/rio/features.py
@@ -1,6 +1,8 @@
import functools
import json
import logging
+from math import ceil, floor
+import os
import sys
import warnings
@@ -8,13 +10,15 @@ import click
from cligj import (
precision_opt, indent_opt, compact_opt, projection_geographic_opt,
projection_projected_opt, sequence_opt, use_rs_opt,
- geojson_type_feature_opt, geojson_type_bbox_opt)
+ geojson_type_feature_opt, geojson_type_bbox_opt, files_inout_arg,
+ format_opt)
import rasterio
from rasterio.transform import Affine
from rasterio.rio.cli import cli, coords, write_features
+logger = logging.getLogger('rio')
warnings.simplefilter('default')
@@ -142,3 +146,233 @@ def shapes(
except Exception:
logger.exception("Failed. Exception caught")
sys.exit(1)
+
+
+# Rasterize command.
+@cli.command(short_help='Rasterize features.')
+@files_inout_arg
+@format_opt
+@click.option('--like', type=click.Path(exists=True),
+ help='Raster dataset to use as a template for obtaining affine '
+ 'transform (bounds and resolution), crs, data type, and driver '
+ 'used to create the output. Only a single band will be output.')
+@click.option('--bounds', nargs=4, type=float, default=None,
+ help='Output bounds: left, bottom, right, top.')
+@click.option('--dimensions', nargs=2, type=int, default=None,
+ help='Output dataset width, height in number of pixels.')
+@click.option('--res', multiple=True, type=float, default=None,
+ help='Output dataset resolution in units of coordinate '
+ 'reference system. Pixels assumed to be square if this option is
'
+ 'used once, otherwise use: '
+ '--res pixel_width --res pixel_height')
+@click.option('--src_crs', default='EPSG:4326',
+ help='Source coordinate reference system. Limited to EPSG '
+ 'codes for now. Used as output coordinate system if output does
'
+ 'not exist or --like option is not used. Default: EPSG:4326')
+@click.option('--all_touched', is_flag=True, default=False)
+@click.option('--default_value', type=float, default=1, help='Default value '
+ 'for rasterized pixels')
+@click.option('--fill', type=float, default=0, help='Fill value for all pixels
'
+ 'not overlapping features. Will be evaluated as NoData pixels '
+ 'for output. Default: 0')
+@click.option('--property', type=str, default=None, help='Property in '
+ 'GeoJSON features to use for rasterized values. Any features '
+ 'that lack this property will be given --default_value instead.')
+@click.pass_context
+def rasterize(
+ ctx,
+ files,
+ driver,
+ like,
+ bounds,
+ dimensions,
+ res,
+ src_crs,
+ all_touched,
+ default_value,
+ fill,
+ property):
+
+ """Rasterize GeoJSON into a new or existing raster.
+
+ If the output raster exists, rio-rasterize will rasterize feature values
+ into all bands of that raster. The GeoJSON is assumed to be in the same
+ coordinate reference system as the output.
+
+ --default_value or property values when using --property must be using a
+ data type valid for the data type of that raster.
+
+
+ If a template raster is provided using the --like option, the affine
+ transform, coordinate reference system, and data type from that raster will
+ be used to create the output.
+
+ --default_value or property values when using --property must be using a
+ data type valid for the data type of that raster.
+
+ --driver, --bounds, --dimensions, --res, and --src_crs are ignored when
+ output exists or --like raster is provided
+
+ The GeoJSON must be within the bounds of the existing output or --like
+ raster, or an error will be raised.
+
+
+ If the output does not exist and --like raster is not provided, the input
+ GeoJSON will be used to determine the bounds of the output unless
+ provided using --bounds.
+
+ --dimensions or --res are required in this case.
+
+ If --res is provided, the bottom and right coordinates of bounds are
+ ignored.
+
+
+ Note:
+ The GeoJSON is not projected to match the coordinate reference system
+ of the output or --like rasters at this time. This functionality may be
+ added in the future.
+ """
+
+ import numpy as np
+ from rasterio.features import rasterize
+ from rasterio.features import bounds as calculate_bounds
+
+ verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+
+ files = list(files)
+ output = files.pop()
+ input = click.open_file(files.pop(0) if files else '-')
+
+ # If values are actually meant to be integers, we need to cast them
+ # as such or rasterize creates floating point outputs
+ if default_value == int(default_value):
+ default_value = int(default_value)
+ if fill == int(fill):
+ fill = int(fill)
+
+ with rasterio.drivers(CPL_DEBUG=verbosity > 2):
+
+ def feature_value(feature):
+ if property and 'properties' in feature:
+ return feature['properties'].get(property, default_value)
+ return default_value
+
+ def disjoint_bounds(bounds1, bounds2):
+ return (bounds1[0] > bounds2[2] or
+ bounds1[2] < bounds2[0] or
+ bounds1[1] > bounds2[3] or
+ bounds1[3] < bounds2[1])
+
+ geojson = json.loads(input.read())
+ if 'features' in geojson:
+ geometries = []
+ for f in geojson['features']:
+ geometries.append((f['geometry'], feature_value(f)))
+ elif 'geometry' in geojson:
+ geometries = ((geojson['geometry'], feature_value(geojson)), )
+ else:
+ raise click.BadParameter('Invalid GeoJSON', param=input,
+ param_hint='input')
+
+ geojson_bounds = geojson.get('bbox', calculate_bounds(geojson))
+
+ if os.path.exists(output):
+ with rasterio.open(output, 'r+') as out:
+ if disjoint_bounds(geojson_bounds, out.bounds):
+ raise click.BadParameter('GeoJSON outside bounds of '
+ 'existing output raster',
+ param=input, param_hint='input')
+
+ meta = out.meta.copy()
+
+ result = rasterize(
+ geometries,
+ out_shape=(meta['height'], meta['width']),
+ transform=meta.get('affine', meta['transform']),
+ all_touched=all_touched,
+ dtype=meta.get('dtype', None),
+ default_value=default_value,
+ fill = fill)
+
+ for bidx in range(1, meta['count'] + 1):
+ data = out.read_band(bidx, masked=True)
+ # Burn in any non-fill pixels, and update mask accordingly
+ ne = result != fill
+ data[ne] = result[ne]
+ data.mask[ne] = False
+ out.write_band(bidx, data)
+
+ else:
+ if like is not None:
+ template_ds = rasterio.open(like)
+ if disjoint_bounds(geojson_bounds, template_ds.bounds):
+ raise click.BadParameter('GeoJSON outside bounds of '
+ '--like raster', param=input,
+ param_hint='input')
+
+ kwargs = template_ds.meta.copy()
+ kwargs['count'] = 1
+ template_ds.close()
+
+ else:
+ bounds = bounds or geojson_bounds
+
+ if src_crs == 'EPSG:4326':
+ if (bounds[0] < -180 or bounds[2] > 180 or
+ bounds[1] < -80 or bounds[3] > 80):
+ raise click.BadParameter(
+ 'Bounds are beyond the valid extent for
EPSG:4326.',
+ ctx, param=bounds, param_hint='--bounds')
+
+ if dimensions:
+ width, height = dimensions
+ res = (
+ (bounds[2] - bounds[0]) / float(width),
+ (bounds[3] - bounds[1]) / float(height)
+ )
+
+ else:
+ if not res:
+ raise click.BadParameter(
+ 'pixel dimensions are required',
+ ctx, param=res, param_hint='--res')
+ elif len(res) == 1:
+ res = (res[0], res[0])
+
+ width = max(int(ceil((bounds[2] - bounds[0]) /
+ float(res[0]))), 1)
+ height = max(int(ceil((bounds[3] - bounds[1]) /
+ float(res[1]))), 1)
+
+ src_crs = src_crs.upper()
+ if not src_crs.count('EPSG:'):
+ raise click.BadParameter(
+ 'invalid CRS. Must be an EPSG code.',
+ ctx, param=src_crs, param_hint='--src_crs')
+
+ kwargs = {
+ 'count': 1,
+ 'crs': src_crs,
+ 'width': width,
+ 'height': height,
+ 'transform': Affine(res[0], 0, bounds[0], 0, -res[1],
+ bounds[3]),
+ 'driver': driver
+ }
+
+ result = rasterize(
+ geometries,
+ out_shape=(kwargs['height'], kwargs['width']),
+ transform=kwargs.get('affine', kwargs['transform']),
+ all_touched=all_touched,
+ dtype=kwargs.get('dtype', None),
+ default_value=default_value,
+ fill = fill)
+
+ if 'dtype' not in kwargs:
+ kwargs['dtype'] = result.dtype
+
+ kwargs['nodata'] = fill
+
+ with rasterio.open(output, 'w', **kwargs) as out:
+ out.write_band(1, result)
diff --git a/rasterio/rio/info.py b/rasterio/rio/info.py
index 9c34a9d..427f672 100644
--- a/rasterio/rio/info.py
+++ b/rasterio/rio/info.py
@@ -61,18 +61,29 @@ def env(ctx, key):
@click.option('--bounds', 'meta_member', flag_value='bounds',
help="Print the boundary coordinates "
"(left, bottom, right, top).")
+@click.option('--res', 'meta_member', flag_value='res',
+ help="Print pixel width and height.")
+@click.option('--lnglat', 'meta_member', flag_value='lnglat',
+ help="Print longitude and latitude at center.")
+@click.option('--stats', 'meta_member', flag_value='stats',
+ help="Print statistics (min, max, mean) of a single band "
+ "(use --bidx).")
+@click.option('-v', '--tell-me-more', '--verbose', is_flag=True,
+ help="Output extra information.")
+@click.option('--bidx', type=int, default=1,
+ help="Input file band index (default: 1).")
@click.pass_context
-def info(ctx, input, aspect, indent, namespace, meta_member):
+def info(ctx, input, aspect, indent, namespace, meta_member, verbose, bidx):
"""Print metadata about the dataset as JSON.
Optionally print a single metadata item as a string.
"""
verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
logger = logging.getLogger('rio')
- stdout = click.get_text_stream('stdout')
+ mode = 'r' if (verbose or meta_member == 'stats') else 'r-'
try:
with rasterio.drivers(CPL_DEBUG=(verbosity > 2)):
- with rasterio.open(input, 'r-') as src:
+ with rasterio.open(input, mode) as src:
info = src.meta
del info['affine']
del info['transform']
@@ -82,19 +93,30 @@ def info(ctx, input, aspect, indent, namespace,
meta_member):
if proj4.startswith('+init=epsg'):
proj4 = proj4.split('=')[1].upper()
info['crs'] = proj4
+ info['res'] = src.res
+ info['lnglat'] = src.lnglat()
+ if verbose:
+ stats = [{'min': float(b.min()),
+ 'max': float(b.max()),
+ 'mean': float(b.mean())} for b in src.read()]
+ info['stats'] = stats
if aspect == 'meta':
- if meta_member:
+ if meta_member == 'stats':
+ band = src.read(bidx)
+ click.echo('%f %f %f' % (
+ float(band.min()),
+ float(band.max()),
+ float(band.mean())))
+ elif meta_member:
if isinstance(info[meta_member], (list, tuple)):
- print(" ".join(map(str, info[meta_member])))
+ click.echo(" ".join(map(str, info[meta_member])))
else:
- print(info[meta_member])
+ click.echo(info[meta_member])
else:
- stdout.write(json.dumps(info, indent=indent))
- stdout.write("\n")
+ click.echo(json.dumps(info, indent=indent))
elif aspect == 'tags':
- stdout.write(json.dumps(src.tags(ns=namespace),
+ click.echo(json.dumps(src.tags(ns=namespace),
indent=indent))
- stdout.write("\n")
sys.exit(0)
except Exception:
logger.exception("Failed. Exception caught")
diff --git a/rasterio/rio/main.py b/rasterio/rio/main.py
index 1bccb93..bd17c85 100644
--- a/rasterio/rio/main.py
+++ b/rasterio/rio/main.py
@@ -2,7 +2,8 @@
from rasterio.rio.cli import cli
from rasterio.rio.bands import stack
-from rasterio.rio.features import shapes
+from rasterio.rio.features import shapes, rasterize
from rasterio.rio.info import env, info
from rasterio.rio.merge import merge
from rasterio.rio.rio import bounds, insp, transform
+from rasterio.rio.sample import sample
diff --git a/rasterio/rio/sample.py b/rasterio/rio/sample.py
new file mode 100644
index 0000000..fd4684a
--- /dev/null
+++ b/rasterio/rio/sample.py
@@ -0,0 +1,94 @@
+import json
+import logging
+import sys
+import warnings
+
+import click
+
+import rasterio
+from rasterio.rio.cli import cli
+
+
+warnings.simplefilter('default')
+
+
+@cli.command(short_help="Sample a dataset.")
+@click.argument('files', nargs=-1, required=True, metavar='FILE "[x, y]"')
+@click.option('--bidx', default=None, help="Indexes of input file bands.")
+@click.pass_context
+def sample(ctx, files, bidx):
+ """Sample a dataset at one or more points
+
+ Sampling points (x, y) encoded as JSON arrays, in the coordinate
+ reference system of the dataset, are read from the second
+ positional argument or stdin. Values of the dataset's bands
+ are also encoded as JSON arrays and are written to stdout.
+
+ Example:
+
+ \b
+ $ cat << EOF | rio sample tests/data/RGB.byte.tif
+ > [220650, 2719200]
+ > [219650, 2718200]
+ > EOF
+ [28, 29, 27]
+ [25, 29, 19]
+
+ By default, rio-sample will sample all bands. Optionally, bands
+ may be specified using a simple syntax:
+
+ --bidx N samples the Nth band (first band is 1).
+
+ --bidx M,N,0 samples bands M, N, and O.
+
+ --bidx M..O samples bands M-O, inclusive.
+
+ --bidx ..N samples all bands up to and including N.
+
+ --bidx N.. samples all bands from N to the end.
+
+ Example:
+
+ \b
+ $ cat << EOF | rio sample tests/data/RGB.byte.tif --bidx ..2
+ > [220650, 2719200]
+ > [219650, 2718200]
+ > EOF
+ [28, 29]
+ [25, 29]
+
+ """
+ verbosity = (ctx.obj and ctx.obj.get('verbosity')) or 1
+ logger = logging.getLogger('rio')
+
+ files = list(files)
+ source = files.pop(0)
+ input = files.pop(0) if files else '-'
+
+ # Handle the case of file, stream, or string input.
+ try:
+ points = click.open_file(input).readlines()
+ except IOError:
+ points = [input]
+
+ try:
+ with rasterio.drivers(CPL_DEBUG=verbosity>2):
+ with rasterio.open(source) as src:
+ if bidx is None:
+ indexes = src.indexes
+ elif '..' in bidx:
+ start, stop = map(
+ lambda x: int(x) if x else None, bidx.split('..'))
+ if start is None:
+ start = 1
+ indexes = src.indexes[slice(start-1, stop)]
+ else:
+ indexes = list(map(int, bidx.split(',')))
+ for vals in src.sample(
+ (json.loads(line) for line in points),
+ indexes=indexes):
+ click.echo(json.dumps(vals.tolist()))
+ sys.exit(0)
+ except Exception:
+ logger.exception("Failed. Exception caught")
+ sys.exit(1)
diff --git a/rasterio/warp.py b/rasterio/warp.py
index 25fd5a7..dd971d8 100644
--- a/rasterio/warp.py
+++ b/rasterio/warp.py
@@ -1,6 +1,7 @@
"""Raster warping and reprojection"""
-from rasterio._warp import _reproject, _transform, _transform_geom, RESAMPLING
+from rasterio._base import _transform
+from rasterio._warp import _transform_geom, _reproject, RESAMPLING
from rasterio.transform import guard_transform
diff --git a/requirements-dev.txt b/requirements-dev.txt
index fb826ea..4d9ff53 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -6,5 +6,5 @@ delocate
enum34
numpy>=1.8.0
pytest
-setuptools
+setuptools>=0.9.8
wheel
diff --git a/setup.py b/setup.py
index cd491cf..1788e52 100755
--- a/setup.py
+++ b/setup.py
@@ -11,16 +11,21 @@
import logging
import os
+import pprint
import shutil
import subprocess
import sys
-from setuptools import setup
-from distutils.extension import Extension
+from setuptools import setup
+from setuptools.extension import Extension
logging.basicConfig()
log = logging.getLogger()
+# python -W all setup.py ...
+if 'all' in sys.warnoptions:
+ log.level = logging.DEBUG
+
# Parse the version from the fiona module.
with open('rasterio/__init__.py') as f:
for line in f:
@@ -87,10 +92,9 @@ try:
except Exception as e:
log.warning("Failed to get options via gdal-config: %s", str(e))
-# Conditionally copy PROJ.4 data. Presumes PROJ.4 is installed locally
-# with --prefix=/usr/local.
+# Conditionally copy PROJ.4 data.
if os.environ.get('PACKAGE_DATA'):
- projdatadir = '/usr/local/share/proj'
+ projdatadir = os.environ.get('PROJ_LIB', '/usr/local/share/proj')
if os.path.exists(projdatadir):
try:
shutil.rmtree('rasterio/proj_data')
@@ -104,6 +108,8 @@ ext_options = dict(
libraries=libraries,
extra_link_args=extra_link_args)
+log.debug('ext_options:\n%s', pprint.pformat(ext_options))
+
# When building from a repo, Cython is required.
if os.path.exists("MANIFEST.in") and "clean" not in sys.argv:
log.info("MANIFEST.in found, presume a repo, cythonizing...")
@@ -126,10 +132,12 @@ if os.path.exists("MANIFEST.in") and "clean" not in
sys.argv:
Extension(
'rasterio._warp', ['rasterio/_warp.pyx'], **ext_options),
Extension(
+ 'rasterio._fill', ['rasterio/_fill.pyx',
'rasterio/rasterfill.cpp'], **ext_options),
+ Extension(
'rasterio._err', ['rasterio/_err.pyx'], **ext_options),
Extension(
'rasterio._example', ['rasterio/_example.pyx'], **ext_options),
- ])
+ ], quiet=True)
# If there's no manifest template, as in an sdist, we just specify .c files.
else:
@@ -147,6 +155,8 @@ else:
Extension(
'rasterio._warp', ['rasterio/_warp.cpp'], **ext_options),
Extension(
+ 'rasterio._fill', ['rasterio/_fill.cpp',
'rasterio/rasterfill.cpp'], **ext_options),
+ Extension(
'rasterio._err', ['rasterio/_err.c'], **ext_options),
Extension(
'rasterio._example', ['rasterio/_example.c'], **ext_options),
diff --git a/tests/conftest.py b/tests/conftest.py
index e55a25e..c5a9933 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -4,6 +4,8 @@ import os
import sys
import pytest
+from click.testing import CliRunner
+
if sys.version_info > (3,):
reduce = functools.reduce
@@ -11,6 +13,7 @@ if sys.version_info > (3,):
test_files = [os.path.join(os.path.dirname(__file__), p) for p in [
'data/RGB.byte.tif', 'data/float.tif', 'data/float_nan.tif',
'data/shade.tif']]
+
def pytest_cmdline_main(config):
# Bail if the test raster data is not present. Test data is not
# distributed with sdists since 0.12.
@@ -19,3 +22,8 @@ def pytest_cmdline_main(config):
else:
print("Test data not present. See download directions in
tests/README.txt")
sys.exit(1)
+
+
+@pytest.fixture(scope='function')
+def runner():
+ return CliRunner()
diff --git a/tests/test_features_bounds.py b/tests/test_features_bounds.py
new file mode 100644
index 0000000..a7015cd
--- /dev/null
+++ b/tests/test_features_bounds.py
@@ -0,0 +1,65 @@
+from rasterio.features import bounds
+
+
+# Tests copied from Fiona 1.4.1
+
+def test_bounds_point():
+ g = {'type': 'Point', 'coordinates': [10, 10]}
+ assert bounds(g) == (10, 10, 10, 10)
+
+
+def test_bounds_line():
+ g = {'type': 'LineString', 'coordinates': [[0, 0], [10, 10]]}
+ assert bounds(g) == (0, 0, 10, 10)
+
+
+def test_bounds_polygon():
+ g = {'type': 'Polygon', 'coordinates': [[[0, 0], [10, 10], [10, 0]]]}
+ assert bounds(g) == (0, 0, 10, 10)
+
+
+def test_bounds_z():
+ g = {'type': 'Point', 'coordinates': [10, 10, 10]}
+ assert bounds(g) == (10, 10, 10, 10)
+
+
+# TODO: add these to Fiona with update to bounds
+def test_bounds_existing_bbox():
+ """ Test with existing bbox in geojson, similar to that produced by
+ rasterio. Values specifically modified here for testing, bboxes are not
+ valid as written.
+ """
+
+ fc = {
+ 'bbox': [-107, 40, -105, 41],
+ 'features': [{
+ 'bbox': [-107, 40, -104, 42],
+ 'geometry': {
+ 'coordinates': [
+ [[-107, 40], [-106, 40], [-106, 41], [-107, 41], [-107,
40]]
+ ],
+ 'type': 'Polygon'
+ },
+ 'type': 'Feature'
+ }],
+ 'type': 'FeatureCollection'
+ }
+ assert bounds(fc['features'][0]) == (-107, 40, -104, 42)
+ assert bounds(fc) == (-107, 40, -105, 41)
+
+
+def test_feature_collection():
+ fc = {
+ 'features': [{
+ 'geometry': {
+ 'coordinates': [
+ [[-107, 40], [-106, 40], [-106, 41], [-107, 41], [-107,
40]]
+ ],
+ 'type': 'Polygon'
+ },
+ 'type': 'Feature'
+ }],
+ 'type': 'FeatureCollection'
+ }
+ assert bounds(fc['features'][0]) == (-107, 40, -106, 41)
+ assert bounds(fc) == (-107, 40, -106, 41)
diff --git a/tests/test_fillnodata.py b/tests/test_fillnodata.py
new file mode 100644
index 0000000..2ccee1e
--- /dev/null
+++ b/tests/test_fillnodata.py
@@ -0,0 +1,45 @@
+import logging
+import sys
+import numpy
+import pytest
+
+import rasterio
+from rasterio.fill import fillnodata
+
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+def test_fillnodata():
+ """Test filling nodata values in an ndarray"""
+ # create a 5x5 array, with some missing data
+ a = numpy.ones([3, 3]) * 42
+ a[1][1] = 0
+ # find the missing data
+ mask = ~(a == 0)
+ # fill the missing data using interpolation from the edges
+ result = fillnodata(a, mask)
+ assert(numpy.all((numpy.ones([3, 3]) * 42) == result))
+
+def test_fillnodata_invalid_types():
+ a = numpy.ones([3, 3])
+ with pytest.raises(ValueError):
+ fillnodata(None, a)
+ with pytest.raises(ValueError):
+ fillnodata(a, 42)
+
+def test_fillnodata_mask_ones():
+ # when mask is all ones, image should be unmodified
+ a = numpy.ones([3, 3]) * 42
+ a[1][1] = 0
+ mask = numpy.ones([3, 3])
+ result = fillnodata(a, mask)
+ assert(numpy.all(a == result))
+
+'''
+def test_fillnodata_smooth():
+ a = numpy.array([[1,3,3,1],[2,0,0,2],[2,0,0,2],[1,3,3,1]],
dtype=numpy.float64)
+ mask = ~(a == 0)
+ result = fillnodata(a, mask, max_search_distance=1, smoothing_iterations=0)
+ assert(result[1][1] == 3)
+ result = fillnodata(a, mask, max_search_distance=1, smoothing_iterations=1)
+ assert(round(result[1][1], 1) == 2.2)
+'''
diff --git a/tests/test_rio_features.py b/tests/test_rio_features.py
index c886515..e524cd7 100644
--- a/tests/test_rio_features.py
+++ b/tests/test_rio_features.py
@@ -1,14 +1,81 @@
import logging
+import os
import re
import sys
import click
from click.testing import CliRunner
+import rasterio
from rasterio.rio import features
-logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+# logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+TEST_FEATURES = """{
+ "geometry": {
+ "coordinates": [
+ [
+ [-110, 40],
+ [-100, 40],
+ [-100, 45],
+ [-105, 45],
+ [-110, 40]
+ ]
+ ],
+ "type": "Polygon"
+ },
+ "properties": {
+ "val": 15
+ },
+ "type": "Feature"
+}"""
+
+
+# > rio shapes tests/data/shade.tif --mask --sampling 500 --projected
--precision 0
+TEST_MERC_FEATURECOLLECTION = """{
+ "bbox": [-11858135.0, 4803914.0, -11848351.0, 4813698.0],
+ "features": [{
+ "bbox": [-11853357.504145855, 4808920.97837715,
+ -11848580.189878704, 4813698.2926443005],
+ "geometry": {
+ "coordinates": [
+ [
+ [-11853357.504145855, 4813698.2926443005],
+ [-11853357.504145855, 4808920.97837715],
+ [-11848580.189878704, 4808920.97837715],
+ [-11848580.189878704, 4813698.2926443005],
+ [-11853357.504145855, 4813698.2926443005]
+ ]
+ ],
+ "type": "Polygon"
+ },
+ "properties": {
+ "val": 2
+ },
+ "type": "Feature"
+ }, {
+ "bbox": [-11858134.818413004, 4804143.66411,
+ -11853357.504145855, 4808920.97837715],
+ "geometry": {
+ "coordinates": [
+ [
+ [-11858134.818413004, 4808920.97837715],
+ [-11858134.818413004, 4804143.66411],
+ [-11853357.504145855, 4804143.66411],
+ [-11853357.504145855, 4808920.97837715],
+ [-11858134.818413004, 4808920.97837715]
+ ]
+ ],
+ "type": "Polygon"
+ },
+ "properties": {
+ "val": 3
+ },
+ "type": "Feature"
+ }],
+ "type": "FeatureCollection"
+}"""
def test_err():
@@ -18,24 +85,21 @@ def test_err():
assert result.exit_code == 1
-def test_shapes():
- runner = CliRunner()
+def test_shapes(runner):
result = runner.invoke(features.shapes, ['tests/data/shade.tif'])
assert result.exit_code == 0
assert result.output.count('"FeatureCollection"') == 1
assert result.output.count('"Feature"') == 232
-def test_shapes_sequence():
- runner = CliRunner()
+def test_shapes_sequence(runner):
result = runner.invoke(features.shapes, ['tests/data/shade.tif',
'--sequence'])
assert result.exit_code == 0
assert result.output.count('"FeatureCollection"') == 0
assert result.output.count('"Feature"') == 232
-def test_shapes_sequence_rs():
- runner = CliRunner()
+def test_shapes_sequence_rs(runner):
result = runner.invoke(
features.shapes, [
'tests/data/shade.tif',
@@ -47,24 +111,21 @@ def test_shapes_sequence_rs():
assert result.output.count(u'\u001e') == 232
-def test_shapes_with_nodata():
- runner = CliRunner()
+def test_shapes_with_nodata(runner):
result = runner.invoke(features.shapes, ['tests/data/shade.tif',
'--with-nodata'])
assert result.exit_code == 0
assert result.output.count('"FeatureCollection"') == 1
assert result.output.count('"Feature"') == 288
-def test_shapes_indent():
- runner = CliRunner()
+def test_shapes_indent(runner):
result = runner.invoke(features.shapes, ['tests/data/shade.tif',
'--indent', '2'])
assert result.exit_code == 0
assert result.output.count('"FeatureCollection"') == 1
assert result.output.count('\n') == 70139
-def test_shapes_compact():
- runner = CliRunner()
+def test_shapes_compact(runner):
result = runner.invoke(features.shapes, ['tests/data/shade.tif',
'--compact'])
assert result.exit_code == 0
assert result.output.count('"FeatureCollection"') == 1
@@ -72,8 +133,7 @@ def test_shapes_compact():
assert result.output.count(': ') == 0
-def test_shapes_sampling():
- runner = CliRunner()
+def test_shapes_sampling(runner):
result = runner.invoke(
features.shapes, ['tests/data/shade.tif', '--sampling', '10'])
assert result.exit_code == 0
@@ -81,8 +141,7 @@ def test_shapes_sampling():
assert result.output.count('"Feature"') == 124
-def test_shapes_precision():
- runner = CliRunner()
+def test_shapes_precision(runner):
result = runner.invoke(
features.shapes, ['tests/data/shade.tif', '--precision', '1'])
assert result.exit_code == 0
@@ -91,9 +150,194 @@ def test_shapes_precision():
assert re.search(r'\d*\.\d{2,}', result.output) is None
-def test_shapes_mask():
- runner = CliRunner()
+def test_shapes_mask(runner):
result = runner.invoke(features.shapes, ['tests/data/RGB.byte.tif',
'--mask'])
assert result.exit_code == 0
assert result.output.count('"FeatureCollection"') == 1
assert result.output.count('"Feature"') == 9
+
+
+def test_rasterize_err(tmpdir, runner):
+ output = str(tmpdir.join('test.tif'))
+ # Test invalid stdin
+ result = runner.invoke(features.rasterize, [output], input='BOGUS')
+ assert result.exit_code == -1
+
+ # Test invalid GeoJSON
+ result = runner.invoke(features.rasterize, [output],
+ input='{"foo": "bar"}')
+ assert result.exit_code == 2
+
+ # Test invalid res
+ result = runner.invoke(features.rasterize, [output], input=TEST_FEATURES)
+ assert result.exit_code == 2
+
+ # Test invalid CRS for bounds
+ result = runner.invoke(features.rasterize, [output, '--res', 1],
+ input=TEST_MERC_FEATURECOLLECTION)
+ assert result.exit_code == 2
+
+ # Test invalid CRS value
+ result = runner.invoke(features.rasterize, [output,
+ '--res', 1,
+ '--src_crs', 'BOGUS'],
+ input=TEST_MERC_FEATURECOLLECTION)
+ assert result.exit_code == 2
+
+
+def test_rasterize(tmpdir, runner):
+ # Test dimensions
+ output = str(tmpdir.join('test.tif'))
+ result = runner.invoke(features.rasterize,
+ [output, '--dimensions', 20, 10],
+ input=TEST_FEATURES)
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+ with rasterio.open(output) as out:
+ assert out.count == 1
+ assert out.meta['width'] == 20
+ assert out.meta['height'] == 10
+ data = out.read_band(1, masked=False)
+ assert (data == 0).sum() == 55
+ assert (data == 1).sum() == 145
+
+ # Test dimensions and bounds
+ output = str(tmpdir.join('test2.tif'))
+ result = runner.invoke(features.rasterize,
+ [output,
+ '--dimensions', 40, 20,
+ '--bounds', -120, 30, -90, 50
+ ], input=TEST_FEATURES)
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+ with rasterio.open(output) as out:
+ assert out.count == 1
+ assert out.meta['width'] == 40
+ assert out.meta['height'] == 20
+ data = out.read_band(1, masked=False)
+ assert (data == 0).sum() == 748
+ assert (data == 1).sum() == 52
+
+ # Test resolution
+ output = str(tmpdir.join('test3.tif'))
+ result = runner.invoke(features.rasterize,
+ [output, '--res', 0.5], input=TEST_FEATURES)
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+ with rasterio.open(output) as out:
+ assert out.count == 1
+ assert out.meta['width'] == 20
+ assert out.meta['height'] == 10
+ data = out.read_band(1, masked=False)
+ assert (data == 0).sum() == 55
+ assert (data == 1).sum() == 145
+
+
+def test_rasterize_existing_output(tmpdir, runner):
+ output = str(tmpdir.join('test.tif'))
+ result = runner.invoke(features.rasterize,
+ [output, '--res', 0.5], input=TEST_FEATURES)
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+
+ geojson = """{
+ "geometry": {
+ "coordinates": [
+ [
+ [-102, 40],
+ [-98, 40],
+ [-98, 45],
+ [-100, 45],
+ [-102, 40]
+ ]
+ ],
+ "type": "Polygon"
+ },
+ "type": "Feature"
+ }"""
+
+ result = runner.invoke(features.rasterize, [output, '--default_value', 2],
+ input=geojson)
+
+ with rasterio.open(output) as out:
+ assert out.count == 1
+ data = out.read_band(1, masked=False)
+ assert (data == 0).sum() == 55
+ assert (data == 1).sum() == 125
+ assert (data == 2).sum() == 20
+
+
+def test_rasterize_like(tmpdir, runner):
+ output = str(tmpdir.join('test.tif'))
+ result = runner.invoke(features.rasterize,
+ [output, '--like', 'tests/data/shade.tif'],
+ input=TEST_MERC_FEATURECOLLECTION)
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+ with rasterio.open(output) as out:
+ assert out.count == 1
+ data = out.read_band(1, masked=False)
+ assert (data == 0).sum() == 548576
+ assert (data == 1).sum() == 500000
+
+ # Test invalid like raster
+ output = str(tmpdir.join('test2.tif'))
+ result = runner.invoke(features.rasterize,
+ [output, '--like', 'foo.tif'], input=TEST_FEATURES)
+ assert result.exit_code == 2
+
+
+def test_rasterize_property_value(tmpdir, runner):
+ # Test feature collection property values
+ output = str(tmpdir.join('test.tif'))
+ result = runner.invoke(features.rasterize,
+ [output,
+ '--res', 1000,
+ '--property', 'val',
+ '--src_crs', 'EPSG:3857'
+ ],
+ input=TEST_MERC_FEATURECOLLECTION)
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+ with rasterio.open(output) as out:
+ assert out.count == 1
+ data = out.read_band(1, masked=False)
+ assert (data == 0).sum() == 50
+ assert (data == 2).sum() == 25
+ assert (data == 3).sum() == 25
+
+ # Test feature property values
+ output = str(tmpdir.join('test2.tif'))
+ result = runner.invoke(features.rasterize,
+ [output, '--res', 0.5, '--property', 'val'],
+ input=TEST_FEATURES)
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+ with rasterio.open(output) as out:
+ assert out.count == 1
+ data = out.read_band(1, masked=False)
+ assert (data == 0).sum() == 55
+ assert (data == 15).sum() == 145
+
+def test_rasterize_out_of_bounds(tmpdir, runner):
+ output = str(tmpdir.join('test.tif'))
+
+ # Test out of bounds of --like raster
+ result = runner.invoke(features.rasterize,
+ [output, '--like', 'tests/data/shade.tif'],
+ input=TEST_FEATURES)
+ assert result.exit_code == 2
+
+ # Test out of bounds of existing output raster (first have to create one)
+ result = runner.invoke(features.rasterize,
+ [output,
+ '--res', 1000,
+ '--property', 'val',
+ '--src_crs', 'EPSG:3857'
+ ],
+ input=TEST_MERC_FEATURECOLLECTION)
+ assert result.exit_code == 0
+ assert os.path.exists(output)
+
+ result = runner.invoke(features.rasterize, [output], input=TEST_FEATURES)
+ assert result.exit_code == 2
diff --git a/tests/test_rio_info.py b/tests/test_rio_info.py
index fce04eb..0dd5c65 100644
--- a/tests/test_rio_info.py
+++ b/tests/test_rio_info.py
@@ -71,3 +71,45 @@ def test_info_tags():
['tests/data/RGB.byte.tif', '--tags'])
assert result.exit_code == 0
assert result.output == '{"AREA_OR_POINT": "Area"}\n'
+
+
+def test_info_res():
+ runner = CliRunner()
+ result = runner.invoke(
+ info.info,
+ ['tests/data/RGB.byte.tif', '--res'])
+ assert result.exit_code == 0
+ assert result.output.startswith('300.037')
+
+
+def test_info_lnglat():
+ runner = CliRunner()
+ result = runner.invoke(
+ info.info,
+ ['tests/data/RGB.byte.tif', '--lnglat'])
+ assert result.exit_code == 0
+ assert result.output.startswith('-77.757')
+
+
+def test_mo_info():
+ runner = CliRunner()
+ result = runner.invoke(info.info, ['tests/data/RGB.byte.tif'])
+ assert result.exit_code == 0
+ assert '"res": [300.037' in result.output
+ assert '"lnglat": [-77.757' in result.output
+
+
+def test_info_stats():
+ runner = CliRunner()
+ result = runner.invoke(info.info, ['tests/data/RGB.byte.tif',
'--tell-me-more'])
+ assert result.exit_code == 0
+ assert '"max": 255.0' in result.output
+ assert '"min": 1.0' in result.output
+ assert '"mean": 44.4344' in result.output
+
+
+def test_info_stats_only():
+ runner = CliRunner()
+ result = runner.invoke(info.info, ['tests/data/RGB.byte.tif', '--stats',
'--bidx', '2'])
+ assert result.exit_code == 0
+ assert result.output.startswith('1.000000 255.000000 66.02')
diff --git a/tests/test_rio_sample.py b/tests/test_rio_sample.py
new file mode 100644
index 0000000..2c6a329
--- /dev/null
+++ b/tests/test_rio_sample.py
@@ -0,0 +1,81 @@
+import logging
+import sys
+
+import click
+from click.testing import CliRunner
+
+import rasterio
+from rasterio.rio import sample
+
+
+logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
+
+
+def test_sample_err():
+ runner = CliRunner()
+ result = runner.invoke(
+ sample.sample,
+ ['bogus.tif'],
+ "[220650.0, 2719200.0]")
+ assert result.exit_code == 1
+
+
+def test_sample_stdin():
+ runner = CliRunner()
+ result = runner.invoke(
+ sample.sample,
+ ['tests/data/RGB.byte.tif'],
+ "[220650.0, 2719200.0]\n[220650.0, 2719200.0]",
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ assert result.output.strip() == '[28, 29, 27]\n[28, 29, 27]'
+
+
+def test_sample_arg():
+ runner = CliRunner()
+ result = runner.invoke(
+ sample.sample,
+ ['tests/data/RGB.byte.tif', "[220650.0, 2719200.0]"],
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ assert result.output.strip() == '[28, 29, 27]'
+
+
+def test_sample_bidx():
+ runner = CliRunner()
+ result = runner.invoke(
+ sample.sample,
+ ['tests/data/RGB.byte.tif', '--bidx', '1,2', "[220650.0, 2719200.0]"],
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ assert result.output.strip() == '[28, 29]'
+
+
+def test_sample_bidx2():
+ runner = CliRunner()
+ result = runner.invoke(
+ sample.sample,
+ ['tests/data/RGB.byte.tif', '--bidx', '1..2', "[220650.0, 2719200.0]"],
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ assert result.output.strip() == '[28, 29]'
+
+
+def test_sample_bidx3():
+ runner = CliRunner()
+ result = runner.invoke(
+ sample.sample,
+ ['tests/data/RGB.byte.tif', '--bidx', '..2', "[220650.0, 2719200.0]"],
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ assert result.output.strip() == '[28, 29]'
+
+
+def test_sample_bidx4():
+ runner = CliRunner()
+ result = runner.invoke(
+ sample.sample,
+ ['tests/data/RGB.byte.tif', '--bidx', '3', "[220650.0, 2719200.0]"],
+ catch_exceptions=False)
+ assert result.exit_code == 0
+ assert result.output.strip() == '[27]'
diff --git a/tests/test_sampling.py b/tests/test_sampling.py
new file mode 100644
index 0000000..c4030cc
--- /dev/null
+++ b/tests/test_sampling.py
@@ -0,0 +1,17 @@
+import rasterio
+
+
+def test_sampling():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ data = next(src.sample([(220650.0, 2719200.0)]))
+ assert list(data) == [28, 29, 27]
+
+def test_sampling_beyond_bounds():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ data = next(src.sample([(-10, 2719200.0)]))
+ assert list(data) == [0, 0, 0]
+
+def test_sampling_indexes():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ data = next(src.sample([(220650.0, 2719200.0)], indexes=[2]))
+ assert list(data) == [29]
diff --git a/tests/test_tags.py b/tests/test_tags.py
index 4431055..514d068 100644
--- a/tests/test_tags.py
+++ b/tests/test_tags.py
@@ -54,3 +54,11 @@ def test_tags_update_twice():
assert dst.tags() == {'a': '1', 'b': '2'}
dst.update_tags(c=3)
assert dst.tags() == {'a': '1', 'b': '2', 'c': '3'}
+
+
+def test_tags_eq():
+ with rasterio.open(
+ 'test.tif', 'w',
+ 'GTiff', 3, 4, 1, dtype=rasterio.ubyte) as dst:
+ dst.update_tags(a="foo=bar")
+ assert dst.tags() == {'a': "foo=bar"}
diff --git a/tests/test_transform.py b/tests/test_transform.py
index b17feab..acd0d54 100644
--- a/tests/test_transform.py
+++ b/tests/test_transform.py
@@ -1,4 +1,3 @@
-
import rasterio
def test_window():
@@ -9,3 +8,16 @@ def test_window():
assert src.window(left, top-src.res[1], left+src.res[0], top) == (
(0, 1), (0, 1))
+
+def test_window_transform():
+ with rasterio.open('tests/data/RGB.byte.tif') as src:
+ assert src.window_transform(((0, None), (0, None))) == src.affine
+ assert src.window_transform(((None, None), (None, None))) == src.affine
+ assert src.window_transform(
+ ((1, None), (1, None))).c == src.bounds.left + src.res[0]
+ assert src.window_transform(
+ ((1, None), (1, None))).f == src.bounds.top - src.res[1]
+ assert src.window_transform(
+ ((-1, None), (-1, None))).c == src.bounds.left - src.res[0]
+ assert src.window_transform(
+ ((-1, None), (-1, None))).f == src.bounds.top + src.res[1]
--
Alioth's /usr/local/bin/git-commit-notice on
/srv/git.debian.org/git/pkg-grass/rasterio.git
_______________________________________________
Pkg-grass-devel mailing list
Pkg-grass-devel@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
