Add RasterIO backend #1260

fmaussion · 2017-02-10T20:59:31Z

Follow-up to #1070

This is my first backend so this is going to be a bit more work than I expected, but with your help we should be able to get through this.

A long todo list:

closes ENH: Optional Read-Only RasterIO backend #790
add tests
passes git diff upstream/master | flake8 --diff
whatsnew entry
wrap getitem for lazy data read
~~make coordinates variables lazy too~~

fmaussion · 2017-02-10T21:00:16Z

Before I'll get more into details with what needs to be done with rasterio itself, I'd like to get some xarray internals ready first. No need to do a full review yet, but I'd appreciate help with the following points:

The order of the variables and coordinates is random. The current __repr__ can look like this:

<xarray.Dataset>
Dimensions:  (band: 1, x: 4, y: 3)
Coordinates:
  * x        (x) float64 1.0 1.5 2.0 2.5
  * band     (band) int64 1
  * y        (y) float64 2.0 1.0 0.0
Data variables:
    lat      (y, x) float64 2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0
    raster   (band, y, x) float32 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 ...
    lon      (y, x) float64 1.0 1.5 2.0 2.5 1.0 1.5 2.0 2.5 1.0 1.5 2.0 2.5
Attributes:
    crs: CRS({'wktext': True, 'proj': 'longlat', 'no_defs': True, 'ellps': 'WGS84'})
    transform: [1.0, 0.5, 0.0, 2.0, 0.0, -1.0]

How is this possible?

Is it possible to make the x and y coords lazy? Could you point me on a place in the code where this has been done already?
I'd like to have the lon and lat variables listed as optional coordinates, and also make them lazy. Any hint?

Thanks a lot for your help, I'm afraid this is going to need a few iterations ;-)

shoyer · 2017-02-10T21:28:02Z

The order of the variables and coordinates is random.

Most likely you're using a dict instead of an OrderedDict somewhere.

Is it possible to make the x and y coords lazy?

See here for an example of lazily computing arrays:

xarray/xarray/conventions.py

Line 314 in 6233320

class MaskedAndScaledArray(utils.NDArrayMixin):

There are other examples in that file, for decoding CF conventions.

I'd like to have the lon and lat variables listed as optional coordinates

Can you clarify what you mean by an optional coordinate?

fmaussion · 2017-02-10T21:37:45Z

Can you clarify what you mean by an optional coordinate?

Yes sorry, I meant the have them listed as coordinates without * instead of Data variables.

fmaussion · 2017-02-10T21:52:35Z

thanks @shoyer , I think I can work on this a bit further now and I'll get back to you if I have more questions.

shoyer · 2017-02-11T00:14:36Z

Getting optional coordinates is hard is because the current "backends" system in xarray was really only designed for handling netCDF-like libraries.

For example, every backend gets passed through xarray.conventions.decode_cf to use CF conventions to decode its data. You would actually need to add a "global attribute" specifying coords='lat lon' to set coordinates variables with this system.

One alternative would be to simply write a function for reading files with rasterio straight into xarray Dataset objects.

There are probably hacks that can get this working with data stores but this really calls out for the bigger refactor like #1087.

fmaussion · 2017-02-14T11:43:27Z

I made some progress with the lazy indexing, I'd be glad to have a first rough feedback on whether this is going in the right direction or not.

We have a decision to make regarding the API: I think that creating the optional lon and lat coords automatically is not a good idea:

