Skip to content

Center the coordinates to pixels for rasterio backend #1468

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
shoyer merged 5 commits into from
Jul 5, 2017
Merged

Center the coordinates to pixels for rasterio backend #1468

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1a8c6f1
Center the coordinates for rasterio backend
gbrener Jun 29, 2017
2bc35a8
Add documentation to whats-new.rst + docstring
gbrener Jun 30, 2017
55a2a27
Minor aesthetic correction
gbrener Jun 30, 2017
8691766
Merge branch 'master' into center_coords
gbrener Jul 5, 2017
aeeef43
Remove unnecessary line breaks
gbrener Jul 5, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions doc/whats-new.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,10 @@ v0.9.7 (unreleased)
Enhancements
~~~~~~~~~~~~

- :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>`_.

@fmaussion fmaussion Jul 1, 2017

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This change should rather appear in "Bug fixes" I think


Bug fixes
~~~~~~~~~

Expand Down
13 changes: 10 additions & 3 deletions xarray/backends/rasterio_.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>`_

@fmaussion fmaussion Jul 1, 2017

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for looking into this. Isn't there a better ref than on the web archive?

@gbrener gbrener Jul 5, 2017

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Happy to help @fmaussion . The web archive link is used by the GeoTIFF website (first link here: http://trac.osgeo.org/geotiff). A quick Google search didn't yield anything better.

for more information).

Parameters
----------
Expand Down Expand Up @@ -132,8 +135,12 @@ 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,

@fmaussion fmaussion Jul 1, 2017

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit: unnecessary linebreak

@gbrener gbrener Jul 5, 2017 •
edited
Loading

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was just to satisfy flake8's line-length requirement. The three lines were for consistency, but I changed it (in my recent commit) to be two lines.

stop=(y0 + (ny - 1) * dy) + dy/2)
coords['x'] = np.linspace(start=x0 + dx/2,
num=nx,

@fmaussion fmaussion Jul 1, 2017

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit: unnecessary linebreak

@gbrener gbrener Jul 5, 2017

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was just to satisfy flake8's line-length requirement. The three lines were for consistency, but I changed it (in my recent commit) to be two lines.

stop=(x0 + (nx - 1) * dx) + dx/2)

# Attributes
attrs = {}
Expand Down
31 changes: 18 additions & 13 deletions xarray/tests/test_backends.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1443,7 +1443,6 @@ class TestPyNioAutocloseTrue(TestPyNio):
class TestRasterio(TestCase):

def test_serialization_utm(self):

import rasterio
from rasterio.transform import from_origin

Expand All @@ -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
Expand Down Expand Up @@ -1498,13 +1499,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)
Expand Down Expand Up @@ -1536,11 +1538,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:
Expand Down Expand Up @@ -1628,11 +1631,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
Expand Down Expand Up @@ -1671,6 +1675,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:
Expand All @@ -1681,8 +1686,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
Expand Down