Add raster intersection functionality

[jdhughes-usgs]

from @zstanko-usgs:

I've been working with some functions to read in vector/raster files and resample them to a model grid. similar to what Jeff Starn did here with the geospatial processing notebook.

I'm wondering if there is any interest in being able to fold some of that into Flopy utilities. The dependency on any geospatial libraries is certainly a concern but I see the export(.tif) requires rasterio so an input(*.tif) feature might not be that much of a stretch.

Basically, for any dataset that could could be entered as an array, it would be sweet if we could import it from a geotiff, ascii grid, or shapefile. What would be the limitations for being able to do this within Flopy code or do you think this would always best be done external to Flopy?

[jlarsen-usgs]

@zstanko-usgs:

I think that having some raster intersection support would be pretty useful especially for working with DEM raster data. It seems that vector file operations are being addressed with issue #536

@jdhughes-usgs has there been any work on the raster intersection issue? I'd be willing to put together a draft object that applies sampling/resampling methods for raster objects.

[jdhughes-usgs]

@jlarsen-usgs Nothing has been done yet. Feel free to put together a draft object.

[jonathanqv]

It´ll be awesome if It could be implemented in Modflow. I´ve worked with some structured grids using scipy.interpolate.griddata and the DEM is processed as an .xyz file, here is my script. It worked pretty well for me.

from scipy.interpolate import griddata
dem = pd.read_csv('../RST/asscigrid.xyz', sep=' ',header=None)
points = dem[[0,1]].values
values = dem[2].values
grid_x = mf.modelgrid.xcellcenters
grid_y = mf.modelgrid.ycellcenters
mtop = griddata(points, values, (grid_x, grid_y), method='nearest')

[jlarsen-usgs]

I've put together a draft class for raster intersection and resampling. Here's a link to an example I put together that shows the class functionality. The source code is located in the same repository.

rasters example

[langevin-usgs]

I like this @jlarsen-usgs. Very sleek and usable. Would be good to hear from @zstanko-usgs and @dbrakenhoff. Any thoughts?

[dbrakenhoff]

This looks really nice! Most of the operations I can think of are included. I might have a project soon in which I'd definitely use/test this.

It might be cool to have a function that allows you to calculate some statistic within a model gridcell based on the rastercells within that model gridcell. Sort of a repeated crop where each model gridcell is a polygon. In that case having the option of not cropping in-place might be useful. Though there might also be much better ways to go about something like that.

[jlarsen-usgs]

@dbrakenhoff, I think I understand what you mean by the statistics function. It would be possible to add a method that given a polygon will return an unordered list or generator object which contains all raster values within that polygon. Alternatively, it could return a small object that has built in statistics methods attached to it. What are your thoughts on that?

[dbrakenhoff]

@jlarsen-usgs, yea I think either of those options would work. I could live with either implementation but I think an unordered list is probably simplest(?) and gives the user the freedom to apply their own statistics. Plus it saves having to create those statistics methods for a new object? On the other hand, some users might find it nice to have an object like that ready at their fingertips without having to think about it. But my preference would be to keep it simple.

[jlarsen-usgs]

@dbrakenhoff, I added a sample_polygon method to the Rasters class this morning. This method returns an unordered list based on the intersection of a user supplied polygon and a raster. I also pushed an update to the notebook example to include this new method.

rasters example

[langevin-usgs]

@jlarsen-usgs what do we need to do to get this into FloPy? Is that where this should live? Would you say that it is developed enough to incorporate?

[jlarsen-usgs]

@langevin-usgs this should be relatively easy to incorporate into FloPy. Since Rasters is a stand alone class it could easily live in the utils folders. I think that the base functionality is developed enough to be included in FloPy. Additional methods could be added on later if there is a need for more advanced features.