in some cases, the x and y coordinates are already lons and lats and the 2D coords are obsolete
for big data files this is going to take ages and take a lot of memory
my initial idea to make them lazily evaluated might work, but in an ugly way: computing both lons and lats needs to be done in one single operation, and I'm not sure how this can be done in an elegant way
additionally, there is no way to make them show up as coordinates (as per @shoyer 's comment: Add RasterIO backend #1260 (comment))

The current implementation delegates this task to a utility function (get_latlon_coords_from_crs). It is currently very rasterio specific, but could be made more general -- API to be defined.

Thoughts?

shoyer · 2017-02-15T02:13:27Z

my initial idea to make them lazily evaluated might work, but in an ugly way: computing both lons and lats needs to be done in one single operation, and I'm not sure how this can be done in an elegant way

This could be done sanely with dask.array (example of multiple outputs from an element-wise function).

additionally, there is no way to make them show up as coordinates (as per @shoyer 's comment: #1260 (comment))

We could make this work, it just looks pretty hacky.

fmaussion · 2017-03-05T12:42:45Z

@shoyer thanks for the tips, I think that getting the lons and lats from dask is probably the most elegant method.

After trying various things I am still struggling with dask, and in particular on how to apply elemwise to functions like np.meshgrid (I get shape broadcasting errors). To get me started, I'd be grateful for an example on how to use dask to replace the code snippet below:

import xarray as xr
import numpy as np

ds = xr.DataArray(np.zeros((2, 3)),
                  coords={'x': np.arange(3), 'y': np.arange(2)},
                  dims=['y', 'x']).to_dataset(name='data')

# non-dask version
lon, lat = np.meshgrid(ds.x, ds.y)
ds['lon'] = (('y', 'x'), lon)
ds['lat'] = (('y', 'x'), lat)
ds.set_coords(['lon', 'lat'], inplace=True)
print(ds)

Thanks a lot!

…updates

…p.linspace now.

shoyer · 2017-05-31T01:03:44Z

xarray/backends/rasterio_.py

+        If True, cache data loaded from the underlying datastore in memory as
+        NumPy arrays when accessed to avoid reading from the underlying data-
+        store multiple times. Defaults to True unless you specify the `chunks`
+        argument to use dask, in which case it defaults to False. Does not


This last sentence can be safely removed.

shoyer · 2017-05-31T01:04:16Z

xarray/backends/rasterio_.py

+    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
+        DataArray into a dask array. This is an experimental feature; see the


We can probably remove the "experimental feature" bit. Dask is not going away.

shoyer · 2017-05-31T01:05:09Z

xarray/backends/rasterio_.py

+    lock : False, True or threading.Lock, optional
+        If chunks is provided, this argument is passed on to
+        :py:func:`dask.array.from_array`. By default, a per-variable lock is
+        used when reading data from netCDF files with the netcdf4 and h5netcdf


Does rasterio/GDAL need a thread lock? Let's remove netcdf4/h5netcdf specific references here and add sane defaults.

You cannot safely share a Rasterio or GDAL dataset object among threads, so a lock may be needed. See https://trac.osgeo.org/gdal/wiki/FAQMiscellaneous#IstheGDALlibrarythread-safe.

OK, I am on it.

@shoyer , I have trouble to understand what's going on in the open_dataset workflow. The documentation says: "By default, a per-variable lock is used when reading data from netCDF files with the netcdf4 and h5netcdf engines", but then the lock is decided once here. This will default to GLOBAL_LOCK and then passed to maybe_decode_store. Where is the "per variable" logic happening? Or is the GLOBAL_LOCK already doing this? (it is defined here)

(excuse my poor knowledge about these issues)

The documentation is out of sync here. We use a global lock for HDF5 (and should probably use one for Rasterio, too).

Now I understand much more what's going on, and also why salem was having troubles when adding diagnostic variables to netcdf objects. Can I update the docs all at once here?

shoyer · 2017-05-31T01:06:25Z

xarray/backends/rasterio_.py

+            import dask
+            if LooseVersion(dask.__version__) < LooseVersion('0.6'):
+                raise ImportError(
+                    'xarray requires dask version 0.6 or newer')


dask v0.6 is pretty old now, maybe better to replace this whole conditional logic with:

import dask from dask.base import tokenize

shoyer · 2017-05-31T15:51:24Z

Yes, please!

…

On Wed, May 31, 2017 at 8:48 AM, Fabien Maussion ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In xarray/backends/rasterio_.py <#1260 (comment)>: > + 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 + DataArray into a dask array. This is an experimental feature; see the + documentation for more details. + cache : bool, optional + If True, cache data loaded from the underlying datastore in memory as + NumPy arrays when accessed to avoid reading from the underlying data- + store multiple times. Defaults to True unless you specify the `chunks` + argument to use dask, in which case it defaults to False. Does not + change the behavior of coordinates corresponding to dimensions, which + always load their data from disk into a ``pandas.Index``. + lock : False, True or threading.Lock, optional + If chunks is provided, this argument is passed on to + :py:func:`dask.array.from_array`. By default, a per-variable lock is + used when reading data from netCDF files with the netcdf4 and h5netcdf Now I understand much more what's going on, and also why salem was having troubles when adding diagnostic variables to netcdf objects. Can I update the docs all at once here? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1260 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABKS1p7_JyTe8jdB1omZm_4V005U_ex_ks5r_YuygaJpZM4L92Pl> .

fmaussion · 2017-05-31T16:23:30Z

OK, all green.

Currently the rasterio tests are running on py36 only. Should I add rasterio to the other test suites as well?

shoyer · 2017-05-31T16:32:27Z

Currently the rasterio tests are running on py36 only. Should I add rasterio to the other test suites as well?

Yes, except for the one that tests bare-minimal dependencies.

gidden · 2017-06-01T07:53:59Z

Hey @fmaussion, is this ready for me to try out again? I wanted to let you and @shoyer iterate first.

…asterio

fmaussion · 2017-06-01T08:39:02Z

@gidden yes absolutely please give it a try, I think it's ready

fmaussion · 2017-06-01T09:53:24Z

@gidden I just updated the documentation recipe to use an accessor to compute the lons and lats, let me know what you think

shoyer

Looks great to me, just a couple more very minor points

shoyer · 2017-06-04T06:00:29Z

doc/gallery/plot_rasterio.py

+
+# Define the accessor
+@xr.register_dataarray_accessor('rasterio')
+class RasterioAccessor(object):


I would lean towards keeping this as just a function -- it keeps the example simpler and an accessor is probably overkill for this problem.

shoyer · 2017-06-04T06:03:23Z

xarray/backends/rasterio_.py

+        return out
+
+
+def open_rasterio(filename, chunks=None, cache=None, lock=False):


The default should be lock=None to use the global lock.

shoyer · 2017-06-04T06:04:09Z

xarray/backends/rasterio_.py

+        token = tokenize(filename, mtime, chunks)
+        name_prefix = 'open_rasterio-%s' % token
+        if lock is None:
+            lock = GLOBAL_LOCK


Actually, can you make a separate global lock for rasterio? We can safely use GDAL and HDF5 at the same time.

shoyer · 2017-06-04T06:05:51Z

xarray/tests/test_backends.py

@@ -1438,6 +1438,265 @@ class TestPyNioAutocloseTrue(TestPyNio):
    autoclose = True


+@requires_rasterio
+class TestRasterio(CFEncodedDataTest, Only32BitTypes, TestCase):


I don't think there's any advantage to inheriting from the other test methods, given that this isn't an actual datastore. I think this can just inherit from TestCase (then you can delete the first two methods).

shoyer · 2017-06-04T06:07:06Z

xarray/tests/test_backends.py

+                    with xr.open_dataarray(tmp_nc_file) as ncds:
+                        assert_identical(rioda, ncds)
+
+        # Create a geotiff file in latlong proj


If it's self-contained, make this a separate test method?

fmaussion · 2017-06-06T10:25:01Z

OK, let's get this one in and see what people will report about it.

Thanks @shoyer for your patience, @NicWayand and @jhamman for the original PR, @gidden for the testing/reviews and @sgillies for rasterio!

jhamman · 2017-06-06T16:43:30Z

Nice work @fmaussion. Thanks for bringing this home.

PeterDSteinberg · 2017-06-06T16:44:43Z

+1 @fmaussion - Thanks!

fmaussion mentioned this pull request Feb 14, 2017

add __getitem__ to rasterio datasets rasterio/rasterio#920

Closed

shoyer mentioned this pull request Mar 9, 2017

dask.array.meshgrid dask/dask#2057

Closed

fmaussion force-pushed the feature-rasterio branch from cc80626 to d877826 Compare May 3, 2017 20:25

Joe Hamman and others added 19 commits May 18, 2017 00:49

rasterio checkin

1643ce5

temp fixes

067dedb

update rasterio reader, no lazy loading, no decoding of coords

7906cfd

keep band dim even for single band. Fix longitude typo

2e1b528

Fix lat/lon decoding. Remove requirment comment

532c5d3

Attr error suppression. DataArray to Variable objects. CI requirment …

8bc3da3

…updates

remove >= requirment

77cc0ca

added conda-forge channel to CI check

776cbd9

add scipy requirement

eb739de

roll back ci requirements. Rename vars

3a394ae

roll back

061b8fd

fixed coord spacing bug where x and y were +1 dim than raster. Uses n…

c0962fa

…p.linspace now.

change test env to py36

7275ffa

first tests

51a60af

other tests

5155634

fix test

09196ee

get the order right

e2b6786

some progress with indexing

228a5a3

cosmetic changes

4b57d1c

Final tweaks

223ce0c

shoyer reviewed May 31, 2017

View reviewed changes

fmaussion added 2 commits May 31, 2017 18:10

Lock-doc tweaks

6cf2ce9

Merge branch 'master' into feature-rasterio

2cd0386

fmaussion added 2 commits June 1, 2017 10:37

Add rasterio to other test suites

9193a2b

Merge remote-tracking branch 'origin/feature-rasterio' into feature-r…

4bf4b6a

…asterio

fmaussion added 2 commits June 1, 2017 11:07

use context managers in tests for windows

c778948

Change example to use an accessor

4299957

shoyer reviewed Jun 4, 2017

View reviewed changes

fmaussion added 3 commits June 5, 2017 11:08

Reviews

fcdd894

Merge branch 'master' into feature-rasterio

1ca6e38

typo

d5c964e

shoyer approved these changes Jun 5, 2017

View reviewed changes

shoyer mentioned this pull request Jun 5, 2017

v0.9.6 release #1441

Closed

fmaussion merged commit 48c7164 into pydata:master Jun 6, 2017

PeterDSteinberg mentioned this pull request Jun 6, 2017

Review / modify earthio/tif.py as it relates to xarray rasterio ContinuumIO/earthio#17

Closed

perrygeo mentioned this pull request Sep 11, 2017

rasterio as a xarray backend rasterio/rasterio#576

Closed

forman mentioned this pull request Sep 25, 2017

Support for GeoTIFF CCI-Tools/cate#173

Open

fmaussion mentioned this pull request Nov 13, 2017

Use rasterio's transform instead of homemade coordinates #1712

Merged

4 tasks

rabernat mentioned this pull request Apr 9, 2018

Anyone working on a to_tiff? Alternatively, how do you write an xarray to a geotiff? #2042

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add RasterIO backend #1260

Add RasterIO backend #1260

fmaussion commented Feb 10, 2017 •

edited

Loading

fmaussion commented Feb 10, 2017

shoyer commented Feb 10, 2017

fmaussion commented Feb 10, 2017

fmaussion commented Feb 10, 2017

shoyer commented Feb 11, 2017

fmaussion commented Feb 14, 2017

shoyer commented Feb 15, 2017

fmaussion commented Mar 5, 2017

shoyer May 31, 2017

shoyer May 31, 2017

shoyer May 31, 2017

sgillies May 31, 2017

fmaussion May 31, 2017

shoyer May 31, 2017

fmaussion May 31, 2017

shoyer May 31, 2017

shoyer commented May 31, 2017 via email

fmaussion commented May 31, 2017

shoyer commented May 31, 2017

gidden commented Jun 1, 2017

fmaussion commented Jun 1, 2017

fmaussion commented Jun 1, 2017

shoyer left a comment

shoyer Jun 4, 2017

shoyer Jun 4, 2017

shoyer Jun 4, 2017

shoyer Jun 4, 2017

shoyer Jun 4, 2017

fmaussion commented Jun 6, 2017

jhamman commented Jun 6, 2017

PeterDSteinberg commented Jun 6, 2017

		return out


		def open_rasterio(filename, chunks=None, cache=None, lock=False):

Add RasterIO backend #1260

Add RasterIO backend #1260

Conversation

fmaussion commented Feb 10, 2017 • edited Loading

fmaussion commented Feb 10, 2017

shoyer commented Feb 10, 2017

fmaussion commented Feb 10, 2017

fmaussion commented Feb 10, 2017

shoyer commented Feb 11, 2017

fmaussion commented Feb 14, 2017

shoyer commented Feb 15, 2017

fmaussion commented Mar 5, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

shoyer commented May 31, 2017 via email

fmaussion commented May 31, 2017

shoyer commented May 31, 2017

gidden commented Jun 1, 2017

fmaussion commented Jun 1, 2017

fmaussion commented Jun 1, 2017

shoyer left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fmaussion commented Jun 6, 2017

jhamman commented Jun 6, 2017

PeterDSteinberg commented Jun 6, 2017

fmaussion commented Feb 10, 2017 •

edited

Loading