[gdal-dev] Map algebra

[Ari Jolma]

08.04.2016, 19:28, Gregory, Matthew kirjoitti:
> I would want any map algebra library to handle three types of issues when dealing with raster data:
>
>    1) Handling mask or nodata values (e.g. an analysis mask should be
>       able to be set on any algebra operation that retains or combines
>       nodata values from its component operands)

GDAL supports nodata values and mask bands to identify cells that do not 
contain valid data. In arithmetics I think it is straightforward that  x 
+ nodata = nodata. In logical operations it depends on the 
interpretation of nodata. If it means unknown, then we are dealing with 
three-valued logics where for example unknown OR true = true. Think 
about for example a raster which shows polluted areas as a result of two 
studies (true, false, unknown). If one study shows a site polluted and 
the other has no data, then the result should be that the site is polluted.

I haven't yet got to the point in my code where I look at nodata and 
mask bands.

>
>    2) Handling different spatial extents (e.g. the ability to specify
>       an output extent -- coming from a specific raster or from the
>       union/intersection of rasters)

This may be desirable in some cases. For example if we have a huge 
raster and want to study some parts of it. In general I would, however, 
leave this for other already existing GDAL tools, especially 
gdal_translate, whose functionality is available in APIs as well.

>
>    3) Handling different cell sizes (e.g. the ability to specify an
>       analysis cell size and resample/reproject(??) all operands based
>       on this parameter)

I think using gdal_translate in this case as a first step makes even 
more sense.

Alex asked about how I iterate over blocks when there are two or more 
bands involved. The basic idea is that all methods operate on one band, 
and other bands (in my case if there is a second band) are adjusted to 
that. The algorithm is the following. "band" is a wrapper (struct) for a 
GDALRasterBand object and block cache. A block cache is simply an array 
of pointers to data, which are obtained from the GDALRasterBand object 
with ReadBlock() and which may be written back with WriteBlock().

band1 = gma_band_initialize(b1);
band2 = gma_band_initialize(b2);

while (iterate) {

         block_index i;
         for (i = blocks in band1) {

                 add block i to cache in band1;
                 block1 = get block i from band1;

                 update cache in band1 to allow desired focal distance;
                 update cache in band2 to allow desired focal distance;

                 call method specific callback with (band1, block1, 
band2, retval, arg);
                 make a note if the callback indicates a need for iteration;

                 based on the return value from the callback

                     return with error,
                     continue with next block from band1, or
                     write block1 to band1 and continue with next block 
from band1;

             }
         }
         if iteration, then in some cases some band level things need to 
be done;
     }
     empty cache in band1;
     empty cache in band2;
}

The cache update function reads all needed blocks into the cache and 
discards those that are not needed. If there were no discards, the cache 
would eventually contain the whole band.

In the method callback values from the band2 are obtained by first 
computing the global cell index and then the respective block and cell 
index within the block.

This may not be the best possible algorithm but to me it is pretty 
understandable, which is a goal in itself.

The code is at

https://github.com/ajolma/gdal/tree/trunk/gdal/map_algebra

The curse of dimensionality comes in at the point, where the above 
function is called. The method callback is a template function and it 
needs the datatype of each band that is given to it as an argument. (I 
see now that I have made also the above function a template function, 
which may not be needed).

Best,

Ari