[gdal-dev] Map algebra

[alex]

> -----Original Message-----
> From: gdal-dev [mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of 
> Ari Jolma
> Sent: 08 April 2016 07:05
> To: gdal-dev at lists.osgeo.org
> Subject: Re: [gdal-dev] Map algebra
> 
> 08.04.2016, 02:21, alex kirjoitti:
> >> That is one aspect, yes. However, so far I've considered only 
> >> operations on one or two bands at time. Operations on more than two 
> >> can of course be dissolved into a series of operations with max two at a time.
> > That will often be possible, but not always. For instance how would 
> > you dissolve the following function? (without resorting to complete 
> > ugliness like c * a  + !c * b)
> >
> > int conditional(bool c, int a, int b) { return c ? a : b; }
> 
> As an assignment that would be
> 
> d = c ? a : b;
> 
> If all variables are spatial then that's four bands. To save space 
> we'll do in-place modifications to bands. If we do not want to change 
> any of the bands a, b and c, we'll begin by making d a copy of b:
> 
> d = b;
> 
> and then use a conditional
> 
> d = a if c;
> 

Ha, that was easier than I thought :)

 But also quite inefficient. And you'll still need to negotiate three different block sizes.

[snip]

 

> >> The approach is also strongly supported by the basic design of GDAL 
> >> and the fundamental algorithm is presented in GDAL documentation (I 
> >> gave the link in earlier email).
> > I didn't see that email and just trawled through the whole thread, 
> > can you give
> it again?
> 
> http://www.gdal.org/classGDALRasterBand.html#ad80cecc562fd48f5783677de
> 625360ac
> 
> it is a simple algorithm but demonstrates the two level approach.

OK, I also use the ReadBlock and WriteBlock functions instead of the RasterIO. But I still don't see how the principle of iterating over blocks and within blocks is an approach that can be used for two or three bands with different block sizes. 

> >
> >> The solution I already have seems to work well for the one and two 
> >> band case even with 3x3 focal region and as far as I can see it 
> >> could be extended to cases involving more bands.
> 
> I'll add the disclaimer about the curse of dimensionality I mention 
> above for which a solution should be found first. It could be a simple 
> one like copy all blocks to same datatype first.

Can you describe the principle of your method? And where your curse of dimensionality comes in. Perhaps somebody has a solution.

[snip]
> 
> It really depends on the block size, which is something GDAL core and 
> drivers handle.

> >>> 3. Zonal operations: operations that calculate indicators for 
> >>> groups
> >>> (zones) of
> >> pixels.
> >>> This seems more trivial than the other two components of map 
> >>> algebra and I
> >> haven't given it much attention, but I think should be supported.
> >>
> >> One complication, which I just today worked on, seems to stem from 
> >> one requirement that I have. I want the API be free of templates 
> >> and thus there is a need for generic classes for, e.g., numbers 
> >> (GDAL supports bands where data is various length integers or floats).
> >> Thus, for example a simple zonal method, which returns which zones 
> >> share a boundary returns a rather complex data structure.
> > I understand zonal methods to be operating on pixels within zones 
> > and not on
> the relationship between zones. What kind of methods do you have in mind?
> 
> I admit that it is perhaps not a "pure" zonal method. For example, 
> think about the question: "What are the neighboring countries of each country?".

Sure, that is useful functionality, but I was trying to find out what is the rationale and purpose of your Map Algebra library. I think I misunderstood your original announcement email. I thought you intended to offer generic support for Map Algebra type analysis. Instead it seems that you are offering a specific set of Map Algebra functions (as I can now see in the RFC).