13.09.2016, 15:40, alex kirjoitti:
Ari Wrote:
At FOSS4G my conclusion was that my first approach was not good for many reasons,
most >importantly because it did not scale to several bands in one operation. So
I've trying with the >following ideas/decisions:
The problem is to compute y=f(x1, x2, ...), where y is a new dataset with one band
or an existing >dataset, into which a new band is added.
x1, x2, ... are existing bands. f is an expression. The goal is to be able to
create an expression >object, with which one can write
Hi Ari,
Perhaps it is worthwhile looking back at the email I wrote earlier in the
process:
http://permalink.gmane.org/gmane.comp.gis.gdal.devel/42609
At the time, I suggested using expression objects and to allow functions with
more than 2 arguments and gave some suggestions on the implementation.
Yes, that has had impact on my thinking.
In my experience the key first step is to develop pixel iterators for raster
bands. Once you have that, raster bands can be wrapped as regular ranges and
you can make use of 'regular' tools and methods to develop your expression
objects, e.g. as in the Range v3 library by Eric Niebler.
Hence, my repeated suggestion for GDAL: create iterators for raster bands :)
It is still my opinion that the iterators must be in practice two level:
block and cell.
Anyway, I might at least cease working on the RFC now, if not withdraw
it. I'll rewrite the RFC text to reflect what I found out.
The main reason is that I took a new look at using PDL (Perl Data
Language) for raster algebra at the block level. I've had some support
for it in the Perl bindings for a long time but improved it lately. I
haven't done a comprehensive review but basically PDL is for Perl
perhaps what NumPy is for Python. I've now improved and documented the
'Piddle' method for Band objects, with which one can read data from Band
to a PDL object and vice versa. The method is surprisingly simple and
doesn't need any support from Swig. That solves mostly my immediate
needs but of course leaves open for example the question how to link
GDAL to for example some terrain analysis algorithms which need
neighborhood operations etc.
For information about PDL go to http://pdl.perl.org/
PDL and some related modules support parallel programming which is
interesting but I haven't yet any experience from it.
Ari
As you may recall, the code below works just fine (but requires C++11).
With kind regards, Alex
#include <blink/map_algebra/map_algebra.h>
namespace ma = blink::map_algebra;
int my_function(int w, int x, int y, int z)
{
return w * x + y * z;
}
int main()
{
auto a = ma::open_read_only<int>("input_1.tif");
auto b = ma::open_read_only<int>("input_2.tif");
auto w = ma::create_from_model<double>("output_1.tif", a);
auto x = ma::create_from_model<double>("output_2.tif", a);
auto y = ma::create_from_model<double>("output_3.tif", a);
auto z = ma::create_from_model<double>("output_4.tif", a);
// Example 1: Using operators
w = a + 3 * b;
// Example 2: Using assigning operators
x = 1;
x *= a;
x += 3 * b;
// Example 3: Map algebra using cell-by-cell functions
y = ma::apply(my_function, 1, a, 3, b);
// Example 4: Combination
z = b + 3 * ma::apply(my_function, 1, a, 2, b);
return 0;
}
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
