diff --git a/doc/io.rst b/doc/io.rst index eac941b336b..c177496f6f2 100644 --- a/doc/io.rst +++ b/doc/io.rst @@ -512,10 +512,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 3acc5551173..c624c1f5ff8 100644 --- a/xarray/backends/rasterio_.py +++ b/xarray/backends/rasterio_.py @@ -1,5 +1,7 @@ import os +import warnings from collections import OrderedDict +from distutils.version import LooseVersion import numpy as np from .. import DataArray @@ -113,7 +115,8 @@ def default(s): return parsed_meta -def open_rasterio(filename, chunks=None, cache=None, lock=None): +def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, + lock=None): """Open a file with rasterio (experimental). This should work with any file that rasterio can open (most often: @@ -123,15 +126,25 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): `_ for more information). + 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['x'], da.sizes['y'] + x, y = np.meshgrid(np.arange(nx)+0.5, np.arange(ny)+0.5) * transform + + Parameters ---------- 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. The default is to automatically + 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 ``{'x': 5, 'y': 5}``. If chunks is provided, it used to load the new @@ -146,6 +159,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 @@ -161,18 +179,38 @@ def open_rasterio(filename, chunks=None, cache=None, lock=None): raise ValueError('Unknown dims') coords['band'] = np.asarray(riods.indexes) - # Get geo coords - 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) + # Get coordinates + 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 + 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 = {} + 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 diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py index 126da0926d4..3ea0fea9a37 100644 --- a/xarray/tests/test_backends.py +++ b/xarray/tests/test_backends.py @@ -2227,30 +2227,96 @@ def test_utm(self): with create_tmp_geotiff() as (tmp_file, expected): 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 + 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 + + # 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) + + # Default is to not parse coords + with xr.open_rasterio(tmp_file) as rioda: + assert 'x' not in rioda.coords + assert 'y' not in rioda.coords + assert 'crs' not in rioda.attrs + assert isinstance(rioda.attrs['res'], tuple) + assert isinstance(rioda.attrs['is_tiled'], np.uint8) + 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): with create_tmp_geotiff(8, 10, 1, transform_args=[1, 2, 0.5, 2.], crs='+proj=latlong') \ as (tmp_file, expected): 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): + # regression test for https://github.com/pydata/xarray/issues/1686 + import rasterio + import warnings + + # Create a geotiff file + with warnings.catch_warnings(): + # 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 + 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 isinstance(rioda.attrs['res'], tuple) + assert isinstance(rioda.attrs['is_tiled'], np.uint8) + assert isinstance(rioda.attrs['transform'], tuple) + def test_indexing(self): with create_tmp_geotiff(8, 10, 3, transform_args=[1, 2, 0.5, 2.], crs='+proj=latlong') as (tmp_file, expected):