Grid features

[mkborregaard]

Hi, when I do geographical stuff, our sampling units are very frequently laid out as a spatial grid. This can be represented very compactly (e.g. https://github.com/EcoJulia/SpatialEcology.jl/blob/mkb/dev/src/DataTypes.jl#L19-L26 / https://github.com/EcoJulia/SpatialEcology.jl/blob/mkb/dev/src/DataTypes.jl#L40-L43 or https://github.com/FugroRoames/SpatialGrids.jl (where they call it a "raster"). I think of this as something distinct from a Raster, which is a Matrix holding a single value in each cell, and essentially we treat them as square polygons.

Does such a type have a place in your model here?

[mauro3]

Would it make sense to use a range instead of hand-coding the start, step, length?

[mkborregaard]

Yes.
But maybe I should clarify - my question was not "do you want this implementation?" but rather "should we include an Abstract Interface to this?". I'm a little unclear how to use this package, just very keen for the spatial types I define to play well with the JuliaGeo ecosystem, and this seems like the entry point.

[visr]

Thanks for bringing this up. Didn't think about this that much yet, but it seems like a good addition to this package.

First we must decide what an AbstractRaster (name tbd) is, to be able to make an interface. Some points:

• An AbstractRaster is a AbstractArray
• An AbstractRaster provides a mapping from indices to x-y coordinates (think CoordinateTransformations.AffineMap)
• An AbstractRaster optionally holds a coordinate reference system

Does that make sense so far, what am I missing? Is there a basis that we can use for this? Would be good to at least be partly compatible with https://www.gdal.org/gdal_datamodel.html.

[mkborregaard]

I think that looks quite complete in terms of base attributes. But, is it worth thinking of grid and raster as different things (as I suggested above)? They classically are, eg in R there's a raster which is an (out-of-core) numeric Matrix with a spatial mapping, and SpatialGrid which is like an AbstractFeature where the AbstractGeometry is SpatialPixels (a Point shapefile where points are regularly spaced on a grid pattern), and each entry/point links into a dataframe row. (I have some code, which might not be very good, for transforming regular spatial points to a grid here: https://github.com/EcoJulia/SpatialEcology.jl/blob/master/src/Constructor_helperfunctions.jl#L135-L168)

But on the other hand maybe not? That distinction derives from the concrete implementation in traditional GIS systems, where the representations of rasters and shapefiles are very different. Since a Julia Matrix can hold anything it might be possible to represent these things the same way? Not sure here.

[meggart]

I think that looks quite complete in terms of base attributes. But, is it worth thinking of grid and raster as different things (as I suggested above)?

I would not make this distinction here. As you already mentioned, the difference is more an implementation detail than a conceptual one. I think it should be possible to implement e.g. an IterableTable interface (https://github.com/queryverse/IterableTables.jl) for a raster Matrix object that simply loops over each grid point and the associated coordinates. On the other hand, it should be possible to make a SpatialDataFrame like object behave as if it was a Spatial matrix.

So my suggestion would be to define an Interface for AbstractRaster similar to what @visr suggested and then allow different implementations, one in which the actual data is stored as a plain 2d array with associated grid metadata, and one where the data is stored as a DataFrame-like object with data, lon and lat values as columns.

[visr]

See also yeesian/ArchGDAL.jl#37 (comment) for some related discussion. Main issue; if we have an AbstractRaster that is a AbstractArray, how do we know which dimension is which? For GDAL it is (col, row, band), though the (col, row) ordering might be confusing. And how many dimensions should we support? The simplest answer might be fixing it to a 2D (row, col), but that's very restricting. It would also be nice to support GDAL rasters without transposing, and to support arbitrary extra dimensions such as time.

[mauro3]

how do we know which dimension is which?

There is https://github.com/JuliaArrays/AxisArrays.jl. (Edit: never mind, you're aware of it as stated in the issue you link)

As for the rest, probably best to start simple and add bits as needed.

[mkborregaard]

The issue about storing time is also something we're discussing a lot in the context of SpatialEcology.jl. WRT the index ordering, that seems like an aspect of implementation rather than interface? At least I always store my grid data in the long format (lats, longs, various values) because that plays more nicely with how we treat the data in analyses.

[visr]

issue about storing time

Good to know that the GDAL folks (Even) are also thinking about adding support for extra dimensions in the GDAL data model in RFC 75:
https://github.com/rouault/gdal/blob/rfc75_text/gdal/doc/source/development/rfc/rfc75_multidimensional_arrays.rst
OSGeo/gdal#1580