Center the coordinates to pixels for rasterio backend (pydata#1468)

* Center the coordinates for rasterio backend

Rasterio uses edge-based coordinates, which contradict the treatment of coordinates in xarray as being centered over the pixels. This centers them with an offset of half the resolution.

* Add documentation to whats-new.rst + docstring

* Minor aesthetic correction

* Remove unnecessary line breaks

Minor cleanup for @fmaussion

doc/whats-new.rst

@@ -28,6 +28,10 @@ Enhancements
 Bug fixes
 ~~~~~~~~~
 
+- :py:func:`~xarray.open_rasterio` method now shifts the rasterio
+  coordinates so that they are centered in each pixel.
+  By `Greg Brener <https://github.com/gbrener>`_.
+
 .. _whats-new.0.9.6:
 
 v0.9.6 (8 June 2017)

xarray/backends/rasterio_.py

@@ -87,7 +87,10 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None):
 
     This should work with any file that rasterio can open (most often:
     geoTIFF). The x and y coordinates are generated automatically from the
-    file's geoinformation.
+    file's geoinformation, shifted to the center of each pixel (see
+    `"PixelIsArea" Raster Space
+    <http://web.archive.org/web/20160326194152/http://remotesensing.org/geotiff/spec/geotiff2.5.html#2.5.2>`_
+    for more information).
 
     Parameters
     ----------
@@ -132,8 +135,10 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None):
     dx, dy = riods.res[0], -riods.res[1]
     x0 = riods.bounds.right if dx < 0 else riods.bounds.left
     y0 = riods.bounds.top if dy < 0 else riods.bounds.bottom
-    coords['y'] = np.linspace(start=y0, num=ny, stop=(y0 + (ny - 1) * dy))
-    coords['x'] = np.linspace(start=x0, num=nx, stop=(x0 + (nx - 1) * dx))
+    coords['y'] = np.linspace(start=y0 + dy/2, num=ny,
+                              stop=(y0 + (ny - 1) * dy) + dy/2)
+    coords['x'] = np.linspace(start=x0 + dx/2, num=nx,
+                              stop=(x0 + (nx - 1) * dx) + dx/2)
 
     # Attributes
     attrs = {}

xarray/tests/test_backends.py

@@ -1443,7 +1443,6 @@ class TestPyNioAutocloseTrue(TestPyNio):
 class TestRasterio(TestCase):
 
     def test_serialization_utm(self):
-
         import rasterio
         from rasterio.transform import from_origin
 
@@ -1462,13 +1461,15 @@ def test_serialization_utm(self):
                     transform=transform,
                     dtype=rasterio.float32) as s:
                 s.write(data)
+                dx, dy = s.res[0], -s.res[1]
 
             # Tests
             expected = DataArray(data, dims=('band', 'y', 'x'),
-                                 coords={'band': [1, 2, 3],
-                                         'y': -np.arange(ny) * 2000 + 80000,
-                                         'x': np.arange(nx) * 1000 + 5000,
-                                         })
+                                 coords={
+                                     'band': [1, 2, 3],
+                                     'y': -np.arange(ny) * 2000 + 80000 + dy/2,
+                                     'x': np.arange(nx) * 1000 + 5000 + dx/2,
+                                 })
             with xr.open_rasterio(tmp_file) as rioda:
                 assert_allclose(rioda, expected)
                 assert 'crs' in rioda.attrs
@@ -1504,13 +1505,14 @@ def test_serialization_platecarree(self):
                     transform=transform,
                     dtype=rasterio.float32) as s:
                 s.write(data, indexes=1)
+                dx, dy = s.res[0], -s.res[1]
 
             # Tests
             expected = DataArray(data[np.newaxis, ...],
                                  dims=('band', 'y', 'x'),
                                  coords={'band': [1],
-                                         'y': -np.arange(ny)*2 + 2,
-                                         'x': np.arange(nx)*0.5 + 1,
+                                         'y': -np.arange(ny)*2 + 2 + dy/2,
+                                         'x': np.arange(nx)*0.5 + 1 + dx/2,
                                          })
             with xr.open_rasterio(tmp_file) as rioda:
                 assert_allclose(rioda, expected)
@@ -1548,11 +1550,12 @@ def test_indexing(self):
                     transform=transform,
                     dtype=rasterio.float32) as s:
                 s.write(data)
+                dx, dy = s.res[0], -s.res[1]
 
             # ref
             expected = DataArray(data, dims=('band', 'y', 'x'),
-                                 coords={'x': np.arange(nx)*0.5 + 1,
-                                         'y': -np.arange(ny)*2 + 2,
+                                 coords={'x': (np.arange(nx)*0.5 + 1) + dx/2,
+                                         'y': (-np.arange(ny)*2 + 2) + dy/2,
                                          'band': [1, 2, 3]})
 
             with xr.open_rasterio(tmp_file, cache=False) as actual:
@@ -1640,11 +1643,12 @@ def test_caching(self):
                     transform=transform,
                     dtype=rasterio.float32) as s:
                 s.write(data)
+                dx, dy = s.res[0], -s.res[1]
 
             # ref
             expected = DataArray(data, dims=('band', 'y', 'x'),
-                                 coords={'x': np.arange(nx)*0.5 + 1,
-                                         'y': -np.arange(ny)*2 + 2,
+                                 coords={'x': (np.arange(nx)*0.5 + 1) + dx/2,
+                                         'y': (-np.arange(ny)*2 + 2) + dy/2,
                                          'band': [1, 2, 3]})
 
             # Cache is the default
@@ -1683,6 +1687,7 @@ def test_chunks(self):
                     transform=transform,
                     dtype=rasterio.float32) as s:
                 s.write(data)
+                dx, dy = s.res[0], -s.res[1]
 
             # Chunk at open time
             with xr.open_rasterio(tmp_file, chunks=(1, 2, 2)) as actual:
@@ -1693,8 +1698,8 @@ def test_chunks(self):
 
                 # ref
                 expected = DataArray(data, dims=('band', 'y', 'x'),
-                                     coords={'x': np.arange(nx)*0.5 + 1,
-                                             'y': -np.arange(ny)*2 + 2,
+                                     coords={'x': np.arange(nx)*0.5 + 1 + dx/2,
+                                             'y': -np.arange(ny)*2 + 2 + dy/2,
                                              'band': [1, 2, 3]})
 
                 # do some arithmetic