Re: [gdal-dev] GDAL raster processing library?

[Rutger]

Great discussion about improving the raster processing. As opposed to the
ideas above, of incorporating new machinery into GDAL, it might also
interesting to consider changing the API (a bit) so GDAL integrates better
with existing tools.

Within the Python ecosystem for example, there have recently emerged some
great tools for serious data processing. A combination of Dask (inc
Dask.distributed) and Numba for example allow relative easy scaling to
multicore (including GPU) and to clusters. 

Dask easily handles windowed functions as discussed above, but also much
more complicated relations between 'blocks' for example with linear algebra,
matrix multiplication etc.

Here is a small experiment i made a while back, to test if i could use GDAL
for the I/O to populate Dask arrays. It computes the slope, a windowed
function, from an elevation map:
http://nbviewer.jupyter.org/gist/RutgerK/a2ad8c074c78d000dd4a1e35cc229dee

You basically only need a Python '__getitem__' interface so a Dataset can be
indexed similar as a Numpy array, and some attributes like the shape of the
data. That's much easier to implement than handling and scheduling blocks
yourself.

The downside is of course that you rely on third partys, and you get a
different experience based on the programming language you use because the
availability of packages/modules will be different. My example uses Python
and therefore would be useless to users of other languages. 

Btw, has anyone thought of using vectorized Numba functions as GDAL pixel
functions? I have no clue whether that would technically work or even make
sense. Perhaps it would be a lot easier, and faster, then using Python
derived functions. 


Regards,
Rutger



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/GDAL-raster-processing-library-tp5275948p5280032.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDAL raster processing library?

[jramm]

Hi Ari,
I began some work to clarify my ideas here: www.github.com/JamesRamm/GeoAlg
perhaps there is potential for merging with your project? For
neighbourhoods I provide 2 iterators - one simple block based and one
'buffered'. The block can be the natural block or user defined. I also
planned to provide a 'mosaic' iterator to support the planned 'reduce'
function. For callbacks I initially planned a simple function, but then
realised that many algorithms may wish to save state, so I changed to a
simple interface class where there is a single function to override. My
final plan was to expose opencv and other functionality (eg maybe taudem)
as a library of ready made callbacks.

If you think it worthwhile and possible, I would be open to porting
relevant parts to your code. Features I see as absolutely essential are:
- 'map' and 'reduce' functions for raster datasets, which handle large data
and user defined callbacks
- ability to save state in callbacks. This is essential for more complex
algorithms.
- A rich set of block iterators

A ready made library of a number of callbacks (preferably by exposing
functionality from existing well known libraries) would also be useful and
help build interest I think.

I will be in Bonn at foss4g from the Tuesday evening, it would be good to
meet!

On 8 Aug 2016 1:56 p.m., "Ari Jolma-2 [via OSGeo.org]" <
ml-node+s1560n5279906...@n6.nabble.com> wrote:

