Rasterio: Wrap ndarry, transform, and CRS into an object

[brendan-ward]

There are a number of places where a rasterio operation takes as input a raster or raster band, OR a numpy ndarray, transform, and sometimes CRS. This complicates the function signatures a bit, as it means having parameters that are meaningful in some cases but not others, and it lacks a bit of cohesion between related entities that may be helpful in many cases.

Perhaps it is time to create a lightweight object (tuple or class) that wraps these together, so that we can pass just one thing around and have it be more equivalent to passing a raster / raster band?

Something like the Raster object in @perrygeo 's python-rasterstats package, but not wrapping rasters / bands; only ndarrays.

This surfaced in #578, among others.

[brendan-ward]

Putting this in a class rather than a tuple may make it easier to implement useful methods like index and window like we have now on rasters.

[perrygeo]

It would be ideal to have an object with an interface like the DatasetReader but without IO. Then functions could accept either an open dataset or this raster object itself and behave in roughly the same way. I haven't wrapped my head around exactly how we would be implement this.

A lightweight named tuple could also provide some value here without being quite as big of a lift.

I'd say that however it is implemented, we should aim to keep the DatasetReader API stable. Since it would be an addition to the API, not a breaking change, I'd like to propose punting this to post-1.0? Is that good with you, @brendan-ward ?

[brendan-ward]

@perrygeo I was thinking this would lead to breaking changes, wherein we standardize on this instead of passing ndarrays, transform, CRS and related parameters as inputs to functions when we are not passing in an array. I realize this has sweeping implications for documentation, tests, etc.

The more I've thought about this, the more I think that we should define the minimal API this object should provide, and implement it as a class so that we have the appropriate methods and properties at hand. I don't think a tuple gets us enough here to be as useful as I'd like.

For example, it seems natural that this should have a bounds property and thus fulfill part of the API needed for something like your recent changes in #795 .

I'd like to avoid having to do a round of deprecation warnings on this, because that would make things worse rather than better (now there are 3 ways to pass in inputs to functions instead of 2), though if we must, we must.

If you and @sgillies are open to this idea at all, I will sketch out a proposal for what this could look like, then we can carefully consider where it fits in before or after 1.0?

[brendan-ward]

Very rough sketch of the interface. Not sure about a few properties / methods. Thinking of calling this GeoArray

Constructor:
(ndarray, transform, crs=None)

Properties:

• transform
• crs
• shape (from ndarray)
• bounds
• res
• nodata (from ndarray fill_value if masked)?
• nodatavals?
• ul?
• height?
• width?
• count?
• indexes?

Methods:

• index
• window
• window_transform
• window_bounds
• lnglat
• get_transform?
• get_crs?

Will start stubbing this out in a PR to see what this might actually look like.

[sgillies]

@brendan-ward I remember seeing another Python project for georeferenced ndarrays, but am unable to track it down at the moment.

[brendan-ward]

@sgillies I didn't find much, and nothing was close to the API we'd like here. Did you mean https://github.com/schaefed/geoarray which you referenced in #578?

[perrygeo]

There's also https://github.com/ozak/georasters ? Many others have done similar things, I'd imagine there lots of prior art here.

It sounds like we're going for basically an object that implements the DatasetReader (Soon to be RasterDataReader) interface but doesn't have a read() method. This might imply a couple things

• RasterData would be a logical classname (RasterDataReader minus the read)
• They could share a common base class. Many of the non-IO methods/properties could be written in python and be inherited by both RasterDataReader and RasterData. This would be a shining example of the usefulness of the public io refactor.

[sgillies]

If you don't mind, I'm going to put this on the post-1.0 milestone.

[brendan-ward]

@sgillies I don't want to hold up 1.0, but this might have significant API implications. Mostly around reducing the number of parameters we pass around. Is there a post-1.0 milestone we should be targeting for other significant API changes?

[perrygeo]

I'm convinced that we can implement this without any API changes - in fact if this new object is intended to implement the same interface as standard datasets, the lack of API changes is implied.

