[gdal-dev] Map algebra

[Ari Jolma]

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.

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

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
>