> Hi,
>
> Sorry for missing this last month due to holidays.
>
> In fact in my approach I do not want to limit to single pixel operations
> but support neighborhood operations too. That is necessary for watershed
> analysis for example. Making sure that the neighborhood is readily
> available for computations is probably what is causing the biggest problem
> for the code. I also have a callback support. See for example
> https://github.com/ajolma/gdal/blob/trunk/gdal/map_algebra/test.cpp where
> on the lines 5 to 12 is a callback function.
>
> When I left the code for my extended holidays I was doing profiling and
> intended to test another approach for caching. Now the caching is based on
> the native GDAL blocks, which is efficient for going through a single
> dataset. Since the block may be different for different datasets this
> approach may not be the best when there are more than one dataset. For
> example the raster algebra in QGIS uses an approach where the block is the
> same for all datasets. That simplifies the code and may be very efficient
> as a whole.
>
> I'll be at the code sprint in FOSS4G 2016 and that would be a good place
> to discuss the RFC too.
>
> Best,
>
> Ari
>
> 13.07.2016, 10:09, James Ramm kirjoitti:
>
> Peter,
> I think 'Grid Algebra' would be what Ari Jolma is proposing here:
> https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra
>
> As Even pointed out, there is some overlap, though my proposal is
> technically very different.
> The key differences I see are:
>
> - Users can submit functions which operate on each sub window of the
> raster, rather than an algebraic expression. This potentially allows for
> much more complicated algorithms to be used (e.g. I dont think it would be
> possible to run a watershed segmentation with a raster algebra
> implementation, or to have algorithms which behave differently depending on
> the 'location' within the raster or the values of surrounding pixels etc
> etc).
>
> - Functions can be chained together for a complex processing toolchain.
> Some overlap with VRT here, although again this introduces a little more
> flexibility.
>
> On 12 July 2016 at 15:47, Peter Halls <[hidden email]
> > wrote:
>
>> James,
>>
>>   in reality, are you not requesting an implementation of Tomlin's
>> 'Grid Algebra' in GDAL?  That defines the whole range of functions from
>> whole raster to pixel and has the distinct advantage of being both
>> published and extremely well known because of other implementations ...
>> which also provide ready-made reference bases for the GDAL implementors ...
>>
>> Best wishes,
>> Peter
>>
>> On 12 July 2016 at 15:39, James Ramm <[hidden email]
>> > wrote:
>>
>>> Hi Even,
>>>
>>> The difference I see with pixel functions is that, as far as I
>>> understand, the pixel function is applied per pixel, so there is no
>>> possibility of e.g. the pixel buffer when have the function apply to
>>> 'blocks'.
>>> I may be way off, but many of the algorithms we deal with require some
>>> kind of neighbourhood search - a polygonise algorithm or flow direction
>>> algorithm being good examples.
>>> I dont think VRT pixel functions allow this?
>>>
>>> So in that sense I'd see a VRT being 'just' another potential input data
>>> source.
>>>
>>> Perhaps VRT pixel functions could be expanded to also allow 'window'
>>> functions?
>>>
>>> A downside is it requires creating a VRT even when you only want to
>>> apply a such a function to a single dataset

Re: [gdal-dev] GDAL raster processing library?

[Ari Jolma]

Hi,

Sorry for missing this last month due to holidays.

In fact in my approach I do not want to limit to single pixel operations 
but support neighborhood operations too. That is necessary for watershed 
analysis for example. Making sure that the neighborhood is readily 
available for computations is probably what is causing the biggest 
problem for the code. I also have a callback support. See for example 
https://github.com/ajolma/gdal/blob/trunk/gdal/map_algebra/test.cpp 
where on the lines 5 to 12 is a callback function.


When I left the code for my extended holidays I was doing profiling and 
intended to test another approach for caching. Now the caching is based 
on the native GDAL blocks, which is efficient for going through a single 
dataset. Since the block may be different for different datasets this 
approach may not be the best when there are more than one dataset. For 
example the raster algebra in QGIS uses an approach where the block is 
the same for all datasets. That simplifies the code and may be very 
efficient as a whole.


I'll be at the code sprint in FOSS4G 2016 and that would be a good place 
to discuss the RFC too.


Best,

Ari


13.07.2016, 10:09, James Ramm kirjoitti:

Peter,
I think 'Grid Algebra' would be what Ari Jolma is proposing here: 
https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra


As Even pointed out, there is some overlap, though my proposal is 
technically very different.

The key differences I see are:

- Users can submit functions which operate on each sub window of the 
raster, rather than an algebraic expression. This potentially allows 
for much more complicated algorithms to be used (e.g. I dont think it 
would be possible to run a watershed segmentation with a raster 
algebra implementation, or to have algorithms which behave differently 
depending on the 'location' within the raster or the values of 
surrounding pixels etc etc).


- Functions can be chained together for a complex processing 
toolchain. Some overlap with VRT here, although again this introduces 
a little more flexibility.


On 12 July 2016 at 15:47, Peter Halls > wrote:


James,

  in reality, are you not requesting an implementation of
Tomlin's 'Grid Algebra' in GDAL?  That defines the whole range of
functions from whole raster to pixel and has the distinct
advantage of being both published and extremely well known because
of other implementations ... which also provide ready-made
reference bases for the GDAL implementors ...

Best wishes,
Peter

On 12 July 2016 at 15:39, James Ramm mailto:jamessr...@gmail.com>> wrote:

Hi Even,

The difference I see with pixel functions is that, as far as I
understand, the pixel function is applied per pixel, so there
is no possibility of e.g. the pixel buffer when have the
function apply to 'blocks'.
I may be way off, but many of the algorithms we deal with
require some kind of neighbourhood search - a polygonise
algorithm or flow direction algorithm being good examples.
I dont think VRT pixel functions allow this?

So in that sense I'd see a VRT being 'just' another potential
input data source.

Perhaps VRT pixel functions could be expanded to also allow
'window' functions?

