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
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
