From c0388976b3e8d99c6a8fe36e10303656675de95e Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sun, 12 Nov 2017 22:49:39 +0100 Subject: [PATCH 01/13] Fix y-coordinates for non-georeferenced data --- xarray/backends/rasterio_.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index b39299f3541..53a1154a383 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -9,6 +9,7 @@ from dask.utils import SerializableLock as Lock except ImportError: from threading import Lock +from affine import Affine RASTERIO_LOCK = Lock() @@ -132,14 +133,12 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): coords['band'] = np.asarray(riods.indexes) # Get geo coords + transfo = riods.transform * Affine.translation(0.5, 0.5) nx, ny = riods.width, riods.height - 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 + 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) + x, _ = (np.arange(nx), np.zeros(nx)) * transfo + _, y = (np.zeros(ny), np.arange(ny)) * transfo + coords['y'] = y + coords['x'] = x # Attributes attrs = {} From 10195d1a74fa31b52b79237e9a908cdd33ff4f47 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sun, 12 Nov 2017 23:11:49 +0100 Subject: [PATCH 02/13] Add regression test, fix imports --- xarray/backends/rasterio_.py | 2 +- xarray/tests/test_backends.py | 37 +++++++++++++++++++++++++++++++++++ 2 files changed, 38 insertions(+), 1 deletion(-) diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index 53a1154a383..4f66fced3b8 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -9,7 +9,6 @@ from dask.utils import SerializableLock as Lock except ImportError: from threading import Lock -from affine import Affine RASTERIO_LOCK = Lock() @@ -120,6 +119,7 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): """ import rasterio + from rasterio.transform import Affine riods = rasterio.open(filename, mode='r') if cache is None: diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py index 316c97eb487..c1fdaef32a8 100644 --- a/xarray/tests/test_backends.py +++ b/xarray/tests/test_backends.py @@ -1923,6 +1923,43 @@ def test_platecarree(self): assert 'transform' in rioda.attrs assert isinstance(rioda.attrs['transform'], tuple) + def test_notransform(self): + # regression test for https://github.com/pydata/xarray/issues/1686 + import rasterio + import warnings + from rasterio.errors import NotGeoreferencedWarning + + # Create a geotiff file + with warnings.catch_warnings(): + # This warning is expected + warnings.filterwarnings('ignore', category=NotGeoreferencedWarning) + with create_tmp_file(suffix='.tif') as tmp_file: + # data + nx, ny, nz = 4, 3, 3 + data = np.arange(nx*ny*nz, + dtype=rasterio.float32).reshape(nz, ny, nx) + with rasterio.open( + tmp_file, 'w', + driver='GTiff', height=ny, width=nx, count=nz, + dtype=rasterio.float32) as s: + s.write(data) + + # Tests + expected = DataArray(data, + dims=('band', 'y', 'x'), + coords={'band': [1, 2, 3], + 'y': [0.5, 1.5, 2.5], + 'x': [0.5, 1.5, 2.5, 3.5], + }) + with xr.open_rasterio(tmp_file) as rioda: + assert_allclose(rioda, expected) + assert 'res' in rioda.attrs + assert isinstance(rioda.attrs['res'], tuple) + assert 'is_tiled' in rioda.attrs + assert isinstance(rioda.attrs['is_tiled'], np.uint8) + assert 'transform' in rioda.attrs + assert isinstance(rioda.attrs['transform'], tuple) + def test_indexing(self): import rasterio From 67a88aec3695c862c96b2929ade31f33cc58cf14 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sun, 12 Nov 2017 23:40:57 +0100 Subject: [PATCH 03/13] version conundrum --- xarray/backends/rasterio_.py | 7 ++++++- xarray/tests/test_backends.py | 3 +-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index 4f66fced3b8..5915af3193d 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -1,5 +1,6 @@ import os from collections import OrderedDict +from distutils.version import LooseVersion import numpy as np from .. import DataArray @@ -133,7 +134,11 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): coords['band'] = np.asarray(riods.indexes) # Get geo coords - transfo = riods.transform * Affine.translation(0.5, 0.5) + if LooseVersion(rasterio.__version__) < '1.0': + transfo = riods.affine * Affine.translation(0.5, 0.5) + else: + transfo = riods.transform * Affine.translation(0.5, 0.5) + nx, ny = riods.width, riods.height x, _ = (np.arange(nx), np.zeros(nx)) * transfo _, y = (np.zeros(ny), np.arange(ny)) * transfo diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py index c1fdaef32a8..6827100b4b8 100644 --- a/xarray/tests/test_backends.py +++ b/xarray/tests/test_backends.py @@ -1927,12 +1927,11 @@ def test_notransform(self): # regression test for https://github.com/pydata/xarray/issues/1686 import rasterio import warnings - from rasterio.errors import NotGeoreferencedWarning # Create a geotiff file with warnings.catch_warnings(): # This warning is expected - warnings.filterwarnings('ignore', category=NotGeoreferencedWarning) + warnings.filterwarnings('ignore', category=UserWarning) with create_tmp_file(suffix='.tif') as tmp_file: # data nx, ny, nz = 4, 3, 3 From 9595a7f6a8a1d2234e1266491fa67b4651dd47a6 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Mon, 13 Nov 2017 00:26:45 +0100 Subject: [PATCH 04/13] What's new --- doc/whats-new.rst | 3 +++ xarray/backends/rasterio_.py | 13 +++++++------ 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/doc/whats-new.rst b/doc/whats-new.rst index ffcc4eece51..c57a13823da 100644 --- a/doc/whats-new.rst +++ b/doc/whats-new.rst @@ -58,6 +58,9 @@ Bug fixes (:issue:`1697`). By `Stephan Hoyer `_. +- Fixed wrong y-coordinates for some rasterio files (:issue:`1686`). + By `Fabien Maussion `_. + Testing ~~~~~~~ diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index 5915af3193d..bce9df13f17 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -133,15 +133,16 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): raise ValueError('Unknown dims') coords['band'] = np.asarray(riods.indexes) - # Get geo coords + # Get coordinates if LooseVersion(rasterio.__version__) < '1.0': - transfo = riods.affine * Affine.translation(0.5, 0.5) + transform = riods.affine else: - transfo = riods.transform * Affine.translation(0.5, 0.5) - + transform = riods.transform + # Xarray's convention is pixel centered + transform = transform * Affine.translation(0.5, 0.5) nx, ny = riods.width, riods.height - x, _ = (np.arange(nx), np.zeros(nx)) * transfo - _, y = (np.zeros(ny), np.arange(ny)) * transfo + x, _ = (np.arange(nx), np.zeros(nx)) * transform + _, y = (np.zeros(ny), np.arange(ny)) * transform coords['y'] = y coords['x'] = x From b6eddc05f67baf0c175e4f6837aa94577363a0bf Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Wed, 22 Nov 2017 09:40:32 +0100 Subject: [PATCH 05/13] Add recti block --- xarray/backends/rasterio_.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index bce9df13f17..9971fb7cd77 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -120,7 +120,6 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): """ import rasterio - from rasterio.transform import Affine riods = rasterio.open(filename, mode='r') if cache is None: @@ -138,11 +137,15 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): transform = riods.affine else: transform = riods.transform - # Xarray's convention is pixel centered - transform = transform * Affine.translation(0.5, 0.5) nx, ny = riods.width, riods.height - x, _ = (np.arange(nx), np.zeros(nx)) * transform - _, y = (np.zeros(ny), np.arange(ny)) * transform + if transform.is_rectilinear: + # Xarray's convention is pixel centered + x, _ = (np.arange(nx)+0.5, np.zeros(nx)+0.5) * transform + _, y = (np.zeros(ny)+0.5, np.arange(ny)+0.5) * transform + else: + # Can this happen? + raise RuntimeError() + coords['y'] = y coords['x'] = x From 2f6c28606cd1a0dfb88fbc06b05ad46273bac429 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 2 Dec 2017 21:56:06 +0100 Subject: [PATCH 06/13] Add 2d coords --- xarray/backends/rasterio_.py | 50 ++++++++++++++++++++--------------- xarray/tests/test_backends.py | 33 ++++++++++++++++++++++- 2 files changed, 61 insertions(+), 22 deletions(-) diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index 9971fb7cd77..95951835c9d 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -84,7 +84,8 @@ def __getitem__(self, key): return out -def open_rasterio(filename, chunks=None, cache=None, lock=None): +def open_rasterio(filename, parse_coordinates=True, chunks=None, cache=None, + lock=None): """Open a file with rasterio (experimental). This should work with any file that rasterio can open (most often: @@ -98,11 +99,10 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): ---------- filename : str Path to the file to open. - - Returns - ------- - data : DataArray - The newly created DataArray. + parse_coordinates : bool, optional + Whether to parse the x and y coordinates out of the file's + ``transform`` attribute or not. It can be useful to switch this off + if your files are large or if you don't need the coordinates. chunks : int, tuple or dict, optional Chunk sizes along each dimension, e.g., ``5``, ``(5, 5)`` or ``{'x': 5, 'y': 5}``. If chunks is provided, it used to load the new @@ -117,6 +117,11 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): :py:func:`dask.array.from_array`. By default, a global lock is used to avoid issues with concurrent access to the same file when using dask's multithreaded backend. + + Returns + ------- + data : DataArray + The newly created DataArray. """ import rasterio @@ -133,21 +138,24 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): coords['band'] = np.asarray(riods.indexes) # Get coordinates - if LooseVersion(rasterio.__version__) < '1.0': - transform = riods.affine - else: - transform = riods.transform - nx, ny = riods.width, riods.height - if transform.is_rectilinear: - # Xarray's convention is pixel centered - x, _ = (np.arange(nx)+0.5, np.zeros(nx)+0.5) * transform - _, y = (np.zeros(ny)+0.5, np.arange(ny)+0.5) * transform - else: - # Can this happen? - raise RuntimeError() - - coords['y'] = y - coords['x'] = x + if parse_coordinates: + if LooseVersion(rasterio.__version__) < '1.0': + transform = riods.affine + else: + transform = riods.transform + nx, ny = riods.width, riods.height + if transform.is_rectilinear: + # for 1d coordinates we need two calls: columns and rows + # xarray's convention is pixel centered + x, _ = (np.arange(nx)+0.5, np.zeros(nx)+0.5) * transform + _, y = (np.zeros(ny)+0.5, np.arange(ny)+0.5) * transform + coords['y'] = y + coords['x'] = x + else: + # 2d coordinates + x, y = np.meshgrid(np.arange(nx), np.arange(ny)) * transform + coords['yy'] = (('y', 'x'), y) + coords['xx'] = (('y', 'x'), x) # Attributes attrs = {} diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py index 6827100b4b8..cad4eb72a3d 100644 --- a/xarray/tests/test_backends.py +++ b/xarray/tests/test_backends.py @@ -1819,7 +1819,7 @@ class TestPyNioAutocloseTrue(TestPyNio): class TestRasterio(TestCase): @requires_scipy_or_netCDF4 - def test_serialization(self): + def test_netcdf_serialization(self): import rasterio from rasterio.transform import from_origin @@ -1885,6 +1885,37 @@ def test_utm(self): assert 'transform' in rioda.attrs assert isinstance(rioda.attrs['transform'], tuple) + def test_non_rectilinear(self): + import rasterio + from rasterio.transform import from_origin + + # Create a geotiff file with 2d coordinates + with create_tmp_file(suffix='.tif') as tmp_file: + # data + nx, ny, nz = 4, 3, 3 + data = np.arange(nx*ny*nz, + dtype=rasterio.float32).reshape(nz, ny, nx) + transform = from_origin(0, 3, 1, 1).rotation(45) + with rasterio.open( + tmp_file, 'w', + driver='GTiff', height=ny, width=nx, count=nz, + transform=transform, + dtype=rasterio.float32) as s: + s.write(data) + + # Not much we can test here. See if things are parsed ok + with xr.open_rasterio(tmp_file) as rioda: + assert rioda['xx'].shape == rioda['yy'].shape + assert 'x' not in rioda.coords + assert 'y' not in rioda.coords + assert 'crs' not in rioda.attrs + assert 'res' in rioda.attrs + assert isinstance(rioda.attrs['res'], tuple) + assert 'is_tiled' in rioda.attrs + assert isinstance(rioda.attrs['is_tiled'], np.uint8) + assert 'transform' in rioda.attrs + assert isinstance(rioda.attrs['transform'], tuple) + def test_platecarree(self): import rasterio From 7ab95c8983305c4a69541c83bfbe57d95a898461 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 2 Dec 2017 21:59:40 +0100 Subject: [PATCH 07/13] whatsnew mixup --- doc/whats-new.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/whats-new.rst b/doc/whats-new.rst index e799844859c..899175af45f 100644 --- a/doc/whats-new.rst +++ b/doc/whats-new.rst @@ -28,6 +28,8 @@ Enhancements Bug fixes ~~~~~~~~~ +.. _whats-new.0.10.0: + v0.10.0 (20 November 2017) -------------------------- From 95016ce5b6390004d565c973c6403867ce29368f Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 2 Dec 2017 22:33:47 +0100 Subject: [PATCH 08/13] center pixels --- xarray/backends/rasterio_.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index 8d091f9c763..8643bd6c620 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -151,7 +151,8 @@ def open_rasterio(filename, parse_coordinates=True, chunks=None, cache=None, coords['x'] = x else: # 2d coordinates - x, y = np.meshgrid(np.arange(nx), np.arange(ny)) * transform + x, y = (np.meshgrid(np.arange(nx)+0.5, np.arange(ny)+0.5) * + transform) coords['yy'] = (('y', 'x'), y) coords['xx'] = (('y', 'x'), x) From 2b4fc4150944dc2fa6a4d3b096e9fedbb609dcbd Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 14 Dec 2017 11:56:17 +0100 Subject: [PATCH 09/13] More sensible defaults --- xarray/backends/rasterio_.py | 43 ++++++++++++++++++++--------------- xarray/tests/test_backends.py | 18 +++++++++++---- 2 files changed, 38 insertions(+), 23 deletions(-) diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index 8643bd6c620..639f6627547 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -1,4 +1,5 @@ import os +import warnings from collections import OrderedDict from distutils.version import LooseVersion import numpy as np @@ -82,7 +83,7 @@ def __getitem__(self, key): return out -def open_rasterio(filename, parse_coordinates=True, chunks=None, cache=None, +def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, lock=None): """Open a file with rasterio (experimental). @@ -99,8 +100,10 @@ def open_rasterio(filename, parse_coordinates=True, chunks=None, cache=None, Path to the file to open. parse_coordinates : bool, optional Whether to parse the x and y coordinates out of the file's - ``transform`` attribute or not. It can be useful to switch this off - if your files are large or if you don't need the coordinates. + ``transform`` attribute or not. The default is to automatically + parse the coordinates if they are rectilinear (1D) and don't parse them + if they aren't (2D). It can be useful to set `parse_coordinates=False` + if your files are very large or if you don't need the coordinates. chunks : int, tuple or dict, optional Chunk sizes along each dimension, e.g., ``5``, ``(5, 5)`` or ``{'x': 5, 'y': 5}``. If chunks is provided, it used to load the new @@ -136,25 +139,29 @@ def open_rasterio(filename, parse_coordinates=True, chunks=None, cache=None, coords['band'] = np.asarray(riods.indexes) # Get coordinates - if parse_coordinates: - if LooseVersion(rasterio.__version__) < '1.0': - transform = riods.affine - else: - transform = riods.transform - nx, ny = riods.width, riods.height - if transform.is_rectilinear: - # for 1d coordinates we need two calls: columns and rows - # xarray's convention is pixel centered + if LooseVersion(rasterio.__version__) < '1.0': + transform = riods.affine + else: + transform = riods.transform + if transform.is_rectilinear: + # 1d coordinates + parse = True if parse_coordinates is None else parse_coordinates + if parse: + nx, ny = riods.width, riods.height + # xarray coordinates are pixel centered x, _ = (np.arange(nx)+0.5, np.zeros(nx)+0.5) * transform _, y = (np.zeros(ny)+0.5, np.arange(ny)+0.5) * transform coords['y'] = y coords['x'] = x - else: - # 2d coordinates - x, y = (np.meshgrid(np.arange(nx)+0.5, np.arange(ny)+0.5) * - transform) - coords['yy'] = (('y', 'x'), y) - coords['xx'] = (('y', 'x'), x) + else: + # 2d coordinates + parse = False if parse_coordinates is None else parse_coordinates + if parse: + warnings.warn("The file coordinates' transformation isn't " + "rectilinear: xarray won't parse the coordinates " + "in this case. Set `parse_coordinates=False` to " + "suppress this warning.", + RuntimeWarning, stacklevel=3) # Attributes attrs = {} diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py index 84007fac836..87ed4e14f0b 100644 --- a/xarray/tests/test_backends.py +++ b/xarray/tests/test_backends.py @@ -1915,6 +1915,11 @@ def test_utm(self): assert 'transform' in rioda.attrs assert isinstance(rioda.attrs['transform'], tuple) + # Check no parse coords + with xr.open_rasterio(tmp_file, parse_coordinates=False) as rioda: + assert 'x' not in rioda.coords + assert 'y' not in rioda.coords + def test_non_rectilinear(self): import rasterio from rasterio.transform import from_origin @@ -1933,19 +1938,22 @@ def test_non_rectilinear(self): dtype=rasterio.float32) as s: s.write(data) - # Not much we can test here. See if things are parsed ok + # Default is to not parse coords with xr.open_rasterio(tmp_file) as rioda: - assert rioda['xx'].shape == rioda['yy'].shape assert 'x' not in rioda.coords assert 'y' not in rioda.coords assert 'crs' not in rioda.attrs - assert 'res' in rioda.attrs assert isinstance(rioda.attrs['res'], tuple) - assert 'is_tiled' in rioda.attrs assert isinstance(rioda.attrs['is_tiled'], np.uint8) - assert 'transform' in rioda.attrs assert isinstance(rioda.attrs['transform'], tuple) + # See if a warning is raised if we force it + with self.assertWarns("transformation isn't rectilinear"): + with xr.open_rasterio(tmp_file, + parse_coordinates=True) as rioda: + assert 'x' not in rioda.coords + assert 'y' not in rioda.coords + def test_platecarree(self): import rasterio From a3663048b0ff834cbf24d8fa75f91fd676fc172a Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 14 Dec 2017 16:22:59 +0100 Subject: [PATCH 10/13] Some doc --- xarray/backends/rasterio_.py | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index 8417429b9ef..b33d2cdfbba 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -125,6 +125,14 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, `_ for more information). + You can generate coordinates from a non-rectilinear transformation with:: + + from affine import Affine + transform = Affine(*da.attrs['transform']) + nx, ny = , da.sizes['y'] + x, y = np.meshgrid(np.arange(nx)+0.5, np.arange(ny)+0.5) * transform + + Parameters ---------- filename : str @@ -132,8 +140,8 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, parse_coordinates : bool, optional Whether to parse the x and y coordinates out of the file's ``transform`` attribute or not. The default is to automatically - parse the coordinates if they are rectilinear (1D) and don't parse them - if they aren't (2D). It can be useful to set `parse_coordinates=False` + parse the coordinates only if they are rectilinear (1D). + It can be useful to set `parse_coordinates=False` if your files are very large or if you don't need the coordinates. chunks : int, tuple or dict, optional Chunk sizes along each dimension, e.g., ``5``, ``(5, 5)`` or @@ -195,7 +203,12 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, RuntimeWarning, stacklevel=3) # Attributes - attrs = {} + attrs = dict() + # Affine transformation matrix (always available) + # This describes coefficients mapping pixel coordinates to CRS + # For serialization store as tuple of 6 floats, the last row being + # always (0, 0, 1) per definition (see https://github.com/sgillies/affine) + attrs['transform'] = tuple(transform)[:6] if hasattr(riods, 'crs') and riods.crs: # CRS is a dict-like object specific to rasterio # If CRS is not None, we convert it back to a PROJ4 string using @@ -208,10 +221,6 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, # Is the TIF tiled? (bool) # We cast it to an int for netCDF compatibility attrs['is_tiled'] = np.uint8(riods.is_tiled) - if hasattr(riods, 'transform'): - # Affine transformation matrix (tuple of floats) - # Describes coefficients mapping pixel coordinates to CRS - attrs['transform'] = tuple(riods.transform) # Parse extra metadata from tags, if supported parsers = {'ENVI': _parse_envi} From 4172e6e1bbb0bce996bdbc7885b780ce15186358 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 14 Dec 2017 16:39:00 +0100 Subject: [PATCH 11/13] More docs --- doc/io.rst | 10 +++++++--- xarray/backends/rasterio_.py | 11 ++++++----- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/doc/io.rst b/doc/io.rst index 14e82d4aacc..4f5ecc84020 100644 --- a/doc/io.rst +++ b/doc/io.rst @@ -480,10 +480,14 @@ rasterio is installed. Here is an example of how to use [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 ... + * y (y) float64 2.827e+06 2.826e+06 2.826e+06 2.826e+06 2.826e+06 ... + * x (x) float64 1.021e+05 1.024e+05 1.027e+05 1.03e+05 1.033e+05 ... Attributes: - crs: +init=epsg:32618 + res: (300.0379266750948, 300.041782729805) + transform: (300.0379266750948, 0.0, 101985.0, 0.0, -300.041782729805, 28... + is_tiled: 0 + 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 diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index b33d2cdfbba..5863aa7ca53 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -125,11 +125,12 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, `_ for more information). - You can generate coordinates from a non-rectilinear transformation with:: + You can generate 2D coordinates from the file's attributes with:: from affine import Affine + da = xr.open_rasterio('path_to_file.tif') transform = Affine(*da.attrs['transform']) - nx, ny = , da.sizes['y'] + nx, ny = da.sizes['x'], da.sizes['y'] x, y = np.meshgrid(np.arange(nx)+0.5, np.arange(ny)+0.5) * transform @@ -141,7 +142,7 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, Whether to parse the x and y coordinates out of the file's ``transform`` attribute or not. The default is to automatically parse the coordinates only if they are rectilinear (1D). - It can be useful to set `parse_coordinates=False` + It can be useful to set ``parse_coordinates=False`` if your files are very large or if you don't need the coordinates. chunks : int, tuple or dict, optional Chunk sizes along each dimension, e.g., ``5``, ``(5, 5)`` or @@ -184,7 +185,7 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, transform = riods.transform if transform.is_rectilinear: # 1d coordinates - parse = True if parse_coordinates is None else parse_coordinates + parse = True if (parse_coordinates is None) else parse_coordinates if parse: nx, ny = riods.width, riods.height # xarray coordinates are pixel centered @@ -194,7 +195,7 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, coords['x'] = x else: # 2d coordinates - parse = False if parse_coordinates is None else parse_coordinates + parse = False if (parse_coordinates is None) else parse_coordinates if parse: warnings.warn("The file coordinates' transformation isn't " "rectilinear: xarray won't parse the coordinates " From 7e284024329c95a69042ad6af6ec25a2af887742 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 14 Dec 2017 16:48:38 +0100 Subject: [PATCH 12/13] Remove redundant tests --- xarray/tests/test_backends.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py index 4960857d125..63c821e08a5 100644 --- a/xarray/tests/test_backends.py +++ b/xarray/tests/test_backends.py @@ -2151,13 +2151,9 @@ def test_utm(self): }) with xr.open_rasterio(tmp_file) as rioda: assert_allclose(rioda, expected) - assert 'crs' in rioda.attrs assert isinstance(rioda.attrs['crs'], basestring) - assert 'res' in rioda.attrs assert isinstance(rioda.attrs['res'], tuple) - assert 'is_tiled' in rioda.attrs assert isinstance(rioda.attrs['is_tiled'], np.uint8) - assert 'transform' in rioda.attrs assert isinstance(rioda.attrs['transform'], tuple) # Check no parse coords @@ -2228,13 +2224,9 @@ def test_platecarree(self): }) with xr.open_rasterio(tmp_file) as rioda: assert_allclose(rioda, expected) - assert 'crs' in rioda.attrs assert isinstance(rioda.attrs['crs'], basestring) - assert 'res' in rioda.attrs assert isinstance(rioda.attrs['res'], tuple) - assert 'is_tiled' in rioda.attrs assert isinstance(rioda.attrs['is_tiled'], np.uint8) - assert 'transform' in rioda.attrs assert isinstance(rioda.attrs['transform'], tuple) def test_notransform(self): @@ -2266,11 +2258,8 @@ def test_notransform(self): }) with xr.open_rasterio(tmp_file) as rioda: assert_allclose(rioda, expected) - assert 'res' in rioda.attrs assert isinstance(rioda.attrs['res'], tuple) - assert 'is_tiled' in rioda.attrs assert isinstance(rioda.attrs['is_tiled'], np.uint8) - assert 'transform' in rioda.attrs assert isinstance(rioda.attrs['transform'], tuple) def test_indexing(self): From 69566f9809f655ccaa43eac6d8e13bc01487e31c Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 13 Jan 2018 19:49:26 +0100 Subject: [PATCH 13/13] Review --- xarray/backends/rasterio_.py | 2 +- xarray/tests/test_backends.py | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py index 5863aa7ca53..46802b84217 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -185,7 +185,7 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, transform = riods.transform if transform.is_rectilinear: # 1d coordinates - parse = True if (parse_coordinates is None) else parse_coordinates + parse = True if parse_coordinates is None else parse_coordinates if parse: nx, ny = riods.width, riods.height # xarray coordinates are pixel centered diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py index 63c821e08a5..5622eb61d76 100644 --- a/xarray/tests/test_backends.py +++ b/xarray/tests/test_backends.py @@ -2236,8 +2236,10 @@ def test_notransform(self): # Create a geotiff file with warnings.catch_warnings(): - # This warning is expected - warnings.filterwarnings('ignore', category=UserWarning) + # rasterio throws a NotGeoreferencedWarning here, which is + # expected since we test rasterio's defaults in this case. + warnings.filterwarnings('ignore', category=UserWarning, + message='Dataset has no geotransform set') with create_tmp_file(suffix='.tif') as tmp_file: # data nx, ny, nz = 4, 3, 3