[gdal-dev] Map algebra

[alex]

>-----Original Message-----
>From: gdal-dev [mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of Ari
>Jolma
>Sent: 10 April 2016 09:20
>To: gdal-dev at lists.osgeo.org
>Subject: Re: [gdal-dev] Map algebra
>
>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;
>}
>

OK, thanks I get the general picture now. 


>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).
>

In the end, all you need to do that is type specific is to apply the function on the pixel (at least in many cases). Would it be an option to only make the pixel level function a template, rather than the block level function? 

>Best,
>
>Ari
>
>_______________________________________________
>gdal-dev mailing list
>gdal-dev at lists.osgeo.org
>http://lists.osgeo.org/mailman/listinfo/gdal-dev