[gdal-dev] Map algebra

[Even Rouault]

Le vendredi 01 avril 2016 10:38:48, Ari Jolma a écrit :
> I've got this into a rather nice working shape. The code needs to be
> reorganized (it is now, horror, all in header files) but that's another
> thing.

From the perspective of people wanting to use the core mechanisms to extend 
with their own operations, that could be interesting I guess.

> 
> I've so far built this around a basic DEM catchment analysis workflow.
> It includes many kinds of steps so it is a pretty good benchmark. In
> pseudo code the workflow is something like
> 
> dem = open_from_file
> dem->fill_depressions
> fd = dem->compute_flow_directions
> fd->route_flat_cells
> ua = fd->compute_upstream_area
> c = new_raster // compute the catchment into this
> c->set_border_cells // we're interested in the outlet cells
> c *= ua
> c *= c > 10000 // we're interested in the big catchments
> outlets = c->get_cells // outlets is an array of cells, and a cell is an
> object containing (x,y,value)
> c = 0
> c->mark_catchment(fd, outlets[0], 1) // assuming the outlet cell is now
> in outlets[0]
> ua2 = new_raster // for visualization log of upstream area is nice
> ua2 += ua
> ua2->log

A few thoughts and questions :
- I guess what would be cool from the perspective of a C++ user is to be 
really to write the above pseudo code with C++ operator overload when 
appropriate !  A += B, A.min(B), min_A= A.min(), etc...
- There might be some numpy reinvention here... Or at least, it could perhaps 
serve as an inspiration. I guess Python users would in priority use numpy for 
the most basic operations. Of course the GIS specific methods would still be of 
value. 
- Saturation arithmetics could be interesting. numpy by default do wrapped 
arithmetics, you have to cast to a larger type and clip down afterwards, 
whereas OpenCV does it by default : 
http://opencvpython.blogspot.fr/2012/06/difference-between-matrix-arithmetic-
in.html
- Are the methods handling 2 bands able to deal with different block sizes ?
- And with bands of different data types ?
- Why the need to reinvent a hashset mechanism ? Isn't std::set good enough ?
- Same for gma_array vs std::vector ?
- various methods take a void* argument, whereas it seems that could be 
rewritten as datatype* for better type safety
- there seems to be a block cache mechanism. What is it for ? How does that 
work ?
- there is some overlap/connection with the pixel functions of VRT ( see 
"Using Derived Bands" of http://gdal.org/gdal_vrttut.html).


> 
> The C code is of course not this nice but it is not very far from this.
> All the basic four method types I mention below are implemented and
> adding new methods is not very difficult. The code makes heavy use of
> templates and macros in an attempt for versatility (GDAL rasters can
> have many datatypes) and easiness to write methods in a few lines of
> code. The whole set of code is now ~2500 lines but it creates rather
> large executables due to templates. The code is C++ but the API is C
> like, I've so far done no changes to existing GDAL code because of this.
> 
> As written below, the code works on two levels of x,y loops: the blocks
> within a band, which is generic and cells within a block, which is
> separate for each method. The methods work on line by line (no recursion
> etc.) and thus for example the fill_depressions is not optimal but
> instead very simple. Also, very large rasters can be processed.
> 
> I think this may be somewhat similar thing in relation to GDAL as the
> network support, so it could be a thing to add to the library if wished.
> 
> I've added RFC 62 "raster algebra" to spark discussion.
> 
> Ari
> 
> 20.03.2016, 10:21, Ari Jolma kirjoitti:
> > Folks,
> > 
> > I played around a bit with the map algebra idea and came up with
> > something that could work nicely as a plugin. The initial code is at
> > https://github.com/ajolma/gdal/tree/trunk/gdal/map_algebra
> > 
> > The idea is based on processing raster bands block by block and
> > calling a given function on each block. Using macros and templates the
> > code should be rather nice and manageable even when supporting
> > multiple data types for raster cells.
> > 
> > Further, the idea is to support three kinds of map algebra methods,
> > mostly working on a raster band in-place.
> > 
> > 1) Simple ones, which do not take any arguments beyond the raster band
> > itself, examples are functions like log, sin, rand.
> > 
> > 2) Those which take one non-spatial argument. The argument can be a
> > scalar like a number, or more complex like a reclassifier.
> > 
> > 3) Those which operate on two rasters. Examples are summation of
> > rasters and computing flow directions on a DEM. The latter is a bit
> > more complex since it is a focal method and requires a 3 x 3 matrix of
> > blocks from the other operand raster.
> > 
> > Maybe a fourth kind of methods are those which compute something from
> > a raster.
> > 
> > This would lead to four functions in the C API and raster band methods
> > in the bindings.
> > 
> > Comments are welcome, the initial code base contains a small demo
> > program. Eventually, if this works, I'll make a RFC from this.
> > 
> > Ari
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com