A downside is it requires creating a VRT even when you only
want to apply a such a function to a single dataset. Small
effort, but still a bit more than throwing in any GDALDataset
to be processed.

I see the overlap with raster algebra, although yes
technically very different.



On 12 July 2016 at 14:55, Even Rouault
mailto:even.roua...@spatialys.com>> wrote:

James,

There's some intersection with Ari's proposal :
https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra . At
least regarding the
overall purposes, since technically this is quite different.

Actually what you propose is very close to the existing
VRT pixel functions of
derived bands. See "Using Derived Bands" in
http://www.gdal.org/gdal_vrttut.html . In the last days,
we've merged
Antonio's work regarding a predefined set of pixel functions.
Perhaps some extension to allow passing user parameters to
the pixel func
could be useful. It is possible to use pixel functions
from Python as shown in

https://github.com/OSGeo/gdal/blob/trunk/autotest/gcore/testnonboundtoswig.py#L303
although this is a bit ugly as it uses ctypes and not
SWIG. But should be
possible through SWIG by introducing proper types
similarly to what is done
for the progress functions or error handler functions.

Even



_

Re: [gdal-dev] GDAL raster processing library?

[jramm]

I can see where some similarities with other new and existing GDAL work could
be a blocker on this, but I also think this adds a a degree more flexibility
allowing potentially any kind of complex processing to be carried out
without worrying/bothering about boilerplate. 

It would be good to find some kind of middle ground where this flexibility
can be added in without duplicating other efforts?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/GDAL-raster-processing-library-tp5275948p5276078.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDAL raster processing library?

[James Ramm]

Peter,
I think 'Grid Algebra' would be what Ari Jolma is proposing here:
https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra

As Even pointed out, there is some overlap, though my proposal is
technically very different.
The key differences I see are:

- Users can submit functions which operate on each sub window of the
raster, rather than an algebraic expression. This potentially allows for
much more complicated algorithms to be used (e.g. I dont think it would be
possible to run a watershed segmentation with a raster algebra
implementation, or to have algorithms which behave differently depending on
the 'location' within the raster or the values of surrounding pixels etc
etc).

- Functions can be chained together for a complex processing toolchain.
Some overlap with VRT here, although again this introduces a little more
flexibility.

On 12 July 2016 at 15:47, Peter Halls  wrote:

