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@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Spatialys - Geospatial professional services http://www.spatialys.com _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
