Add RasterIO backend

ci/requirements-py27-cdat+pynio.yml

@@ -18,6 +18,7 @@ dependencies:
   - scipy
   - seaborn
   - toolz
+  - rasterio
   - pip:
     - coveralls
     - pytest-cov

ci/requirements-py27-windows.yml

@@ -15,3 +15,4 @@ dependencies:
   - scipy
   - seaborn
   - toolz
+  - rasterio

ci/requirements-py35.yml

@@ -15,6 +15,7 @@ dependencies:
   - scipy
   - seaborn
   - toolz
+  - rasterio
   - pip:
     - coveralls
     - pytest-cov

ci/requirements-py36-windows.yml

@@ -15,3 +15,4 @@ dependencies:
   - scipy
   - seaborn
   - toolz
+  - rasterio

ci/requirements-py36.yml

@@ -15,6 +15,7 @@ dependencies:
   - scipy
   - seaborn
   - toolz
+  - rasterio
   - pip:
     - coveralls
     - pytest-cov

doc/api.rst

@@ -417,6 +417,7 @@ Dataset methods
 
    open_dataset
    open_mfdataset
+   open_rasterio
    Dataset.to_netcdf
    save_mfdataset
    Dataset.to_array

doc/gallery/plot_rasterio.py

@@ -0,0 +1,75 @@
+# -*- coding: utf-8 -*-
+"""
+.. _recipes.rasterio:
+
+=================================
+Parsing rasterio's geocoordinates
+=================================
+
+
+The example illustrates how to use an accessor (see :ref:`internals.accessors`)
+to  convert a projection's cartesian coordinates into 2D longitudes and
+latitudes.
+
+These new coordinates might be handy for plotting and indexing, but it should
+be kept in mind that a grid which is regular in projection coordinates will
+likely be irregular in lon/lat. It is often recommended to work in the data's
+original map projection.
+"""
+
+import os
+import urllib.request
+import numpy as np
+import xarray as xr
+import cartopy.crs as ccrs
+import matplotlib.pyplot as plt
+from rasterio.warp import transform
+
+
+# Define the accessor
+@xr.register_dataarray_accessor('rasterio')
+class RasterioAccessor(object):
+    def __init__(self, xarray_obj):
+        self._obj = xarray_obj
+
+    def add_lonlat_coords(self):
+        """Compute the lon/lat coordinates out of the dataset's crs.
+
+        This adds two non-dimension coordinates ('lon' and 'lat') to the
+        original dataarray.
+        """
+
+        ny, nx = len(self._obj['y']), len(self._obj['x'])
+        x, y = np.meshgrid(self._obj['x'], self._obj['y'])
+
+        # Rasterio works with 1D arrays
+        lon, lat = transform(self._obj.crs, {'init': 'EPSG:4326'},
+                             x.flatten(), y.flatten())
+        lon = np.asarray(lon).reshape((ny, nx))
+        lat = np.asarray(lat).reshape((ny, nx))
+        self._obj.coords['lon'] = (('y', 'x'), lon)
+        self._obj.coords['lat'] = (('y', 'x'), lat)
+
+
+# Download the file from rasterio's repository
+url = 'https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif'
+urllib.request.urlretrieve(url, 'RGB.byte.tif')
+
+# Read the data
+rioda = xr.open_rasterio('RGB.byte.tif')
+
+# Compute the coordinates
+rioda.rasterio.add_lonlat_coords()
+
+# Compute a greyscale out of the rgb image
+greyscale = rioda.mean(dim='band')
+
+# Plot on a map
+ax = plt.subplot(projection=ccrs.PlateCarree())
+greyscale.plot(ax=ax, x='lon', y='lat', transform=ccrs.PlateCarree(),
+               cmap='Greys_r', add_colorbar=False)
+ax.coastlines('10m', color='r')
+plt.show()
+
+# Delete the file
+os.remove('RGB.byte.tif')

doc/io.rst

@@ -384,6 +384,48 @@ over the network until we look at particular values:
 
 .. image:: _static/opendap-prism-tmax.png
 
