[gdal-dev] Raster calculation optimalisation

[Paul Meems]

 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170613/a514624c/attachment.html>