Repository: pydata/xarray

Issue:
rasterio.open exposes the GeoTIFF scales attribute, but xarray.open_rasterio not
#### Code Sample, a copy-pastable example if possible

```python
>> import xarray as xr
>> from pprint import pprint
>> arr = xr.open_rasterio('example.tif')
>> pprint(arr)
<xarray.DataArray (band: 2, y: 7201, x: 5896)>
[84914192 values with dtype=uint16]
Coordinates:
  * band     (band) int64 1 2
  * y        (y) float64 -17.18 -17.18 -17.18 -17.19 ... -23.58 -23.58 -23.58
  * x        (x) float64 -52.23 -52.23 -52.23 -52.23 ... -46.99 -46.99 -46.99
Attributes:
    transform:   (0.0008888888888890009, 0.0, -52.22933333332605, 0.0, -0.000...
    crs:         +init=epsg:4326
    res:         (0.0008888888888890009, 0.0008888888888890009)
    is_tiled:    1
    nodatavals:  (65535.0, 65535.0)
```

Would be nice if it had the `scales` attribute.

#### Problem description

rasterio.open exposes the GeoTIFF scales attribute, but xarray not.

```python
>> import rasterio
>> rarr = rasterio.open('example.tif')
>> print(rarr.scales)
(0.001, 0.001)
```

#### Expected Output

```python
>> import xarray as xr
>> from pprint import pprint
>> arr = xr.open_rasterio('example.tif')
>> pprint(arr)
<xarray.DataArray (band: 2, y: 7201, x: 5896)>
[84914192 values with dtype=uint16]
Coordinates:
  * band     (band) int64 1 2
  * y        (y) float64 -17.18 -17.18 -17.18 -17.19 ... -23.58 -23.58 -23.58
  * x        (x) float64 -52.23 -52.23 -52.23 -52.23 ... -46.99 -46.99 -46.99
Attributes:
    transform:   (0.0008888888888890009, 0.0, -52.22933333332605, 0.0, -0.000...
    crs:         +init=epsg:4326
    res:         (0.0008888888888890009, 0.0008888888888890009)
    is_tiled:    1
    nodatavals:  (65535.0, 65535.0)
    scales: (0.001, 0.001)
```

#### xarray version

```python
xr.__version__ == '0.12.1'
```


Discussion:
dcherian: Does rasterio define a set of valid attributes or can you assign arbitrary attributes? If so, should we read and assign all attributes as we do for netCDF files?
ecarrara: > Does rasterio define a set of valid attributes or can you assign arbitrary attributes? If so, should we read and assign all attributes as we do for netCDF files?

@dcherian the attributes are pre-defined. Looking at the source code[1] of rasterio I found this list:

```property nodatavals             # Nodata values for each band
property nodata                 # The dataset's single nodata value
property mask_flag_enums        # Sets of flags describing the sources of band masks.
property crs                    # The dataset's coordinate reference system
property descriptions           # Descriptions for each dataset band
property transform              # The dataset's georeferencing transformation matrix
property offsets                # Raster offset for each dataset band
property scales                 # Raster scale for each dataset band
property units                  # A list of str: one units string for each dataset band
property bounds                 # Returns the lower left and upper right bounds of the dataset
property res                    # Returns the (width, height) of pixels in the units of its
```

I belive that the attributes not handled yet by xarray are: **scales, offsets and units**.

[1] - https://github.com/mapbox/rasterio/blob/1.0.24/rasterio/_base.pyx
dcherian: I think we should read and expose all those attributes. 

@ecarrara are you up for expanding the scope of your PR? Also, I see you're a "first-time contributor". Welcome! and thanks for contributing!

ping @fmaussion for feedback.
ecarrara: thank you @dcherian!

I added `offsets`, `units`, `descriptions` attributes. The other attributes (`res`, `nodata`, `crs`,  `transform`, `bounds` and `res`) were already available.

`mask_flag_enums` is a list of tuples, I have tried to add this attribute too but the `TestRasterio::test_serialization` fails with the message:

