<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 class="h5"><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 class="HOEnZb"><div class="h5">> _______________________________________________<br>
> gdal-dev mailing list<br>
> <a href="mailto:gdal-dev@lists.osgeo.org">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>
<br>
</div></div><span class="HOEnZb"><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></blockquote></div><br></div>