+.. _io.rasterio:
+
+Rasterio
+--------
+
+GeoTIFFs and other gridded raster datasets can be opened using `rasterio`_, if
+rasterio is installed. Here is an example of how to use
+:py:func:`~xarray.open_rasterio` to read one of rasterio's `test files`_:
+
+.. ipython::
+    :verbatim:
+
+    In [7]: rio = xr.open_rasterio('RGB.byte.tif')
+
+    In [8]: rio
+    Out[8]:
+    <xarray.DataArray (band: 3, y: 718, x: 791)>
+    [1703814 values with dtype=uint8]
+    Coordinates:
+      * band     (band) int64 1 2 3
+      * y        (y) float64 2.827e+06 2.827e+06 2.826e+06 2.826e+06 2.826e+06 ...
+      * x        (x) float64 1.02e+05 1.023e+05 1.026e+05 1.029e+05 1.032e+05 ...
+    Attributes:
+        crs:      +init=epsg:32618
+
+The ``x`` and ``y`` coordinates are generated out of the file's metadata
+(``bounds``, ``width``, ``height``), and they can be understood as cartesian
+coordinates defined in the file's projection provided by the ``crs`` attribute.
+``crs`` is a PROJ4 string which can be parsed by e.g. `pyproj`_ or rasterio.
+See :ref:`recipes.rasterio` for an example of how to convert these to
+longitudes and latitudes.
+
+.. warning::
+
+    This feature has been added in xarray v0.9.6 and should still be
+    considered as being experimental. Please report any bug you may find
+    on xarray's github repository.
+
+.. _rasterio: https://mapbox.github.io/rasterio/
+.. _test files: https://github.com/mapbox/rasterio/blob/master/tests/data/RGB.byte.tif
+.. _pyproj: https://github.com/jswhit/pyproj
+
 .. _io.pynio:
 
 Formats supported by PyNIO

doc/whats-new.rst

@@ -37,11 +37,17 @@ Enhancements
   controlled via ``--run-network-tests`` command line argument
   to ``py.test`` (:issue:`1393`).
   By `Matthew Gidden <https://github.com/gidden>`_.
-
+ 
 - :py:func:`~xarray.align` now supports ``join='exact'``, which raises
   an error instead of aligning when indexes to be aligned are not equal.
   By `Stephan Hoyer <https://github.com/shoyer>`_.
 
+- New backend to open raster files with the
+  `rasterio <https://mapbox.github.io/rasterio/>`_ library.
+  By `Joe Hamman <https://github.com/jhamman>`_,
+  `Nic Wayand <https://github.com/NicWayand>`_ and
+  `Fabien Maussion <https://github.com/fmaussion>`_
+
 Bug fixes
 ~~~~~~~~~
 

xarray/__init__.py

@@ -16,6 +16,8 @@
 
 from .backends.api import (open_dataset, open_dataarray, open_mfdataset,
                            save_mfdataset)
+from .backends.rasterio_ import open_rasterio
+
 from .conventions import decode_cf
 
 try:

xarray/backends/api.py

@@ -181,11 +181,10 @@ def open_dataset(filename_or_obj, group=None, decode_cf=True,
     chunks : int or dict, optional
         If chunks is provided, it used to load the new dataset into dask
         arrays. ``chunks={}`` loads the dataset with dask using a single
-        chunk for all arrays. This is an experimental feature; see the
-        documentation for more details.
+        chunk for all arrays.
     lock : False, True or threading.Lock, optional
         If chunks is provided, this argument is passed on to
-        :py:func:`dask.array.from_array`. By default, a per-variable lock is
+        :py:func:`dask.array.from_array`. By default, a global lock is
         used when reading data from netCDF files with the netcdf4 and h5netcdf
         engines to avoid issues with concurrent access when using dask's
         multithreaded backend.
@@ -228,17 +227,7 @@ def maybe_decode_store(store, lock=False):
         _protect_dataset_variables_inplace(ds, cache)
 
         if chunks is not None:
-            try:
-                from dask.base import tokenize
-            except ImportError:
-                # raise the usual error if dask is entirely missing
-                import dask
-                if LooseVersion(dask.__version__) < LooseVersion('0.6'):
-                    raise ImportError(
-                        'xarray requires dask version 0.6 or newer')
-                else:
-                    raise
-
+            from dask.base import tokenize
             # if passed an actual file path, augment the token with
             # the file modification time
             if (isinstance(filename_or_obj, basestring) and
@@ -369,11 +358,10 @@ def open_dataarray(*args, **kwargs):
         'netcdf4'.
     chunks : int or dict, optional
         If chunks is provided, it used to load the new dataset into dask
-        arrays. This is an experimental feature; see the documentation for more
-        details.
+        arrays.
     lock : False, True or threading.Lock, optional
         If chunks is provided, this argument is passed on to
-        :py:func:`dask.array.from_array`. By default, a per-variable lock is
+        :py:func:`dask.array.from_array`. By default, a global lock is
         used when reading data from netCDF files with the netcdf4 and h5netcdf
         engines to avoid issues with concurrent access when using dask's
         multithreaded backend.