```python
>   ???
E   ValueError: multi-dimensional array attributes not supported

netCDF4/_netCDF4.pyx:1529: ValueError
```

@fmaussion any feedback is welcome!

pep8speaks: Hello @ecarrara! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:



* In the file [`xarray/tests/test_backends.py`](https://github.com/pydata/xarray/blob/bb34639b359187cb1275ef5e194017d9af9c3351/xarray/tests/test_backends.py):

> [Line 3076:80](https://github.com/pydata/xarray/blob/bb34639b359187cb1275ef5e194017d9af9c3351/xarray/tests/test_backends.py#L3076): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (81 > 79 characters)
> [Line 3121:80](https://github.com/pydata/xarray/blob/bb34639b359187cb1275ef5e194017d9af9c3351/xarray/tests/test_backends.py#L3121): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (81 > 79 characters)
> [Line 3148:80](https://github.com/pydata/xarray/blob/bb34639b359187cb1275ef5e194017d9af9c3351/xarray/tests/test_backends.py#L3148): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (81 > 79 characters)


fmaussion: Thanks @ecarrara !

>  but the TestRasterio::test_serialization fails

This is not the case anymore, right? This test is here to ensure that the dataarrays we are creating with `open_rasterio` are compatible with the NetCDF model. I think it is good to try to keep it that way.
dcherian: Yeah I think @ecarrara found a workaround. 

@ecarrara You should also add a note to whats-new and give yourself credit.
dcherian: OK let's remove the mask_flag_enums bits then. I noticed that there's a binary .idx file that was added to. Is that needed?
shoyer: I fixed up the merge. Let's just make sure tests pass...
codecov[bot]: # [Codecov](https://codecov.io/gh/pydata/xarray/pull/3013?src=pr&el=h1) Report
> Merging [#3013](https://codecov.io/gh/pydata/xarray/pull/3013?src=pr&el=desc) into [master](https://codecov.io/gh/pydata/xarray/commit/681893ed96685052e5f5822cc867245aaf138b32?src=pr&el=desc) will **decrease** coverage by `25.58%`.
> The diff coverage is `0%`.


```diff
@@             Coverage Diff             @@
##           master    #3013       +/-   ##
===========================================
- Coverage   92.16%   66.58%   -25.59%     
===========================================
  Files          63       63               
  Lines       12787    12795        +8     
  Branches     2906     2910        +4     
===========================================
- Hits        11785     8519     -3266     
- Misses        688     3975     +3287     
+ Partials      314      301       -13
```


Patch:
diff --git a/.coveragerc b/.coveragerc
index 40369c31..2bf6e08f 100644
--- a/.coveragerc
+++ b/.coveragerc
@@ -1,5 +1,4 @@
 [run]
-branch = True
 omit =
     xarray/tests/*
     xarray/core/dask_array_compat.py
diff --git a/ci/requirements/py36.yml b/ci/requirements/py36.yml
index 313eafd4..59ce883d 100644
--- a/ci/requirements/py36.yml
+++ b/ci/requirements/py36.yml
@@ -8,7 +8,9 @@ dependencies:
   - distributed
   - h5py
   - h5netcdf
-  - matplotlib
+  # pin matplotlib for now, to ensure that we have at least one CI job
+  # testing pseudonetcdf: https://github.com/barronh/pseudonetcdf/issues/69
+  - matplotlib=3.0
   - netcdf4
   - pytest
   - pytest-cov
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 7d730f51..896f38fb 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -30,12 +30,18 @@ Enhancements
 - Renaming variables and dimensions independently:
   Datasets with coordinate dimensions can now have only their dimension 
   (using rename_dim) or only their coordinate (using rename_vars) renamed 
-  instead of the rename function applyingto both. (:issue:`3026`)
+  instead of the rename function applying to both. (:issue:`3026`)
   By `Julia Kent <https://github.com/jukent>`_.
-  
+
+- Add ``scales``, ``offsets``, ``units`` and ``descriptions``
+  attributes to :py:class:`~xarray.DataArray` returned by
+  :py:func:`~xarray.open_rasterio`. (:issue:`3013`)
+  By `Erle Carrara <https://github.com/ecarrara>`_.
+
 Bug fixes
 ~~~~~~~~~
 
+- Resolved deprecation warnings from newer versions of matplotlib and dask.
 - Compatibility fixes for the upcoming pandas 0.25 and NumPy 1.17 releases.
   By `Stephan Hoyer <https://github.com/shoyer>`_.
 
diff --git a/xarray/backends/common.py b/xarray/backends/common.py
index a433c805..6f0859f6 100644
--- a/xarray/backends/common.py
+++ b/xarray/backends/common.py
@@ -122,7 +122,7 @@ class AbstractDataStore(Mapping):
         return variables, attributes
 
     @property
-    def variables(self):
+    def variables(self):  # pragma: no cover
         warnings.warn('The ``variables`` property has been deprecated and '
                       'will be removed in xarray v0.11.',
                       FutureWarning, stacklevel=2)
@@ -130,7 +130,7 @@ class AbstractDataStore(Mapping):
         return variables
 
     @property
-    def attrs(self):
+    def attrs(self):  # pragma: no cover
         warnings.warn('The ``attrs`` property has been deprecated and '
                       'will be removed in xarray v0.11.',
                       FutureWarning, stacklevel=2)
@@ -138,7 +138,7 @@ class AbstractDataStore(Mapping):
         return attrs
 
     @property
-    def dimensions(self):
+    def dimensions(self):  # pragma: no cover
         warnings.warn('The ``dimensions`` property has been deprecated and '
                       'will be removed in xarray v0.11.',
                       FutureWarning, stacklevel=2)
diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py
index 2f00702f..cc533cc3 100644
--- a/xarray/backends/rasterio_.py
+++ b/xarray/backends/rasterio_.py
@@ -301,6 +301,18 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None,
         attrs['nodatavals'] = tuple(
             np.nan if nodataval is None else nodataval
             for nodataval in riods.nodatavals)
+    if hasattr(riods, 'scales'):
+        # The scale values for the raster bands
+        attrs['scales'] = riods.scales
+    if hasattr(riods, 'offsets'):
+        # The offset values for the raster bands
+        attrs['offsets'] = riods.offsets
+    if hasattr(riods, 'descriptions') and any(riods.descriptions):
+        # Descriptions for each dataset band
+        attrs['descriptions'] = riods.descriptions
+    if hasattr(riods, 'units') and any(riods.units):
+        # A list of units string for each dataset band
+        attrs['units'] = riods.units
 
     # Parse extra metadata from tags, if supported
     parsers = {'ENVI': _parse_envi}
diff --git a/xarray/coding/variables.py b/xarray/coding/variables.py
index c23e45e4..cc173f78 100644
--- a/xarray/coding/variables.py
+++ b/xarray/coding/variables.py
@@ -33,12 +33,12 @@ class VariableCoder:
     variables in the underlying store.
     """
 
-    def encode(self, variable, name=None):
+    def encode(self, variable, name=None):  # pragma: no cover
         # type: (Variable, Any) -> Variable
         """Convert an encoded variable to a decoded variable."""
         raise NotImplementedError
 
-    def decode(self, variable, name=None):
+    def decode(self, variable, name=None):  # pragma: no cover
         # type: (Variable, Any) -> Variable
         """Convert an decoded variable to a encoded variable."""
         raise NotImplementedError
diff --git a/xarray/core/computation.py b/xarray/core/computation.py
index 9602dcb3..7a898a3b 100644
--- a/xarray/core/computation.py
+++ b/xarray/core/computation.py
@@ -553,7 +553,7 @@ def apply_variable_ufunc(
             numpy_func = func
 
             def func(*arrays):
-                return _apply_with_dask_atop(
+                return _apply_blockwise(
                     numpy_func, arrays, input_dims, output_dims,
                     signature, output_dtypes, output_sizes)
         elif dask == 'allowed':
@@ -601,9 +601,10 @@ def apply_variable_ufunc(
         return tuple(output)
 
 
-def _apply_with_dask_atop(func, args, input_dims, output_dims, signature,
-                          output_dtypes, output_sizes=None):
+def _apply_blockwise(func, args, input_dims, output_dims, signature,
+                     output_dtypes, output_sizes=None):
     import dask.array as da
+    from .dask_array_compat import blockwise
 
     if signature.num_outputs > 1:
         raise NotImplementedError('multiple outputs from apply_ufunc not yet '
@@ -648,15 +649,15 @@ def _apply_with_dask_atop(func, args, input_dims, output_dims, signature,
 
     (out_ind,) = output_dims
 
-    atop_args = []
+    blockwise_args = []
     for arg, dims in zip(args, input_dims):
         # skip leading dimensions that are implicitly added by broadcasting
         ndim = getattr(arg, 'ndim', 0)
         trimmed_dims = dims[-ndim:] if ndim else ()
-        atop_args.extend([arg, trimmed_dims])
+        blockwise_args.extend([arg, trimmed_dims])
 
-    return da.atop(func, out_ind, *atop_args, dtype=dtype, concatenate=True,
-                   new_axes=output_sizes)
+    return blockwise(func, out_ind, *blockwise_args, dtype=dtype,
+                     concatenate=True, new_axes=output_sizes)
 
 
 def apply_array_ufunc(func, *args, dask='forbidden'):
diff --git a/xarray/core/concat.py b/xarray/core/concat.py
index 5698596d..cd59d878 100644
--- a/xarray/core/concat.py
+++ b/xarray/core/concat.py
@@ -94,7 +94,7 @@ def concat(objs, dim=None, data_vars='all', coords='different',
                       FutureWarning, stacklevel=2)
         dim = 'concat_dims'
 
-    if indexers is not None:  # pragma: nocover
+    if indexers is not None:  # pragma: no cover
         warnings.warn('indexers has been renamed to positions; the alias '
                       'will be removed in a future version of xarray',
                       FutureWarning, stacklevel=2)
diff --git a/xarray/core/dask_array_compat.py b/xarray/core/dask_array_compat.py
index 1f2209e8..83822213 100644
--- a/xarray/core/dask_array_compat.py
+++ b/xarray/core/dask_array_compat.py
@@ -4,6 +4,13 @@ import dask.array as da
 import numpy as np
 from dask import __version__ as dask_version
 
+
+try:
+    blockwise = da.blockwise
+except AttributeError:
+    blockwise = da.atop
+
+
 try:
     from dask.array import isin
 except ImportError:  # pragma: no cover
@@ -20,13 +27,13 @@ except ImportError:  # pragma: no cover
         test_elements = da.asarray(test_elements)
         element_axes = tuple(range(element.ndim))
         test_axes = tuple(i + element.ndim for i in range(test_elements.ndim))
-        mapped = da.atop(_isin_kernel, element_axes + test_axes,
-                         element, element_axes,
-                         test_elements, test_axes,
-                         adjust_chunks={axis: lambda _: 1
-                                        for axis in test_axes},
-                         dtype=bool,
-                         assume_unique=assume_unique)
+        mapped = blockwise(_isin_kernel, element_axes + test_axes,
+                           element, element_axes,
+                           test_elements, test_axes,
+                           adjust_chunks={axis: lambda _: 1
+                                          for axis in test_axes},
+                           dtype=bool,
+                           assume_unique=assume_unique)
         result = mapped.any(axis=test_axes)
         if invert:
             result = ~result
diff --git a/xarray/core/utils.py b/xarray/core/utils.py
index 21e0e121..db2fa816 100644
--- a/xarray/core/utils.py
+++ b/xarray/core/utils.py
@@ -319,18 +319,7 @@ def ordered_dict_intersection(first_dict: Mapping[K, V],
     return new_dict
 
 
-class SingleSlotPickleMixin:
-    """Mixin class to add the ability to pickle objects whose state is defined
-    by a single __slots__ attribute. Only necessary under Python 2.
-    """
-    def __getstate__(self):
-        return getattr(self, self.__slots__[0])
-
-    def __setstate__(self, state):
-        setattr(self, self.__slots__[0], state)
-
-
-class Frozen(Mapping[K, V], SingleSlotPickleMixin):
+class Frozen(Mapping[K, V]):
     """Wrapper around an object implementing the mapping interface to make it
     immutable. If you really want to modify the mapping, the mutable version is
     saved under the `mapping` attribute.
@@ -360,7 +349,7 @@ def FrozenOrderedDict(*args, **kwargs) -> Frozen:
     return Frozen(OrderedDict(*args, **kwargs))
 
 
-class SortedKeysDict(MutableMapping[K, V], SingleSlotPickleMixin):
+class SortedKeysDict(MutableMapping[K, V]):
     """An wrapper for dictionary-like objects that always iterates over its
     items in sorted order by key but is otherwise equivalent to the underlying
     mapping.
diff --git a/xarray/plot/facetgrid.py b/xarray/plot/facetgrid.py
index 6ea867ab..2e351f7c 100644
--- a/xarray/plot/facetgrid.py
+++ b/xarray/plot/facetgrid.py
@@ -445,7 +445,7 @@ class FacetGrid:
             ax.yaxis.set_major_locator(y_major_locator)
             for tick in itertools.chain(ax.xaxis.get_major_ticks(),
                                         ax.yaxis.get_major_ticks()):
-                tick.label.set_fontsize(fontsize)
+                tick.label1.set_fontsize(fontsize)
 
         return self
 
diff --git a/xarray/util/print_versions.py b/xarray/util/print_versions.py
index c34faa74..8295259a 100755
--- a/xarray/util/print_versions.py
+++ b/xarray/util/print_versions.py
@@ -1,8 +1,4 @@
-'''utility functions for printing version information
-
-
-see pandas/pandas/util/_print_versions.py'''
-
+"""Utility functions for printing version information."""
 import codecs
 import importlib
 import locale
@@ -78,7 +74,7 @@ def netcdf_and_hdf5_versions():
     return [('libhdf5', libhdf5_version), ('libnetcdf', libnetcdf_version)]
 
 
-def show_versions(as_json=False):
+def show_versions(file=sys.stdout):
     sys_info = get_sys_info()
 
     sys_info.extend(netcdf_and_hdf5_versions())
@@ -98,7 +94,7 @@ def show_versions(as_json=False):
         ("zarr", lambda mod: mod.__version__),
         ("cftime", lambda mod: mod.__version__),
         ("nc_time_axis", lambda mod: mod.__version__),
-        ("PseudonetCDF", lambda mod: mod.__version__),
+        ("PseudoNetCDF", lambda mod: mod.__version__),
         ("rasterio", lambda mod: mod.__version__),
         ("cfgrib", lambda mod: mod.__version__),
         ("iris", lambda mod: mod.__version__),
@@ -135,49 +131,16 @@ def show_versions(as_json=False):
             except Exception:
                 deps_blob.append((modname, 'installed'))
 
-    if (as_json):
-        try:
-            import json
-        except Exception:
-            import simplejson as json
-
-        j = dict(system=dict(sys_info), dependencies=dict(deps_blob))
-
-        if as_json is True:
-            print(j)
-        else:
-            with codecs.open(as_json, "wb", encoding='utf8') as f:
-                json.dump(j, f, indent=2)
-
-    else:
-
-        print("\nINSTALLED VERSIONS")
-        print("------------------")
-
-        for k, stat in sys_info:
-            print("%s: %s" % (k, stat))
-
-        print("")
-        for k, stat in deps_blob:
-            print("%s: %s" % (k, stat))
-
-
-def main():
-    from optparse import OptionParser
-    parser = OptionParser()
-    parser.add_option("-j", "--json", metavar="FILE", nargs=1,
-                      help="Save output as JSON into file, pass in "
-                      "'-' to output to stdout")
-
-    (options, args) = parser.parse_args()
-
-    if options.json == "-":
-        options.json = True
+    print("\nINSTALLED VERSIONS", file=file)
+    print("------------------", file=file)
 
-    show_versions(as_json=options.json)
+    for k, stat in sys_info:
+        print("%s: %s" % (k, stat), file=file)
 
-    return 0
+    print("", file=file)
+    for k, stat in deps_blob:
+        print("%s: %s" % (k, stat), file=file)
 
 
-if __name__ == "__main__":
-    sys.exit(main())
+if __name__ == '__main__':
+    show_versions()


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index ad8df191..73796209 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -1426,6 +1426,7 @@ class ZarrBase(CFEncodedBase):
                 # chunk size should be the same as original
                 assert v.chunks == original[k].chunks
 
+    @pytest.mark.filterwarnings('ignore:Specified Dask chunks')
     def test_manual_chunk(self):
         original = create_test_data().chunk({'dim1': 3, 'dim2': 4, 'dim3': 3})
 
@@ -3285,7 +3286,8 @@ def create_tmp_geotiff(nx=4, ny=3, nz=3,
                        transform_args=[5000, 80000, 1000, 2000.],
                        crs={'units': 'm', 'no_defs': True, 'ellps': 'WGS84',
                             'proj': 'utm', 'zone': 18},
-                       open_kwargs=None):
+                       open_kwargs=None,
+                       additional_attrs=None):
     # yields a temporary geotiff file and a corresponding expected DataArray
     import rasterio
     from rasterio.transform import from_origin
@@ -3308,6 +3310,11 @@ def create_tmp_geotiff(nx=4, ny=3, nz=3,
             *data_shape)
         if transform is None:
             transform = from_origin(*transform_args)
+        if additional_attrs is None:
+            additional_attrs = {
+                'descriptions': tuple('d{}'.format(n + 1) for n in range(nz)),
+                'units': tuple('u{}'.format(n + 1) for n in range(nz)),
+            }
         with rasterio.open(
                 tmp_file, 'w',
                 driver='GTiff', height=ny, width=nx, count=nz,
@@ -3315,6 +3322,8 @@ def create_tmp_geotiff(nx=4, ny=3, nz=3,
                 transform=transform,
                 dtype=rasterio.float32,
                 **open_kwargs) as s:
+            for attr, val in additional_attrs.items():
+                setattr(s, attr, val)
             s.write(data, **write_kwargs)
             dx, dy = s.res[0], -s.res[1]
 
@@ -3334,7 +3343,7 @@ class TestRasterio:
 
     @requires_scipy_or_netCDF4
     def test_serialization(self):
-        with create_tmp_geotiff() as (tmp_file, expected):
+        with create_tmp_geotiff(additional_attrs={}) as (tmp_file, expected):
             # Write it to a netcdf and read again (roundtrip)
             with xr.open_rasterio(tmp_file) as rioda:
                 with create_tmp_file(suffix='.nc') as tmp_nc_file:
@@ -3346,6 +3355,10 @@ class TestRasterio:
         with create_tmp_geotiff() as (tmp_file, expected):
             with xr.open_rasterio(tmp_file) as rioda:
                 assert_allclose(rioda, expected)
+                assert rioda.attrs['scales'] == (1.0, 1.0, 1.0)
+                assert rioda.attrs['offsets'] == (0.0, 0.0, 0.0)
+                assert rioda.attrs['descriptions'] == ('d1', 'd2', 'd3')
+                assert rioda.attrs['units'] == ('u1', 'u2', 'u3')
                 assert isinstance(rioda.attrs['crs'], str)
                 assert isinstance(rioda.attrs['res'], tuple)
                 assert isinstance(rioda.attrs['is_tiled'], np.uint8)
@@ -3369,6 +3382,10 @@ class TestRasterio:
                 assert 'x' not in rioda.coords
                 assert 'y' not in rioda.coords
                 assert 'crs' not in rioda.attrs
+                assert rioda.attrs['scales'] == (1.0, 1.0, 1.0)
+                assert rioda.attrs['offsets'] == (0.0, 0.0, 0.0)
+                assert rioda.attrs['descriptions'] == ('d1', 'd2', 'd3')
+                assert rioda.attrs['units'] == ('u1', 'u2', 'u3')
                 assert isinstance(rioda.attrs['res'], tuple)
                 assert isinstance(rioda.attrs['is_tiled'], np.uint8)
                 assert isinstance(rioda.attrs['transform'], tuple)
@@ -3389,6 +3406,10 @@ class TestRasterio:
                 as (tmp_file, expected):
             with xr.open_rasterio(tmp_file) as rioda:
                 assert_allclose(rioda, expected)
+                assert rioda.attrs['scales'] == (1.0,)
+                assert rioda.attrs['offsets'] == (0.0,)
+                assert isinstance(rioda.attrs['descriptions'], tuple)
+                assert isinstance(rioda.attrs['units'], tuple)
                 assert isinstance(rioda.attrs['crs'], str)
                 assert isinstance(rioda.attrs['res'], tuple)
                 assert isinstance(rioda.attrs['is_tiled'], np.uint8)
@@ -3417,6 +3438,8 @@ class TestRasterio:
                         tmp_file, 'w',
                         driver='GTiff', height=ny, width=nx, count=nz,
                         dtype=rasterio.float32) as s:
+                    s.descriptions = ('nx', 'ny', 'nz')
+                    s.units = ('cm', 'm', 'km')
                     s.write(data)
 
                 # Tests
@@ -3428,6 +3451,10 @@ class TestRasterio:
                                              })
                 with xr.open_rasterio(tmp_file) as rioda:
                     assert_allclose(rioda, expected)
+                    assert rioda.attrs['scales'] == (1.0, 1.0, 1.0)
+                    assert rioda.attrs['offsets'] == (0.0, 0.0, 0.0)
+                    assert rioda.attrs['descriptions'] == ('nx', 'ny', 'nz')
+                    assert rioda.attrs['units'] == ('cm', 'm', 'km')
                     assert isinstance(rioda.attrs['res'], tuple)
                     assert isinstance(rioda.attrs['is_tiled'], np.uint8)
                     assert isinstance(rioda.attrs['transform'], tuple)
diff --git a/xarray/tests/test_backends_common.py b/xarray/tests/test_backends_common.py
new file mode 100644
index 00000000..a6bf118c
--- /dev/null
+++ b/xarray/tests/test_backends_common.py
@@ -0,0 +1,32 @@
+import pytest
+
+from xarray.backends.common import robust_getitem
+
+
+class DummyFailure(Exception):
+    pass
+
+
+class DummyArray:
+    def __init__(self, failures):
+        self.failures = failures
+
+    def __getitem__(self, key):
+        if self.failures:
+            self.failures -= 1
+            raise DummyFailure
+        return 'success'
+
+
+def test_robust_getitem():
+    array = DummyArray(failures=2)
+    with pytest.raises(DummyFailure):
+        array[...]
+    result = robust_getitem(
+        array, ..., catch=DummyFailure, initial_delay=1)
+    assert result == 'success'
+
+    array = DummyArray(failures=3)
+    with pytest.raises(DummyFailure):
+        robust_getitem(array, ..., catch=DummyFailure, initial_delay=1,
+                       max_retries=2)
diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py
index e3f49fd8..f2f4be2e 100644
--- a/xarray/tests/test_computation.py
+++ b/xarray/tests/test_computation.py
@@ -684,7 +684,10 @@ def test_apply_dask_parallelized_errors():
                     output_core_dims=[('y',)])
 
 
+# it's currently impossible to silence these warnings from inside dask.array:
+# https://github.com/dask/dask/issues/3245
 @requires_dask
+@pytest.mark.filterwarnings('ignore:Mean of empty slice')
 def test_apply_dask_multiple_inputs():
     import dask.array as da
 
diff --git a/xarray/tests/test_print_versions.py b/xarray/tests/test_print_versions.py
new file mode 100644
index 00000000..b1755ac2
--- /dev/null
+++ b/xarray/tests/test_print_versions.py
@@ -0,0 +1,9 @@
+import io
+
+import xarray
+
+
+def test_show_versions():
+    f = io.StringIO()
+    xarray.show_versions(file=f)
+    assert 'INSTALLED VERSIONS' in f.getvalue()