> James,
>
>   in reality, are you not requesting an implementation of Tomlin's
> 'Grid Algebra' in GDAL?  That defines the whole range of functions from
> whole raster to pixel and has the distinct advantage of being both
> published and extremely well known because of other implementations ...
> which also provide ready-made reference bases for the GDAL implementors ...
>
> Best wishes,
> Peter
>
> On 12 July 2016 at 15:39, James Ramm  wrote:
>
>> Hi Even,
>>
>> The difference I see with pixel functions is that, as far as I
>> understand, the pixel function is applied per pixel, so there is no
>> possibility of e.g. the pixel buffer when have the function apply to
>> 'blocks'.
>> I may be way off, but many of the algorithms we deal with require some
>> kind of neighbourhood search - a polygonise algorithm or flow direction
>> algorithm being good examples.
>> I dont think VRT pixel functions allow this?
>>
>> So in that sense I'd see a VRT being 'just' another potential input data
>> source.
>>
>> Perhaps VRT pixel functions could be expanded to also allow 'window'
>> functions?
>>
>> A downside is it requires creating a VRT even when you only want to apply
>> a such a function to a single dataset. Small effort, but still a bit more
>> than throwing in any GDALDataset to be processed.
>>
>> I see the overlap with raster algebra, although yes technically very
>> different.
>>
>>
>>
>> On 12 July 2016 at 14:55, Even Rouault 
>> wrote:
>>
>>> James,
>>>
>>> There's some intersection with Ari's proposal :
>>> https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra . At least
>>> regarding the
>>> overall purposes, since technically this is quite different.
>>>
>>> Actually what you propose is very close to the existing VRT pixel
>>> functions of
>>> derived bands. See "Using Derived Bands" in
>>> http://www.gdal.org/gdal_vrttut.html . In the last days, we've merged
>>> Antonio's work regarding a predefined set of pixel functions.
>>> Perhaps some extension to allow passing user parameters to the pixel func
>>> could be useful. It is possible to use pixel functions from Python as
>>> shown in
>>>
>>> https://github.com/OSGeo/gdal/blob/trunk/autotest/gcore/testnonboundtoswig.py#L303
>>> although this is a bit ugly as it uses ctypes and not SWIG. But should be
>>> possible through SWIG by introducing proper types similarly to what is
>>> done
>>> for the progress functions or error handler functions.
>>>
>>> Even
>>>
>>> Le mardi 12 juillet 2016 14:40:27, jramm a écrit :
>>> > I wonder if there is a use case/interest in a small library within
>>> GDAL to
>>> > facilitate generic raster processing.
>>> >
>>> > My idea would be to have a user-extensible framework to run processing
>>> > functions on whole rasters, with a growing library of common-such
>>> functions
>>> > within GDAL.
>>> >
>>> > Something along the lines of this:
>>> >
>>> > /***/
>>> > typedef CPLErr (*GDALRasterProcessFn)(double *padfInArray, double
>>> > *padfOutArray,
>>> > int nWindowXSize, int nWindowYSize, void *pData, double
>>> > *pdfNoDataValue);
>>> > /**
>>> > * \brief Definition of a raster processing function.
>>> > *
>>> > * A GDALRasterProcessFn accepts an array of data as input, applies
>>> custom
>>> > logic and writes the output to padfOutArray.
>>> > * Such a function can be passed to GDALRunRasterProcess to apply custom
>>> > processing to a GDALDataset in chunks and create
>>> > * a new GDALDataset.
>>> > *
>>> > * @param padfInArray The input array of data.
>>> > *
>>> > * @param padfOutArray The output array of data. On first call (via
>>> > GDALRunRasterProcess) this will be an empty, initialised array,
>>> > *which should be populated with the result of calculations on
>>> > padfInArray. In subsequent calls it will contain the result of the
>>> > *previous window.
>>> > *
>>> > * @param nWindowXSize the actual x size (width) of the read window.
>>> > *
>>> > * @param nWindowYSize the actual y size (height) of the read window.
>>> The
>>> > length of padf

Re: [gdal-dev] GDAL raster processing library?

[Even Rouault]

Le mardi 12 juillet 2016 16:39:34, James Ramm a écrit :
> Hi Even,
> 
> The difference I see with pixel functions is that, as far as I understand,
> the pixel function is applied per pixel, so there is no possibility of e.g.
> the pixel buffer when have the function apply to 'blocks'.
> I may be way off, but many of the algorithms we deal with require some kind
> of neighbourhood search - a polygonise algorithm or flow direction
> algorithm being good examples.
> I dont think VRT pixel functions allow this?

True, they don't. Although there's the KernelFilteredSource that is close

> 
> So in that sense I'd see a VRT being 'just' another potential input data
> source.
> 
> Perhaps VRT pixel functions could be expanded to also allow 'window'
> functions?

Would make sense

> 
> A downside is it requires creating a VRT even when you only want to apply a
> such a function to a single dataset. Small effort, but still a bit more
> than throwing in any GDALDataset to be processed.

You could have an helper function to create the VRT dataset from the input 
dataset, the processing function and its parameters. A bit similarly to what 
Julien is doing in 
https://github.com/OSGeo/gdal/pull/142/files#diff-0c9bf560508edf2453cd2f776d72f905R121

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDAL raster processing library?

[Peter Halls]

James,

  in reality, are you not requesting an implementation of Tomlin's
'Grid Algebra' in GDAL?  That defines the whole range of functions from
whole raster to pixel and has the distinct advantage of being both
published and extremely well known because of other implementations ...
which also provide ready-made reference bases for the GDAL implementors ...

Best wishes,
Peter

On 12 July 2016 at 15:39, James Ramm  wrote:

> Hi Even,
>
> The difference I see with pixel functions is that, as far as I understand,
> the pixel function is applied per pixel, so there is no possibility of e.g.
> the pixel buffer when have the function apply to 'blocks'.
> I may be way off, but many of the algorithms we deal with require some
> kind of neighbourhood search - a polygonise algorithm or flow direction
> algorithm being good examples.
> I dont think VRT pixel functions allow this?
>
> So in that sense I'd see a VRT being 'just' another potential input data
> source.
>
> Perhaps VRT pixel functions could be expanded to also allow 'window'
> functions?
>
> A downside is it requires creating a VRT even when you only want to apply
> a such a function to a single dataset. Small effort, but still a bit more
> than throwing in any GDALDataset to be processed.
>
> I see the overlap with raster algebra, although yes technically very
> different.
>
>
>
> On 12 July 2016 at 14:55, Even Rouault  wrote:
>
>> James,
>>
>> There's some intersection with Ari's proposal :
>> https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra . At least
>> regarding the
>> overall purposes, since technically this is quite different.
>>
>> Actually what you propose is very close to the existing VRT pixel
>> functions of
>> derived bands. See "Using Derived Bands" in
>> http://www.gdal.org/gdal_vrttut.html . In the last days, we've merged
>> Antonio's work regarding a predefined set of pixel functions.
>> Perhaps some extension to allow passing user parameters to the pixel func
>> could be useful. It is possible to use pixel functions from Python as
>> shown in
>>
>> https://github.com/OSGeo/gdal/blob/trunk/autotest/gcore/testnonboundtoswig.py#L303
>> although this is a bit ugly as it uses ctypes and not SWIG. But should be
>> possible through SWIG by introducing proper types similarly to what is
>> done
>> for the progress functions or error handler functions.
>>
>> Even
>>
>> Le mardi 12 juillet 2016 14:40:27, jramm a écrit :
>> > I wonder if there is a use case/interest in a small library within GDAL
>> to
>> > facilitate generic raster processing.
>> >
>> > My idea would be to have a user-extensible framework to run processing
>> > functions on whole rasters, with a growing library of common-such
>> functions
>> > within GDAL.
>> >
>> > Something along the lines of this:
>> >
>> > /***/
>> > typedef CPLErr (*GDALRasterProcessFn)(double *padfInArray, double
>> > *padfOutArray,
>> > int nWindowXSize, int nWindowYSize, void *pData, double
>> > *pdfNoDataValue);
>> > /**
>> > * \brief Definition of a raster processing function.
>> > *
>> > * A GDALRasterProcessFn accepts an array of data as input, applies
>> custom
>> > logic and writes the output to padfOutArray.
>> > * Such a function can be passed to GDALRunRasterProcess to apply custom
>> > processing to a GDALDataset in chunks and create
>> > * a new GDALDataset.
>> > *
>> > * @param padfInArray The input array of data.
>> > *
>> > * @param padfOutArray The output array of data. On first call (via
>> > GDALRunRasterProcess) this will be an empty, initialised array,
>> > *which should be populated with the result of calculations on
>> > padfInArray. In subsequent calls it will contain the result of the
>> > *previous window.
>> > *
>> > * @param nWindowXSize the actual x size (width) of the read window.
>> > *
>> > * @param nWindowYSize the actual y size (height) of the read window. The
>> > length of padfInArray == padfOutArray == nWindowXSize * nWindowYSize
>> > *
>> > * @param pData Process-specific data. This data is passed straight
>> through
>> > to the GDALRasterProcessFn and may contain e.g user defined parameters.
>> > * The GDALRasterProcessFn definition would define the
>> structure/type of
>> > such data.
>> > *
>> > * @param pdfNoDataValue The no data value of the dataset
>> > */
>> >
>> > CPLErr GDALRunRasterProcess(GDALRasterProcessFn processFn, GDALDataset
>> > *poSrcDataset,
>> > GDALDataset *poDstDataset, void *pData, int
>> > *pnWindowXSize,
>> > int *pnWindowYSize, int *pnPixelBuffer);
>> > /**
>> > * \brief Apply a raster processing function to each sub-window of a
>> raster.
>> > *
>> > * The input raster dataset is read in chunks of nWindowXSize *
>> nWindowYSize
>> > and each chunk is passed to the processing
>> > * function. The output array from the function is written to the
>> > destination dataset.
>> > * An optional 'pixel 

Re: [gdal-dev] GDAL raster processing library?

[James Ramm]

Hi Even,

The difference I see with pixel functions is that, as far as I understand,
the pixel function is applied per pixel, so there is no possibility of e.g.
the pixel buffer when have the function apply to 'blocks'.
I may be way off, but many of the algorithms we deal with require some kind
of neighbourhood search - a polygonise algorithm or flow direction
algorithm being good examples.
I dont think VRT pixel functions allow this?

So in that sense I'd see a VRT being 'just' another potential input data
source.

Perhaps VRT pixel functions could be expanded to also allow 'window'
functions?

A downside is it requires creating a VRT even when you only want to apply a
such a function to a single dataset. Small effort, but still a bit more
than throwing in any GDALDataset to be processed.

I see the overlap with raster algebra, although yes technically very
different.



On 12 July 2016 at 14:55, Even Rouault  wrote:

> James,
>
> There's some intersection with Ari's proposal :
> https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra . At least
> regarding the
> overall purposes, since technically this is quite different.
>
> Actually what you propose is very close to the existing VRT pixel
> functions of
> derived bands. See "Using Derived Bands" in
> http://www.gdal.org/gdal_vrttut.html . In the last days, we've merged
> Antonio's work regarding a predefined set of pixel functions.
> Perhaps some extension to allow passing user parameters to the pixel func
> could be useful. It is possible to use pixel functions from Python as
> shown in
>
> https://github.com/OSGeo/gdal/blob/trunk/autotest/gcore/testnonboundtoswig.py#L303
> although this is a bit ugly as it uses ctypes and not SWIG. But should be
> possible through SWIG by introducing proper types similarly to what is done
> for the progress functions or error handler functions.
>
> Even
>
> Le mardi 12 juillet 2016 14:40:27, jramm a écrit :
> > I wonder if there is a use case/interest in a small library within GDAL
> to
> > facilitate generic raster processing.
> >
> > My idea would be to have a user-extensible framework to run processing
> > functions on whole rasters, with a growing library of common-such
> functions
> > within GDAL.
> >
> > Something along the lines of this:
> >
> > /***/
> > typedef CPLErr (*GDALRasterProcessFn)(double *padfInArray, double
> > *padfOutArray,
> > int nWindowXSize, int nWindowYSize, void *pData, double
> > *pdfNoDataValue);
> > /**
> > * \brief Definition of a raster processing function.
> > *
> > * A GDALRasterProcessFn accepts an array of data as input, applies custom
> > logic and writes the output to padfOutArray.
> > * Such a function can be passed to GDALRunRasterProcess to apply custom
> > processing to a GDALDataset in chunks and create
> > * a new GDALDataset.
> > *
> > * @param padfInArray The input array of data.
> > *
> > * @param padfOutArray The output array of data. On first call (via
> > GDALRunRasterProcess) this will be an empty, initialised array,
> > *which should be populated with the result of calculations on
> > padfInArray. In subsequent calls it will contain the result of the
> > *previous window.
> > *
> > * @param nWindowXSize the actual x size (width) of the read window.
> > *
> > * @param nWindowYSize the actual y size (height) of the read window. The
> > length of padfInArray == padfOutArray == nWindowXSize * nWindowYSize
> > *
> > * @param pData Process-specific data. This data is passed straight
> through
> > to the GDALRasterProcessFn and may contain e.g user defined parameters.
> > * The GDALRasterProcessFn definition would define the structure/type
> of
> > such data.
> > *
> > * @param pdfNoDataValue The no data value of the dataset
> > */
> >
> > CPLErr GDALRunRasterProcess(GDALRasterProcessFn processFn, GDALDataset
> > *poSrcDataset,
> > GDALDataset *poDstDataset, void *pData, int
> > *pnWindowXSize,
> > int *pnWindowYSize, int *pnPixelBuffer);
> > /**
> > * \brief Apply a raster processing function to each sub-window of a
> raster.
> > *
> > * The input raster dataset is read in chunks of nWindowXSize *
> nWindowYSize
> > and each chunk is passed to the processing
> > * function. The output array from the function is written to the
> > destination dataset.
> > * An optional 'pixel buffer' can be specified to allow overlaps between
> > successive windows. This is useful for
> > * some algorithms, e.g. blob extraction, watershed/stream flow analysis,
> > convolution etc.
> > * Process specific data can be passed (e.g. configuration parameters).
> This
> > data is simply passed straight through to the processing
> > * function on each call.
> > *
> > * @param processFn A GDALRasterProcessFn to apply to each sub window of
> the
> > raster.
> > *
> > * @param poSrcDataset The source raster dataset from which pixel values
> are
> > read
> > *
> > 

Re: [gdal-dev] GDAL raster processing library?

[Even Rouault]

James,

There's some intersection with Ari's proposal :
https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra . At least regarding the 
overall purposes, since technically this is quite different.

Actually what you propose is very close to the existing VRT pixel functions of 
derived bands. See "Using Derived Bands" in 
http://www.gdal.org/gdal_vrttut.html . In the last days, we've merged 
Antonio's work regarding a predefined set of pixel functions.
Perhaps some extension to allow passing user parameters to the pixel func 
could be useful. It is possible to use pixel functions from Python as shown in 
https://github.com/OSGeo/gdal/blob/trunk/autotest/gcore/testnonboundtoswig.py#L303
 
although this is a bit ugly as it uses ctypes and not SWIG. But should be 
possible through SWIG by introducing proper types similarly to what is done 
for the progress functions or error handler functions.

Even

Le mardi 12 juillet 2016 14:40:27, jramm a écrit :
> I wonder if there is a use case/interest in a small library within GDAL to
> facilitate generic raster processing.
> 
> My idea would be to have a user-extensible framework to run processing
> functions on whole rasters, with a growing library of common-such functions
> within GDAL.
> 
> Something along the lines of this:
> 
> /***/
> typedef CPLErr (*GDALRasterProcessFn)(double *padfInArray, double
> *padfOutArray,
> int nWindowXSize, int nWindowYSize, void *pData, double
> *pdfNoDataValue);
> /**
> * \brief Definition of a raster processing function.
> *
> * A GDALRasterProcessFn accepts an array of data as input, applies custom
> logic and writes the output to padfOutArray.
> * Such a function can be passed to GDALRunRasterProcess to apply custom
> processing to a GDALDataset in chunks and create
> * a new GDALDataset.
> *
> * @param padfInArray The input array of data.
> *
> * @param padfOutArray The output array of data. On first call (via
> GDALRunRasterProcess) this will be an empty, initialised array,
> *which should be populated with the result of calculations on
> padfInArray. In subsequent calls it will contain the result of the
> *previous window.
> *
> * @param nWindowXSize the actual x size (width) of the read window.
> *
> * @param nWindowYSize the actual y size (height) of the read window. The
> length of padfInArray == padfOutArray == nWindowXSize * nWindowYSize
> *
> * @param pData Process-specific data. This data is passed straight through
> to the GDALRasterProcessFn and may contain e.g user defined parameters.
> * The GDALRasterProcessFn definition would define the structure/type of
> such data.
> *
> * @param pdfNoDataValue The no data value of the dataset
> */
> 
> CPLErr GDALRunRasterProcess(GDALRasterProcessFn processFn, GDALDataset
> *poSrcDataset,
> GDALDataset *poDstDataset, void *pData, int
> *pnWindowXSize,
> int *pnWindowYSize, int *pnPixelBuffer);
> /**
> * \brief Apply a raster processing function to each sub-window of a raster.
> *
> * The input raster dataset is read in chunks of nWindowXSize * nWindowYSize
> and each chunk is passed to the processing
> * function. The output array from the function is written to the
> destination dataset.
> * An optional 'pixel buffer' can be specified to allow overlaps between
> successive windows. This is useful for
> * some algorithms, e.g. blob extraction, watershed/stream flow analysis,
> convolution etc.
> * Process specific data can be passed (e.g. configuration parameters). This
> data is simply passed straight through to the processing
> * function on each call.
> *
> * @param processFn A GDALRasterProcessFn to apply to each sub window of the
> raster.
> *
> * @param poSrcDataset The source raster dataset from which pixel values are
> read
> *
> * @param poDstDataset The destination raster dataset to which pixel values
> are written. Must support RasterIO in write mode.
> *
> * @param pData Process-specific data. This is passed straight through to
> the GDALRasterProcessFn on each call.
> *
> * @param pnWindowXSize The desired width of each read window. If NULL it
> defaults to the 'natural' block size of the raster
> *
> * @param pnWindowYSize The desired height of each read window. If NULL it
> defaults to the 'natural' block size.
> *
> * @param pnPixelBuffer A pixel buffer to apply to the read window. The read
> window is expanded by pnPixelBuffer pixels in all directions such that
> *each window overlaps by pnPixelBuffer pixels. Ignored if NULL or 0
> *
> * @return a CPLErr enum indicating whether the function succesfully ran.
> */
> 
> 
> CPLErr GDALRunRasterProcesses(GDALRasterProcessFn *paProcessFn, int
> nProcesses, GDALDataset *poSrcDataset,
> GDALDataset *poDstDataset, void **paData, int
> *pnWindowXSize,
> int *pnWindowYSize, int *pnPixelBuffer);
> 
> /**
> * \brief Apply multiple ras

[gdal-dev] GDAL raster processing library?

[jramm]

I wonder if there is a use case/interest in a small library within GDAL to
facilitate generic raster processing. 

My idea would be to have a user-extensible framework to run processing
functions on whole rasters, with a growing library of common-such functions
within GDAL. 

Something along the lines of this:

/***/
typedef CPLErr (*GDALRasterProcessFn)(double *padfInArray, double
*padfOutArray, 
int nWindowXSize, int nWindowYSize, void *pData, double
*pdfNoDataValue);
/**
* \brief Definition of a raster processing function. 
*
* A GDALRasterProcessFn accepts an array of data as input, applies custom
logic and writes the output to padfOutArray.
* Such a function can be passed to GDALRunRasterProcess to apply custom
processing to a GDALDataset in chunks and create 
* a new GDALDataset. 
*
* @param padfInArray The input array of data. 
*
* @param padfOutArray The output array of data. On first call (via
GDALRunRasterProcess) this will be an empty, initialised array,
*which should be populated with the result of calculations on
padfInArray. In subsequent calls it will contain the result of the 
*previous window. 
*
* @param nWindowXSize the actual x size (width) of the read window.
*
* @param nWindowYSize the actual y size (height) of the read window. The
length of padfInArray == padfOutArray == nWindowXSize * nWindowYSize
*
* @param pData Process-specific data. This data is passed straight through
to the GDALRasterProcessFn and may contain e.g user defined parameters.
* The GDALRasterProcessFn definition would define the structure/type of
such data.
*
* @param pdfNoDataValue The no data value of the dataset
*/