What breaking changes do you anticipate @brendan-ward ?

There are the plotting functions but those could be safely deprecated to use our "pure-function plus dataset method" pattern. IOW I think any breaking changes involved in this class could be safely handled by deprecation until 2.0.

[sgillies]

Taking reproject as an example

@ensure_env
def reproject(
        source,
        destination,
        src_transform=None,
        src_crs=None,
        src_nodata=None,
        dst_transform=None,
        dst_crs=None,
        dst_nodata=None,
        resampling=Resampling.nearest,
        init_dest_nodata=True,
        **kwargs):

Because we've used keyword arguments for everything except the source and destination rasters, we don't need to break the API. Same goes for methods in features and fill, which is good.

The output of merge() on the other hand,

    -------
    tuple

        Two elements:

            dest: numpy ndarray
                Contents of all input rasters in single array.

            out_transform: affine.Affine()
                Information for mapping pixel coordinates in `dest` to another
                coordinate system

could be problematic.

[perrygeo]

Ah, yes if we're _returning_ instances of this class instead of arrays or array+something tuples, that would change the API. The mere existence of this class, making it available for developers who want it, does not.

[brendan-ward]

I was thinking that if we did this, that we could drop the *_transform, *_crs parameters from the reproject function signature (hence the breaking changes), rather than making them optional. They are only there for the case where we pass in ndarrays, and I'm not really keen on adding a 3rd variant to those: raster band, ndarray w/ these params, or geoarray w/o those params.

We could provide deprecation of those parameters for a while if we'd like, but that means handling more rather than less variants.

Perhaps I'm not thinking about this correctly? My goal was more encapsulation of related properties, so that we could make ndarrays more like rasters, and have somewhat fewer cases where we have to handle them fundamentally differently at the API level for rasterio (inside the functions we'll have to handle them differently, for sure).

[perrygeo]

make ndarrays more like rasters

That sounds right to me. But there's a big difference between providing a wrapper class and changing the API across the board to use that class.

I think it's important to think about this wrapper class as something that will behave more or less exactly like dataset objects. So if we're trying to implement the dataset object interface, we should be aiming for the same API with minimal changes.

The discussion of changing function signatures to accept and return dataset objects is a separate concern. I'm -1 on that for reasons outline in https://github.com/mapbox/rasterio/issues/803#issuecomment-226761437

[brendan-ward]

@perrygeo maybe it is just me, but I see a big difference between what we are suggesting here (primary inputs are rasters or ndarrays), versus those indicated in #803 (primary inputs are properties that can be managed independently - bounds, crs, transform). The latter benefit a lot from the approach outlined in #803; the former I think benefit a lot from the approach outlined here.

For instance, reproject above doesn't really benefit from having *_crs, *_transform as specific parameters that are only used in some cases (when other inputs are ndarrays).

[perrygeo]

Good points, this needs more detailed discussion for sure. I submit that we really have two PRs here:

1. Creating an array wrapper that implements the dataset object interface
2. Changing function signatures to accept objects which implement the dataset object interface

I think it's important to think about them as distinct concerns since they can be developed completely independently. I don't want concerns about (2) to hold up development of (1) or vice versa.

[brendan-ward]

@perrygeo excellent suggestions. Let's focus on 1 and prove it helps, then do 2. Thank you for helping break these down into more logical steps.

[brendan-ward]

Sorry for dragging this out, I have very little time for rasterio right now and this is a bit involved.

First up - for showing value of this for reproject, I think I can refactor _reproject to use InMemoryRaster (which should simplify code there), then that should help inform how we want GeoArray and InMemoryRaster to work together in the Cython layer.

[brendan-ward]

@sgillies given that we have not made any progress on this in quite some time, are we still interested in pursuing this?

[jdmcbr]

FWIW, this is of interest to me, though I don't anticipate having enough time in the next three months to lead this.

[sgillies]

I've moved this to the project board and will give it some attention in the next week.

[sgillies]

@brendan-ward @perrygeo I think I see a plan. How about we break this into 2 other tasks:

