[gdal-dev] Map algebra

[alex]

> -----Original Message-----
> From: gdal-dev [mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of
> Gregory, Matthew
> Sent: 08 April 2016 17:29
> To: gdal-dev at lists.osgeo.org
> Subject: Re: [gdal-dev] Map algebra
> 
> Hi all,
> 
> Alex wrote:
> > What exactly do you mean by Map Algebra? I understand it to cover the
> > following:
> >
> > 1. Local operations: pixel-by-pixel operations on raster bands.
> > 2. Focal operations: operations that calculate an indicator for each
> > pixel based on the values in surrounding pixels 3. Zonal operations:
> > operations that calculate indicators for groups
> > (zones) of pixels.
> > 4. Global operations: operations that calculate indicators over all
> > pixels.
> 
> I've been following this discussion with limited understanding of the C++ needed
> to implement such a library, but I echo the functionality that Alex lists above.  I
> *would* like to add my two cents in that 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)
> 

My preferred approach for nodata values is to wrap them as boost::optional<T> (or forthcoming std:: optional<T>). One nice thing about representing as ranges is that operations on ranges can be chained. It is quite trivial to transform a range of T's to a range of boost:::optional<T>'s, e.g. as follows: 

auto band_checked = ma::apply(nodata_checker<int>{nodata_value}, band);

where: 

template<class T>
struct nodata_checker
{
  nodata_checker(const T& nodata) : m_nodata(nodata){}

  boost::optional<T> operator()(const T& value)
  {
    if(value != m_nodata) return value; // else return boost::none;
  }
  T m_nodata;
};

To complement this, the function object that is working on the values in the bands can be wrapped in another function object that intercepts nodata values. 

e.g. in simplified form for a single argument:

boost::optional<f_return_type> f_checked(Function&& f, boost::optional<T>& arg)
{
	if(arg) return f(arg);  // else return boost::none
} 

So even though I have not tried it, I would think that supporting nodata values in range-based map algebra is fairly simple to implement.

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

Again, I haven't actually done it, but I think working with ranges makes it relatively easy to support this. The framework requires ranges to have the same extent and cell-size.  It would be pretty straightforward to define a range that iterates over a given extent that is different from the raster containing it (especially if cell centres align). From then on, all the same methods can be applied. The Adobe Generic Image Library uses a similar approach.

It will be more tricky when the ranges have a different orientation, but still just a linear transformation. 

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

Different cell-sizes, like different extents should be supportable through a range transform. 

> So if raster 'a' was 1000x1000 30m in UTM and raster 'b' was 500x500 90m in
> WGS84, one might be able to do:
> 
>   c = a + b
> 
> with reasonable defaults chosen for output 'c' as well as the option to set output
> options explicitly (I fully admit the reprojection on-the-fly is a big ask ...).
>

I wouldn't know how to do this. But somebody more versed in on-the-fly reproject might just...

 
> Sorry if this hijacks the discussion, but I appreciate your efforts to envision such
> a system.
> 

I don't think you are hijacking the discussion. If GDAL is to support Map Algebra it is necessary to discuss what kind of expectations we have.

> thanks, matt
>