[gdal-dev] Raster calculation optimalisation

[Jan Heckman]

Maybe it does not compile the formula but runs parser etc. as a kind of
interpreter at each function call.
Googled and saw http://beltoforion.de/article.php?a=muparsersse. Did not
try, so no recommendation, let alone guarantee ;)
You could check by defining a test function explicitly, if it is much
faster when compiled directly, that should be it, I guess.
Good luck,
Jan

On Tue, Jun 13, 2017 at 1:20 PM, Paul Meems <bontepaarden at gmail.com> wrote:

> I'm using GDAL v2.1.3 with the SWIG bindings in my custom C# application.
>
> I've created a method to do some calculation of my raster file:
> public bool GdalCalculate(string input, string output, string formula,
> double? minValue = null)
> {
>     if (!File.Exists(input))
>         throw new FileNotFoundException("Can't find the input file", new
> Exception("Working with " + input));
>
>     // 1. First copy the input file to create the output file
>     try
>     {
>         /* --------------------------------------------------------------------
> */
>         /*      Get driver
>      */
>         /* --------------------------------------------------------------------
> */
>         using (var drv = Gdal.GetDriverByName("GTiff"))
>         {
>             if (drv == null)
>             {
>                 throw new Exception("Can't get GTiff driver");
>             }
>
>             /* ------------------------------
> -------------------------------------- */
>             /*      Open dataset.
>           */
>             /* ------------------------------
> -------------------------------------- */
>             using (var ds = Gdal.Open(input, Access.GA_ReadOnly))
>             {
>                 if (ds == null)
>                 {
>                     throw new Exception("Can't open GDAL dataset: " +
> input);
>                 }
>
>                 var options = new[] { "" };
>                 using (var newDataset = drv.CreateCopy(output, ds, 0,
> options, null, "Sample Data"))
>                 {
>                     if (newDataset == null)
>                     {
>                         throw new Exception("Can't create destination
> dataset: " + output);
>                     }
>
>                     // Close input dataset:
>                     ds.Dispose();
>
>                     // 2. Loop through all pixels and perform formula on
> it:
>                     using (var band = newDataset.GetRasterBand(1))
>                     {
>                         double noData;
>                         int hasValue;
>                         band.GetNoDataValue(out noData, out hasValue);
>                         var sizeX = band.XSize;
>                         var numLines = band.YSize;
>                         var func = new Function("f(A) = " + formula);
>
>                         // Loop through each line in turn.
>                         for (var line = 0; line < numLines; line++)
>                         {
>                             var scanline = new float[sizeX];
>                             // Read in data for the current line
>                             var cplReturn = band.ReadRaster(0, line,
> sizeX, 1, scanline, sizeX, 1, 0, 0);
>                             if (cplReturn != CPLErr.CE_None)
>                             {
>                                 throw new Exception("band.ReadRaster
> failed: " + Gdal.GetLastErrorMsg());
>                             }
>
>                             var outputLine = new List<float>();
>                             foreach (var f in scanline)
>                             {
>                                 var pixelValue = f;
>                                 if ((float)f != (float)noData)
>                                 {
>                                     // Calculate
>                                     pixelValue = (float)func.calculate(f);
>                                     if (minValue.HasValue)
>                                         pixelValue =
> (float)Math.Max(pixelValue, minValue.GetValueOrDefault());
>                                 }
>                                 outputLine.Add(pixelValue);
>                             }
>
>                             // Rewrite line:
>                             cplReturn = band.WriteRaster(0, line, sizeX,
> 1, outputLine.ToArray(), sizeX, 1, 0, 0);
>                             if (cplReturn != CPLErr.CE_None)
>                             {
>                                 throw new Exception("band.WriteRaster
> failed: " + Gdal.GetLastErrorMsg());
>                             }
>                         }
>
>                         // 3. Save changes to file:
>                         band.FlushCache();
>                         newDataset.FlushCache();
>                     }
>                 }
>             }
>         }
>     }
>     catch (Exception e)
>     {
>         throw new Exception("Can't open GDAL dataset: " + input, e);
>     }
>
>     return true;
> }
>
> This method is working fine, but is very slow.
> I'm using a 3621x4466 tiff file, which has a size of 67MB.
> It is taking around 10 minutes to process.
> (float)func.calculate(f); is using http://mathparser.org/
> An example of a function I'm using: 9.7E-05-e^(3.1 + 0.9 * A) where A is
> the pixel value.
>
> How can I optimize my function?
>
> Thanks,
>
> Paul
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170613/d0949938/attachment.html>