CPLErr GDALRunRasterProcess(GDALRasterProcessFn processFn, GDALDataset
*poSrcDataset, 
GDALDataset *poDstDataset, void *pData, int
*pnWindowXSize, 
int *pnWindowYSize, int *pnPixelBuffer);
/**
* \brief Apply a raster processing function to each sub-window of a raster.
* 
* The input raster dataset is read in chunks of nWindowXSize * nWindowYSize
and each chunk is passed to the processing
* function. The output array from the function is written to the destination
dataset. 
* An optional 'pixel buffer' can be specified to allow overlaps between
successive windows. This is useful for 
* some algorithms, e.g. blob extraction, watershed/stream flow analysis,
convolution etc.
* Process specific data can be passed (e.g. configuration parameters). This
data is simply passed straight through to the processing
* function on each call. 
*
* @param processFn A GDALRasterProcessFn to apply to each sub window of the
raster. 
* 
* @param poSrcDataset The source raster dataset from which pixel values are
read
*
* @param poDstDataset The destination raster dataset to which pixel values
are written. Must support RasterIO in write mode.
*
* @param pData Process-specific data. This is passed straight through to the
GDALRasterProcessFn on each call. 
*
* @param pnWindowXSize The desired width of each read window. If NULL it
defaults to the 'natural' block size of the raster
*
* @param pnWindowYSize The desired height of each read window. If NULL it
defaults to the 'natural' block size.
*
* @param pnPixelBuffer A pixel buffer to apply to the read window. The read
window is expanded by pnPixelBuffer pixels in all directions such that 
*each window overlaps by pnPixelBuffer pixels. Ignored if NULL or 0
*
* @return a CPLErr enum indicating whether the function succesfully ran. 
*/


CPLErr GDALRunRasterProcesses(GDALRasterProcessFn *paProcessFn, int
nProcesses, GDALDataset *poSrcDataset, 
GDALDataset *poDstDataset, void **paData, int
*pnWindowXSize, 
int *pnWindowYSize, int *pnPixelBuffer);

/**
* \brief Apply multiple raster processing functions to each sub-window of a
raster
*
* For each window, the functions defined by the paProcessFn array are called
in turn, with the array output of the previous function forming the input
* to the next function. This allows processing 'toolchains' to be built
without having to create intermediate datasets, which can be less efficient
in time and space. 
*
*
* @param paProcessFn An array of GDALRasterProcessFn to apply to each sub
window of the raster
*
* @param nProcesses The size of paProcessFn
*
* @param poSrcDataset The source raster dataset from which pixel values are
read
*
* @param poDstDataset The destination raster dataset to which pixel values
are written. Must support RasterIO in write mode.
*
* @param paData an array of process-specific data objects of size
nProcesses. Each data object will be passed to the corresponding
GDALRasterProcessFn
*
* @param pnWindowXSize The desired width of each read window. If NULL it
defaults to the 'natural' block size of the raster
*
* @param pnWindowYSize The desired height of each read window. If NULL it