<div dir="ltr">Peter,<div>I think 'Grid Algebra' would be what Ari Jolma is proposing here: <a href="https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra" rel="noreferrer" target="_blank" style="font-size:12.8px">https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra</a></div><div><br></div><div>As Even pointed out, there is some overlap, though my proposal is technically very different. </div><div>The key differences I see are:</div><div><br></div><div>- 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).</div><div><br></div><div>- Functions can be chained together for a complex processing toolchain. Some overlap with VRT here, although again this introduces a little more flexibility. </div></div><div class="gmail_extra"><br><div class="gmail_quote">On 12 July 2016 at 15:47, Peter Halls <span dir="ltr"><<a href="mailto:p.halls@york.ac.uk" target="_blank">p.halls@york.ac.uk</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div>James,<br><br></div> 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 ...<br><br></div>Best wishes,<br>Peter<br></div><div class="gmail_extra"><div><div class="h5"><br><div class="gmail_quote">On 12 July 2016 at 15:39, James Ramm <span dir="ltr"><<a href="mailto:jamessramm@gmail.com" target="_blank">jamessramm@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi Even,<div><br></div><div>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'. </div><div>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. </div><div>I dont think VRT pixel functions allow this? </div><div><br></div><div>So in that sense I'd see a VRT being 'just' another potential input data source. </div><div><br></div><div>Perhaps VRT pixel functions could be expanded to also allow 'window' functions?</div><div><br></div><div>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. </div><div><br></div><div>I see the overlap with raster algebra, although yes technically very different. </div><div><br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On 12 July 2016 at 14:55, Even Rouault <span dir="ltr"><<a href="mailto:even.rouault@spatialys.com" target="_blank">even.rouault@spatialys.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">James,<br>
<br>
There's some intersection with Ari's proposal :<br>
<a href="https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra" rel="noreferrer" target="_blank">https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra</a> . At least regarding the<br>
overall purposes, since technically this is quite different.<br>
<br>
Actually what you propose is very close to the existing VRT pixel functions of<br>
derived bands. See "Using Derived Bands" in<br>
<a href="http://www.gdal.org/gdal_vrttut.html" rel="noreferrer" target="_blank">http://www.gdal.org/gdal_vrttut.html</a> . In the last days, we've merged<br>
Antonio's work regarding a predefined set of pixel functions.<br>
Perhaps some extension to allow passing user parameters to the pixel func<br>
could be useful. It is possible to use pixel functions from Python as shown in<br>
<a href="https://github.com/OSGeo/gdal/blob/trunk/autotest/gcore/testnonboundtoswig.py#L303" rel="noreferrer" target="_blank">https://github.com/OSGeo/gdal/blob/trunk/autotest/gcore/testnonboundtoswig.py#L303</a><br>
although this is a bit ugly as it uses ctypes and not SWIG. But should be<br>
possible through SWIG by introducing proper types similarly to what is done<br>
for the progress functions or error handler functions.<br>
<br>
Even<br>
<div><div><br>
Le mardi 12 juillet 2016 14:40:27, jramm a écrit :<br>
> I wonder if there is a use case/interest in a small library within GDAL to<br>
> facilitate generic raster processing.<br>
><br>
> My idea would be to have a user-extensible framework to run processing<br>
> functions on whole rasters, with a growing library of common-such functions<br>
> within GDAL.<br>
><br>
> Something along the lines of this:<br>
><br>
> /*******************************************************************/<br>
> typedef CPLErr (*GDALRasterProcessFn)(double *padfInArray, double<br>
> *padfOutArray,<br>
> int nWindowXSize, int nWindowYSize, void *pData, double<br>
> *pdfNoDataValue);<br>
> /**<br>
> * \brief Definition of a raster processing function.<br>
> *<br>
> * A GDALRasterProcessFn accepts an array of data as input, applies custom<br>
> logic and writes the output to padfOutArray.<br>
> * Such a function can be passed to GDALRunRasterProcess to apply custom<br>
> processing to a GDALDataset in chunks and create<br>
> * a new GDALDataset.<br>
> *<br>
> * @param padfInArray The input array of data.<br>
> *<br>
> * @param padfOutArray The output array of data. On first call (via<br>
> GDALRunRasterProcess) this will be an empty, initialised array,<br>
> * which should be populated with the result of calculations on<br>
> padfInArray. In subsequent calls it will contain the result of the<br>
> * previous window.<br>
> *<br>
> * @param nWindowXSize the actual x size (width) of the read window.<br>
> *<br>
> * @param nWindowYSize the actual y size (height) of the read window. The<br>
> length of padfInArray == padfOutArray == nWindowXSize * nWindowYSize<br>
> *<br>
> * @param pData Process-specific data. This data is passed straight through<br>
> to the GDALRasterProcessFn and may contain e.g user defined parameters.<br>
> * The GDALRasterProcessFn definition would define the structure/type of<br>
> such data.<br>
> *<br>
> * @param pdfNoDataValue The no data value of the dataset<br>
> */<br>
><br>
> CPLErr GDALRunRasterProcess(GDALRasterProcessFn processFn, GDALDataset<br>
> *poSrcDataset,<br>
> GDALDataset *poDstDataset, void *pData, int<br>
> *pnWindowXSize,<br>
> int *pnWindowYSize, int *pnPixelBuffer);<br>
> /**<br>
> * \brief Apply a raster processing function to each sub-window of a raster.<br>
> *<br>
> * The input raster dataset is read in chunks of nWindowXSize * nWindowYSize<br>
> and each chunk is passed to the processing<br>
> * function. The output array from the function is written to the<br>
> destination dataset.<br>
> * An optional 'pixel buffer' can be specified to allow overlaps between<br>
> successive windows. This is useful for<br>
> * some algorithms, e.g. blob extraction, watershed/stream flow analysis,<br>
> convolution etc.<br>
> * Process specific data can be passed (e.g. configuration parameters). This<br>
> data is simply passed straight through to the processing<br>
> * function on each call.<br>
> *<br>
> * @param processFn A GDALRasterProcessFn to apply to each sub window of the<br>
> raster.<br>
> *<br>
> * @param poSrcDataset The source raster dataset from which pixel values are<br>
> read<br>
> *<br>
> * @param poDstDataset The destination raster dataset to which pixel values<br>
> are written. Must support RasterIO in write mode.<br>
> *<br>
> * @param pData Process-specific data. This is passed straight through to<br>
> the GDALRasterProcessFn on each call.<br>
> *<br>
> * @param pnWindowXSize The desired width of each read window. If NULL it<br>
> defaults to the 'natural' block size of the raster<br>
> *<br>
> * @param pnWindowYSize The desired height of each read window. If NULL it<br>
> defaults to the 'natural' block size.<br>
> *<br>
> * @param pnPixelBuffer A pixel buffer to apply to the read window. The read<br>
> window is expanded by pnPixelBuffer pixels in all directions such that<br>
> * each window overlaps by pnPixelBuffer pixels. Ignored if NULL or 0<br>
> *<br>
> * @return a CPLErr enum indicating whether the function succesfully ran.<br>
> */<br>
><br>
><br>
> CPLErr GDALRunRasterProcesses(GDALRasterProcessFn *paProcessFn, int<br>
> nProcesses, GDALDataset *poSrcDataset,<br>
> GDALDataset *poDstDataset, void **paData, int<br>
> *pnWindowXSize,<br>
> int *pnWindowYSize, int *pnPixelBuffer);<br>
><br>
> /**<br>
> * \brief Apply multiple raster processing functions to each sub-window of a<br>
> raster<br>
> *<br>
> * For each window, the functions defined by the paProcessFn array are<br>
> called in turn, with the array output of the previous function forming the<br>
> input * to the next function. This allows processing 'toolchains' to be<br>
> built without having to create intermediate datasets, which can be less<br>
> efficient in time and space.<br>
> *<br>
> *<br>
> * @param paProcessFn An array of GDALRasterProcessFn to apply to each sub<br>
> window of the raster<br>
> *<br>
> * @param nProcesses The size of paProcessFn<br>
> *<br>
> * @param poSrcDataset The source raster dataset from which pixel values are<br>
> read<br>
> *<br>
> * @param poDstDataset The destination raster dataset to which pixel values<br>
> are written. Must support RasterIO in write mode.<br>
> *<br>
> * @param paData an array of process-specific data objects of size<br>
> nProcesses. Each data object will be passed to the corresponding<br>
> GDALRasterProcessFn<br>
> *<br>
> * @param pnWindowXSize The desired width of each read window. If NULL it<br>
> defaults to the 'natural' block size of the raster<br>
> *<br>
> * @param pnWindowYSize The desired height of each read window. If NULL it<br>
> defaults to the 'natural' block size.<br>
> *<br>
> * @param pnPixelBuffer A pixel buffer to apply to the read window. The read<br>
> window is expanded by pnPixelBuffer pixels in all directions such that<br>
> * each window overlaps by pnPixelBuffer pixels. If NULL, it is ignored.<br>
> *<br>
> * @return a CPLErr enum indicating whether the function succesfully ran.<br>
> */<br>
> /*******************************************************************/<br>
><br>
> So GDALRunRasterProcess would take care of I/O and looping through the<br>
> dataset in chunks, which can be configure by the user. Each chunk is<br>
> submitted to the processing function, which does 'stuff' and populates the<br>
> output array. GDALRunRasterProcess would write it to the output dataset<br>
> object.<br>
><br>
> An example of a processing function could be something like this:<br>
><br>
> /*******************************************************************/<br>
> typedef struct ReclassArgs {<br>
> int nReclassArgs;<br>
> double *padfMinBound;<br>
> double *padfMaxBound;<br>
> double *padfBinValue;<br>
><br>
> } ReclassArgs;<br>
><br>
><br>
> CPLErr Reclass(double *padfInArray, double *padfOutArray,<br>
> int nWindowXSize, int nWindowYSize, void *pData,<br>
> double *pdfNoDataValue)<br>
> {<br>
> ReclassArgs *poRargs = static_cast<ReclassArgs *>(pData);<br>
> double pixel;<br>
> for (int i = 0; i < nWindowXSize * nWindowYSize - 1; ++i){<br>
> pixel = padfInArray[i];<br>
> for (int k = 0; k < poRargs->nReclassArgs; ++k){<br>
> if ((pixel > poRargs->padfMinBound[k]) && (pixel <=<br>
> poRargs->padfMaxBound[k])){<br>
> padfOutArray[i] = poRargs->padfBinValue[k];<br>
> break;<br>
> }<br>
> }<br>
> }<br>
> return CE_None;<br>
> }<br>
> /*******************************************************************/<br>
><br>
> To use it you would:<br>
><br>
> /*******************************************************************/<br>
> // Open a source dataset<br>
> GDALDataset *poSrc = GDALOpen(...)<br>
><br>
> // Create some output dataset with same dimensions etc..<br>
> GDALDataset *poDst = ....<br>
><br>
> // Setup the reclass args<br>
> ReclassArgs *poRargs = ...<br>
><br>
> // Call the process<br>
> CPLErr retval = GDALRunRasterProcess(Reclass, poSrc, poDst, (void<br>
> *)poRargs, NULL, NULL, NULL);<br>
> /*******************************************************************/<br>
><br>
> GDAL could feature its own library of common functions...reclass, truncate,<br>
> blob extraction etc etc and a user of GDAL could write their own and have<br>
> GDAL process it.<br>
><br>
> It would be useful to expose enough to python via SWIG so that a<br>
> GDALRasterProcessFn could be defined in python - this would then allow<br>
> rapid development in python, but push the expensive work of looping over<br>
> dataset chunks to C++.<br>
> I dont know enough about SWIG yet to see whether this would be feasible.<br>
><br>
> However, even in C++, it allows the user to only think about their<br>
> algorithm and leave the boilerplate to GDAL.<br>
><br>
> I imagine all this being implemented in a 'sub library', something like<br>
> gdpl.h (geospatial/gdal dataset processing library), but sitting in the<br>
> GDAL source tree (then I can happily expose all the useful GDAL goodness<br>
> without being worried about users having to include gdal headers etc, and<br>
> avoid any work in having to abstract it away)<br>
><br>
> This is a first stab at defining something like this and is probably not<br>
> yet comprehensive enough.<br>
><br>
> Any thoughts?<br>
><br>
><br>
><br>
><br>
> --<br>
> View this message in context:<br>
> <a href="http://osgeo-org.1560.x6.nabble.com/GDAL-raster-processing-library-tp52759" rel="noreferrer" target="_blank">http://osgeo-org.1560.x6.nabble.com/GDAL-raster-processing-library-tp52759</a><br>
</div></div>> 48.html Sent from the GDAL - Dev mailing list archive at Nabble.com.<br>
<div><div>> _______________________________________________<br>
> gdal-dev mailing list<br>
> <a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
> <a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" rel="noreferrer" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><span><font color="#888888"><br>
<br>
</font></span></div></div><span><font color="#888888"><span><font color="#888888">--<br>
Spatialys - Geospatial professional services<br>
<a href="http://www.spatialys.com" rel="noreferrer" target="_blank">http://www.spatialys.com</a><br>
</font></span></font></span></blockquote></div><br></div>
<br>_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" rel="noreferrer" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br></blockquote></div><br><br clear="all"><br>-- <br></div></div><div data-smartmail="gmail_signature"><div dir="ltr"><div>----------------------------------------------------------------------------------------------------------------<br>Peter J Halls, PhD Student, Post-war Reconstruction and Development Unit (PRDU),<br> University of York<br><br>Snail mail: c/o Research Centre for the Social Sciences (RCSS),<br> University of York, <br> Heslington, York YO10 5DD<br><br>This message has the status of a private and personal communication<br>----------------------------------------------------------------------------------------------------------------</div></div></div>
</div>
</blockquote></div><br></div>