1. Define dataset interfaces (I'm thinking 3: one with metadata, one reader, one writer) as abstract base classes.
2. Write the wrapper(s) for arrays + anciliary data for internal use only at 1.0.

This is more or less what @perrygeo proposed, but with a formal interface and no need to change method signatures in a breaking way.

The thing we'll do to make the formal interface stick is assert it as soon as possible in our methods:

def do_rasterio_stuff(source, src_crs=None, src_transform=None):
    if isinstance(source, np.array):
        # wrap it up
        source = wrap(source, src_crs, src_transform)

    assert isinstance(source, DatasetABC)

    # After this, only call defined methods.

Or something like that.

What do you think?

[brendan-ward]

@sgillies what do you mean by metadata?

One of the motivations was to reduce the number of parameters we had to pass in where our source is a numpy array, so it seems like we would get the most benefit of this if it were more public than internal use only. Perhaps I don't understand what you have in mind?

Otherwise, in theory this sounds like a reasonable path forward and reduces breaking changes.

[sgillies]

@brendan-ward by "metadata" I mean an interface that would provide crs, transform, count, shape, &c but no reading or writing.

It could be that good documentation for the wrapper wouldn't hold up 1.0, so maybe we could make it public. It's more important to me that the interface be done by 1.0. By adding the ABCs to Rasterio we also enable developers to make conforming wrappers for objects from their own applications and libraries. I think this should help some efforts in the Dask project.

[perrygeo]

By adding the ABCs to Rasterio we also enable developers to make conforming wrappers for objects from their own applications and libraries.

This could be a huge win, I could even imagine a full Dataset implementation that didn't rely on GDAL. For the near-term, the advantages would be more for interface documentation and a clean class hierarchy, correct?

The wrapper would be a great feature but conceptually distinct from the abstract base class. I think the ABC could be implemented first (for purposes of defining the 1.0 interface) and the wrapper implementation could come later.

[sgillies]

After some reading and experimenting, I'm going to shelve the abstract base class I wrote until we're ready to drop support for Python 2.7. There's enough of a difference in the ABC API between 2.7 and 3.5/3.6 that I don't think it's worth trying to do both.

I'm going to close this issue and reopen it down the road. I feel that other issues on the board are more important to address for 1.0 than this one.

The consequence: in 1.0, we will require users to "unpack" their application's objects into separate array, crs, transform/gcp objects before calling functions like warp.reproject().

[sgillies]

Copying the ABCs here for the future:

"""Rasterio's abstract base classes"""

import abc


class GeoreferencedArray(abc.ABC):

    @property
    @abc.abstractmethod
    def crs(self):
        """The coordinate reference system for the array

        Returns
        -------
        rasterio.crs.CRS
        """

    @property
    @abc.abstractmethod
    def bounds(self):
        """The spatial bounding box of the array

        Returns
        -------
        lower left x, lower left y, upper right x, upper right y : float
        """

    @property
    @abc.abstractmethod
    def data(self):
        """The n-d array

        Returns
        -------
        numpy.array
        """


class AffineGeoreferencedArray(GeoreferencedArray):

    @property
    @abc.abstractmethod
    def transform(self):
        """Affine transformation determining array georeferencing

        Returns
        -------
        rasterio.affine.Affine
        """


class GCPGeoreferencedArray(GeoreferencedArray):

    @property
    @abc.abstractmethod
    def gcps(self):
        """Ground control points determining array georeferencing

        Returns
        -------
        list of rasterio.control.GroundControlPoint
        """

[astrojuanlu]

After some reading and experimenting, I'm going to shelve the abstract base class I wrote until we're ready to drop support for Python 2.7.

Is there a roadmap for dropping Python 2 in rasterio? I looked in the issues list and didn't see anything.

Copying the ABCs here for the future:

We've been working on similar APIs in our recently published package telluric, which builds on top of rasterio among others http://telluric.readthedocs.io/en/latest/api.html#telluric.georaster.GeoRaster2 If you think there's anything we can do to contribute back, please ping us :)