Feature/rasterio #1070

NicWayand · 2016-10-31T16:14:55Z

@jhamman started a backend for RasterIO that I have been working on. There are two issues I am stuck on that I could use some help:

Lat/long coords are not being decoded correctly (missing from output dataset). Lat/lon projection are correctly calculated and added here (https://github.com/NicWayand/xray/blob/feature/rasterio/xarray/backends/rasterio_.py#L117). But, it appears (with my limited knowledge of xarray) that the lat/long coords contained within obj are lost at this line (https://github.com/NicWayand/xray/blob/feature/rasterio/xarray/conventions.py#L930).
Lazy-loading needs to be enabled. How can I setup/test this? Are there examples from other backends I could follow?

#790

fmaussion · 2016-10-31T16:27:59Z

This is awesome!

For 2: I think you don't have much more to do than make __getitem__ use the indexer to read the data with a window (example here)

shoyer

You don't need to do anything special to enable lazy loading. That logic is all handled in the xarray.Variable object.

shoyer · 2016-10-31T16:26:06Z

xarray/backends/rasterio_.py

+            except AttributeError:
+                pass
+
+        self.coords = _try_to_get_latlon_coords(self.coords, self._attrs)


You need to return these in the dict returned by get_variables(). DataStore objects don't have a coords attribute.

Thanks. This now pulls in lat/long if the projection was successful, otherwise it just returns the basic x/y indices. Should we also be returning the original projection coords? (currently they are lost).

shoyer · 2016-10-31T16:27:15Z

ci/requirements-py27-cdat+pynio.yml

@@ -3,7 +3,7 @@ channels:
  - conda-forge
 dependencies:
  - python=2.7
-  - cdat-lite
+  # - cdat-lite


This seems unrelated?

Yes, it has been removed.

shoyer

It looks like tests still need some work to pass on Travis-CI.

shoyer · 2016-10-31T17:55:57Z

xarray/backends/rasterio_.py

+        coords = dict(y=coords['y'], x=coords['x'])
+        dims = ('y', 'x')
+
+        coords_out['lat'] = DataArray(


These should be just xarray.Variable objects, not DataArray objects

shoyer · 2016-10-31T18:02:00Z

xarray/backends/rasterio_.py

+        proj = pyproj.Proj(attrs['crs'])
+        x, y = np.meshgrid(coords['x'], coords['y'])
+        proj_out = pyproj.Proj("+init=EPSG:4326", preserve_units=True)
+        xc, yc = _transform_proj(proj, proj_out, x, y)


You might try to do this calculation lazily, e.g., by making a utils.NDArrayMixin subclass in conventions.py and wrapping with indexing.LazilyIndexedArray.

Sorry I don't understand how to do this. Is there an example I can go off of?

We have some examples of making NDArrayMixin subclasses in conventions.py, e.g.,
https://github.com/pydata/xarray/blob/master/xarray/conventions.py#L311

Another option would be to use dask.array to build the coordinate arrays, though the NDArrayMixin would avoid the need for dask.

shoyer · 2016-10-31T18:02:15Z

xarray/backends/rasterio_.py

+            except AttributeError:
+                pass
+
+    # def get_vardata(self, var_id=1):


should be removed

shoyer · 2016-10-31T18:02:51Z

xarray/backends/rasterio_.py

+
+        self._attrs = OrderedDict()
+        for attr_name in ['crs', 'transform', 'proj']:
+            try:


use with pycompat.suppress(AttributeError):

shoyer · 2016-10-31T18:04:46Z

xarray/backends/rasterio_.py

+def _try_to_get_latlon_coords(coords, attrs):
+    coords_out = {}
+    try:
+        import pyproj


I'm a little surprised rasterio doesn't have projections built in.

They do: https://mapbox.github.io/rasterio/topics/reproject.html

I will take a stab at doing it the rasterio way.

That's a reprojection, which does seem useful, but it seems separate from giving you lat and lon coordinates

shoyer · 2016-10-31T18:06:43Z

ci/requirements-py35.yml

@@ -11,6 +11,7 @@ dependencies:
  - pandas
  - seaborn
  - scipy
+  - rasterio


Unless rasterio only supports Python 3, should add it to other test suites, too (especially Python 2.7).

Rasterio supports 2.7 and 3.3-3.5 (https://mapbox.github.io/rasterio/)

I added this requirement to all other test suites, was that correct? Version 1.0 is needed as many function names changed in rasterio at V1.0.

NicWayand · 2016-10-31T18:10:52Z

Tested open_mfdataset() on 100+ geotiffs and lazyloading with rasterio does appear to be working.

…updates

NicWayand · 2016-10-31T18:51:46Z

Travis-ci fail is because it can't find rasterio, which comes through the conda-forge channel (https://github.com/conda-forge/rasterio-feedstock). I think it needs to be added as described here (http://conda.pydata.org/docs/travis.html#additional-steps). But I am new to Travis-ci, so don't want to mess up the current .travis.yml file.

ocefpaf · 2016-10-31T18:55:28Z

ci/requirements-py27-cdat+pynio.yml

@@ -12,5 +12,6 @@ dependencies:
  - scipy
  - pytest-cov
  - cyordereddict
+  - rasterio>=1.0


There is no official release of rasterio >=1.0. However, I did merge the latest alpha, so if you remove the >= you should get it.

shoyer · 2016-10-31T18:56:57Z

Add conda-forge it to the channels section of conda requirements yml files:

xarray/ci/requirements-py27-cdat+pynio.yml

Line 3 in f40d323

- conda-forge

NicWayand · 2016-10-31T19:52:56Z

Any idea why Segmentation faults occur for 3.4 and 4.5?

shoyer · 2016-10-31T19:55:54Z

I would start by adding -v to our invocation of py.test in .travis.yml:

xarray/.travis.yml

Line 72 in a247d93

- py.test xarray --cov=xarray --cov-report term-missing

fmaussion · 2016-10-31T21:02:45Z

xarray/backends/rasterio_.py

+class RasterioArrayWrapper(NDArrayMixin):
+    def __init__(self, ds):
+        self._ds = ds
+        self.array = ds.read()


Unless the behavior of rasterio changed with v1.0, this loads the data in memory. I might be wrong, but I think the call to ds.read() should happen in __getitem__, ideally with a window kwarg.

I also wonder if the call to read shouldn't happen in a with rasterio.Env(): environment.

@shoyer, do we need to maintain the array attribute here? Would it make more sense to just populate the _ds attribute and set the array in __getitem__?

We don't need to set array -- that's merely for the benefit of NDArrayMixin in core/utils.py. You might even skip NDArrayMixin altogether and just use NdimSizeLenMixin, though you'll want to define properties/methods for each of those defined by NDArrayMixin.

@NicWayand -

I think you can basically drop this in replacing the current RasterioArrayWrapper:

class RasterioArrayWrapper(NdimSizeLenMixin): """Mixin class RasterIO datasets for making wrappers of N-dimensional arrays that conform to the ndarray interface required for the data argument to Variable objects. A subclass should set the `array` property and override one or more of `dtype`, `shape` and `__getitem__`. """ def __init__(self, ds): self.ds = ds @property def dtype(self): if len(set(self.ds.dtypes[0])) != 1: raise ValueError( 'Can only handle Rastio dataset with all bands having the same type') return np.dtype(self.ds.dtypes[0]) @property def shape(self): return self.ds.shape def __array__(self, dtype=None): '''Not sure if this will work as is''' return np.asarray(self[...], dtype=dtype) def __getitem__(self, key): band = range(self.shape[0])[key[0]] window = [] for win in key[1:]: if instance(win, slice): window.append((win.start, win.stop)) elif isinstance(win, int): window.append((win, win + 1)) else: # integer ndarray window.append((win.min(), win.max())) raw_data = self.ds.read(band, window=window) # now, fix up raw_data to conform to numpy indexing conventions # - drop axes for integer band/windows # - stride if windows are slices with win.step != 1 # - subset if window is an integer ndarray def __repr__(self): return '%s(array=%r)' % (type(self).__name__, self.ds)

I'm not sure the __array__ will work as it is until we get the __get_item__ fully functional.

It's also worth consider how/if you want to handle automatically masking missing values -- it looks like ds.read() can optionally return a MaskedArray. The standard xarray approach would be to automatically promote the dtype and fill these with NaNs.

Just to confirm, the above code is waiting for the decision from rasterio group?

I would suggest getting something working here, with an eye toward keeping it self-contained. Later, we can try to port it upstream to rasterio.

jhamman · 2016-10-31T21:52:49Z

xarray/backends/rasterio_.py

+    def __getitem__(self, key):
+        if key == () and self.ndim == 0:
+            return self.array.get_value()
+        return self.array[key]


Based on what @fmaussion said above, I think this should be something like:

def __getitem__(self, key): if key == () and self.ndim == 0: return self._ds.read() return self._ds.read(band, window=window)

Where band and window describe the slice of data to be read. Off the top of my head, I'm not exactly sure how to parse the key here though.

By the time indexers get here, they should have already passed through canonicalize_indexer, which means you only need to handle a key that is a tuple of the appropriate number of dimensions (i.e., 3) composed of integers, slices and integer ndarrays.

Based on the docstring for read, we want something like:

def __getitem__(self, key): band = range(self.shape[0])[key[0]] window = [] for win in key[1:]: if instance(win, slice): window.append((win.start, win.stop)) elif isinstance(win, int): window.append((win, win + 1)) else: # integer ndarray window.append((win.min(), win.max())) raw_data = self._ds.read(band, window=window) # now, fix up raw_data to conform to numpy indexing conventions # - drop axes for integer band/windows # - stride if windows are slices with win.step != 1 # - subset if window is an integer ndarray

Honestly, this logic should probably live in rasterio if possible. I'm a little surprised that they have never implemented a __getitem__ method.

See also: https://gist.github.com/lpinner/bd57b54a5c6903e4a6a2 (can't reuse this directly, though, because it doesn't have a license).

Anyways, I would definitely see if the rasterio folks are up for implementing a __getitem__ method. There's no reason why this is xarray specific -- you would need this for just using dask.array with rasterio, as well.

See also rasterio/rasterio#920 for related discussion.

jhamman

@NicWayand - this is coming along nicely. The remaining issues as I see them are:

unit tests
__getitem__ wrapper method
coordinate projection

jhamman · 2016-11-01T15:08:44Z

ci/requirements-py27-min.yml

 dependencies:
  - python=2.7
  - pytest
  - numpy==1.9.3
  - pandas==0.15.0
+  - rasterio


I think we can roll back these changes. This test setup is just for required dependencies which rasterio will not be.

jhamman · 2016-11-01T15:09:08Z

ci/requirements-py27-netcdf4-dev.yml

@@ -8,6 +10,7 @@ dependencies:
  - numpy
  - pandas
  - scipy
+  - rasterio


I think we can also roll back these changes as well

jhamman · 2016-11-01T15:09:54Z

ci/requirements-py33.yml

 dependencies:
  - python=3.3
  - pytest
  - pandas
+  - rasterio


roll back these changes

jhamman · 2016-11-01T15:13:24Z

ci/requirements-py27-pydap.yml

@@ -8,6 +10,7 @@ dependencies:
  - numpy
  - pandas
  - scipy
+  - rasterio


I don't think we need to test rasterio along with pydap so we can roll back these changes.

jhamman · 2016-11-01T15:13:43Z

ci/requirements-py35-dask-dev.yml

@@ -7,6 +9,7 @@ dependencies:
  - pandas
  - scipy
  - toolz
+  - rasterio


roll back these changes

jhamman · 2016-11-01T15:13:56Z

ci/requirements-py35-pandas-dev.yml

@@ -7,6 +9,7 @@ dependencies:
  - netcdf4=1.1.9
  - scipy=0.16.0
  - toolz
+  - rasterio


roll back these changes

jhamman · 2016-11-01T15:19:42Z

xarray/backends/rasterio_.py

+class RasterioArrayWrapper(NDArrayMixin):
+    def __init__(self, ds):
+        self._ds = ds
+        self.array = ds.read()


@shoyer, do we need to maintain the array attribute here? Would it make more sense to just populate the _ds attribute and set the array in __getitem__?

jhamman · 2016-11-01T15:21:52Z

xarray/backends/rasterio_.py

+        self._attrs = OrderedDict()
+        with suppress(AttributeError):
+            for attr_name in ['crs', 'transform', 'proj']:
+                    self._attrs[attr_name] = getattr(self.ds, attr_name)


I think there is some extra indentation here -- should be 4 spaces.

jhamman · 2016-11-01T15:23:00Z

xarray/backends/rasterio_.py

+        coords = _try_to_get_latlon_coords(self.coords, self._attrs)
+        vars = {__rio_varname__: self.open_store_variable(__rio_varname__)}
+        vars.update(coords)
+        return FrozenOrderedDict(vars)


we should use another name other than vars since it is a python builtin.

jhamman · 2016-11-01T15:31:37Z

xarray/test/test_backends.py

+            on_disk = open_dataset(tmp_file, engine='rasterio')
+            actual = on_disk.rename({'foo': 'bar', 'x': 'y'})
+            del on_disk  # trigger garbage collection
+            self.assertDatasetIdentical(actual, expected)


@NicWayand - take a look at the PyNio tests. They are probably the closest analog for what we need to test here.

shoyer · 2016-11-01T15:49:30Z

.travis.yml

@@ -69,7 +69,7 @@ install:
  - python setup.py install

 script:
-  - py.test xarray --cov=xarray --cov-report term-missing
+  - py.test xarray --cov=xarray --cov-report term-missing -v


This change can actually be permanent -- it's valuable for tracking down segfaults on Travis more broadly.

…p.linspace now.

fmaussion · 2017-02-10T13:22:02Z

Hi, unless there are objections I'd like to try to finish this PR. I just forked @NicWayand 's branch and started to look at it, and I hope to get something done during next week or so...

NicWayand · 2017-02-10T17:54:13Z

Hi @fmaussion, no objections here. I got it working just barely for my project, and won't have time in the near future to devote to wrap this up.

byersiiasa · 2017-05-22T07:50:55Z

@gidden you might be interested in this

gidden · 2017-05-22T08:28:53Z

@byersiiasa indeed I am, in fact I'm reviewing #1260 which overtook this one =)

fmaussion · 2017-05-22T08:47:40Z

Closing this in favor of #1260

Joe Hamman and others added 5 commits April 5, 2016 16:34

rasterio checkin

b6080e1

Merge branch 'master' of github.com:pydata/xarray into feature/rasterio

c6fb2b5

temp fixes

abeff81

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

be6094e

keep band dim even for single band. Fix longitude typo

9855d32

shoyer reviewed Oct 31, 2016

View reviewed changes

Fix lat/lon decoding. Remove requirment comment

2a94fa4

shoyer reviewed Oct 31, 2016

View reviewed changes

Nicolas Wayand added 2 commits October 31, 2016 12:31

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

dd71a45

…updates

missed ci/requirements-py35-pandas-dev.yml

1db0533

ocefpaf reviewed Oct 31, 2016

View reviewed changes

Nicolas Wayand added 3 commits October 31, 2016 12:57

remove >= requirment

e094bec

added conda-forge channel to CI check

21f4343

add scipy requirement

b11afaa

temp add verbose output

af902f0

fmaussion reviewed Oct 31, 2016

View reviewed changes

jhamman reviewed Oct 31, 2016

View reviewed changes

jhamman mentioned this pull request Nov 1, 2016

add __getitem__ to rasterio datasets rasterio/rasterio#920

Closed

jhamman requested changes Nov 1, 2016

View reviewed changes

jhamman mentioned this pull request Nov 1, 2016

rasterio as a xarray backend rasterio/rasterio#576

Closed

shoyer reviewed Nov 1, 2016

View reviewed changes

Nicolas Wayand and others added 3 commits November 2, 2016 13:07

roll back ci requirements. Rename vars

0369ef8

roll back

bb73429

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

7037e9d

…p.linspace now.

merge

3e90b99

gidden mentioned this pull request Jan 4, 2017

support multiple bands in rasterio read robintw/XArrayAndRasterio#2

Open

fmaussion mentioned this pull request Feb 10, 2017

Add RasterIO backend #1260

Merged

6 tasks

forman mentioned this pull request Mar 30, 2017

Support for GeoTIFF CCI-Tools/cate#173

Open

fmaussion closed this May 22, 2017

Feature/rasterio #1070

Feature/rasterio #1070

Conversation

NicWayand commented Oct 31, 2016

fmaussion commented Oct 31, 2016 • edited Loading

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

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

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

NicWayand Oct 31, 2016 • edited Loading

Choose a reason for hiding this comment

NicWayand commented Oct 31, 2016

NicWayand commented Oct 31, 2016

Choose a reason for hiding this comment

shoyer commented Oct 31, 2016

NicWayand commented Oct 31, 2016

shoyer commented Oct 31, 2016

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

Choose a reason for hiding this comment

shoyer Oct 31, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jhamman 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

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 Feb 10, 2017

NicWayand commented Feb 10, 2017

byersiiasa commented May 22, 2017

gidden commented May 22, 2017

fmaussion commented May 22, 2017

fmaussion commented Oct 31, 2016 •

edited

Loading

NicWayand Oct 31, 2016 •

edited

Loading

shoyer Oct 31, 2016 •

edited

Loading