[gdal-dev] Raster calculation optimalisation

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

[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


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Raster calculation optimalisation

Damian Dixon
Hi Paul,

What is the time to run your code without applying the calculation to each pixel?

It is usually better to process the pixels in a block rather than walking across each row especially if you are processing a TIFF as these are usually stored as tiles (blocks).

The crossing from C# to C++ will also impact the performance so take a look at the time it takes to set the pixels. There may also be a method to set a block of pixels.

Sorry I can't help you more than some general advice and hints.

Regards
Damian



On 13 June 2017 at 12:20, Paul Meems <[hidden email]> 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
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Raster calculation optimalisation

Damian Dixon
​Hi Paul,

I can only do snippets in C++...

You can get the block size using:

  poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );

then do a read something like:

  for(int band = 0; band <noBands;band++)
  {
    poBand[0]->AdviseRead(startX+xoffset,startY+yoffset,blockWidth,blockHeight,blockWidth,blockHeight,GDT_Int32, NULL);
    CPLErr err = poBand[band]->RasterIO(GF_Read,startX+xoffset,startY+yoffset,blockWidth,blockHeight,buffer[band],blockWidth,blockHeight,GDT_Int32,0,0);
    if(err != CE_None)
    {
      report an error
    }
  }

I don't normally put raster back into GDAL we just use GDAL to read data and do all the processing ourselves.

I suspect all you have to do is change GF_Read to GF_Write and a few other tweaks :>>.

Regards
Damian
The following page may be of help: http://www.gdal.org/gdal_tutorial.html



On 13 June 2017 at 14:19, Paul Meems <[hidden email]> wrote:
Thanks Damian for your hint.
I'll try in a minute without the formula.
You mention process the pixels in a block. How would I do that? Do you have an example?

Thanks,

Paul

2017-06-13 15:05 GMT+02:00 Damian Dixon <[hidden email]>:
Hi Paul,

What is the time to run your code without applying the calculation to each pixel?

It is usually better to process the pixels in a block rather than walking across each row especially if you are processing a TIFF as these are usually stored as tiles (blocks).

The crossing from C# to C++ will also impact the performance so take a look at the time it takes to set the pixels. There may also be a method to set a block of pixels.

Sorry I can't help you more than some general advice and hints.

Regards
Damian



On 13 June 2017 at 12:20, Paul Meems <[hidden email]> 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
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev




_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Raster calculation optimalisation

janhec
In reply to this post by Paul Meems
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 <[hidden email]> 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
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Raster calculation optimalisation

Rutger
In reply to this post by Damian Dixon
Damian Dixon wrote
It is usually better to process the pixels in a block rather than walking
across each row especially if you are processing a TIFF as these are
usually stored as tiles (blocks).
Other layouts are common as well. For example, the Landsat TIFF's provided by the USGS have a row based-layout.  If you can choose it yourself, bocks are prefered in my opinion, since GDAL VRT's have a fixed blocksize of128*128. So when writing TIFF's, setting "TILED=YES" is a good default.

I think your spot on by mentioning the blocks. Don't assume the layout at all, look at the blocksize of the file and use it. If the blocks are relatively small (memoy-wise), using a multiple of the size can increase performance a bit more. So if its row-based and you have plenty of memory to spare, why not read blocks of (xsize, 128). Or if the blocksize is 128x128, use blocks of 256x256 etc.

If the volume of data is really large, increasing GDAL's block cache can be helpful. Although its best to avoid relying on the cache (if possible) by specifying an appropriate blocksize.

Here are a few mini-benchmarks:
http://nbviewer.jupyter.org/gist/RutgerK/4faa43873ee6389873dd4de85de7c450

https://stackoverflow.com/questions/41742162/gdal-readasarray-for-vrt-extremely-slow/41853759#41853759

Regards,
Rutger
Reply | Threaded
Open this post in threaded view
|

Re: Raster calculation optimalisation

Paul Meems
Thanks all for your suggestions.

@Rutger and @Damian:
Thanks for your suggestion about the blocks.
I had a look at Rutger's links,
I create the input file myself so I can add 'TILES=YES' but I'm not sure how to change my calculation code.
I see in this first link xbs, ybs = b1.GetBlockSize() But I don't see when to use the xbs or ybs variables.
I assume I need to change the reading of the data: band.ReadRaster(0, line, sizeX, 1, scanline, sizeX, 1, 0, 0); but I'm not sure how.
And how should I write back the block?

BTW. The main bottle-neck seems to be in the used formula parser.
With the formula it takes 68s, without the formule - just setting the pixelvalue takes 3.3s and with the formula written in code pixelValue = (float)Math.Exp(3.1 + 0.9 * f); it only takes 3.4s.
So not using mathparser has the highest benefits. I will do that first, but I also want to understand the block reading and writing.

@Jan, thanks for your link to the other parser. I had a quick look and it looks very promising. Sadly I couldn't get their C# example working.
I will look at as well.

Just a more general question.
Doesn't it makes sense if gdal would provide a gdal_calculate tool which also is librified like VectorTranslate?
It seems lots of people are implementing it on their own.

Thanks,

Paul


2017-06-14 9:29 GMT+02:00 Rutger <[hidden email]>:
Damian Dixon wrote
> It is usually better to process the pixels in a block rather than walking
> across each row especially if you are processing a TIFF as these are
> usually stored as tiles (blocks).

Other layouts are common as well. For example, the Landsat TIFF's provided
by the USGS have a row based-layout.  If you can choose it yourself, bocks
are prefered in my opinion, since GDAL VRT's have a fixed blocksize
of128*128. So when writing TIFF's, setting "TILED=YES" is a good default.

I think your spot on by mentioning the blocks. Don't assume the layout at
all, look at the blocksize of the file and use it. If the blocks are
relatively small (memoy-wise), using a multiple of the size can increase
performance a bit more. So if its row-based and you have plenty of memory to
spare, why not read blocks of (xsize, 128). Or if the blocksize is 128x128,
use blocks of 256x256 etc.

If the volume of data is really large, increasing GDAL's block cache can be
helpful. Although its best to avoid relying on the cache (if possible) by
specifying an appropriate blocksize.

Here are a few mini-benchmarks:
http://nbviewer.jupyter.org/gist/RutgerK/4faa43873ee6389873dd4de85de7c450

https://stackoverflow.com/questions/41742162/gdal-readasarray-for-vrt-extremely-slow/41853759#41853759

Regards,
Rutger




--
View this message in context: http://osgeo-org.1560.x6.nabble.com/gdal-dev-Raster-calculation-optimalisation-tp5324014p5324120.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Raster calculation optimalisation

Damian Dixon
Personally I would agree with your approach at looking at the cost of running the
formula first.


However if you do change your code to read/write blocks than there are additional
opportunities for speeding the processing up.

Once you have the block size reading working consider using multiple threads.
You will need to use a GDALDataset for each thread.

You will however need to consider how you write the results out again
as you will not be able to write to a single output. You will have to write to
a separate raster for each thread.

This is the approach I have taken. 

I can't help you with writing the block back as once I've got the raster data I
processes it into a format that allows me to load and display the data a lot
quicker than going back to GDAL for the data.

You need to find the equivalent to:
Also these may help...



On 15 June 2017 at 14:04, Paul Meems <[hidden email]> wrote:
Thanks all for your suggestions.

@Rutger and @Damian:
Thanks for your suggestion about the blocks.
I had a look at Rutger's links,
I create the input file myself so I can add 'TILES=YES' but I'm not sure how to change my calculation code.
I see in this first link xbs, ybs = b1.GetBlockSize() But I don't see when to use the xbs or ybs variables.
I assume I need to change the reading of the data: band.ReadRaster(0, line, sizeX, 1, scanline, sizeX, 1, 0, 0); but I'm not sure how.
And how should I write back the block?

BTW. The main bottle-neck seems to be in the used formula parser.
With the formula it takes 68s, without the formule - just setting the pixelvalue takes 3.3s and with the formula written in code pixelValue = (float)Math.Exp(3.1 + 0.9 * f); it only takes 3.4s.
So not using mathparser has the highest benefits. I will do that first, but I also want to understand the block reading and writing.

@Jan, thanks for your link to the other parser. I had a quick look and it looks very promising. Sadly I couldn't get their C# example working.
I will look at as well.

Just a more general question.
Doesn't it makes sense if gdal would provide a gdal_calculate tool which also is librified like VectorTranslate?
It seems lots of people are implementing it on their own.

Thanks,

Paul


2017-06-14 9:29 GMT+02:00 Rutger <[hidden email]>:
Damian Dixon wrote
> It is usually better to process the pixels in a block rather than walking
> across each row especially if you are processing a TIFF as these are
> usually stored as tiles (blocks).

Other layouts are common as well. For example, the Landsat TIFF's provided
by the USGS have a row based-layout.  If you can choose it yourself, bocks
are prefered in my opinion, since GDAL VRT's have a fixed blocksize
of128*128. So when writing TIFF's, setting "TILED=YES" is a good default.

I think your spot on by mentioning the blocks. Don't assume the layout at
all, look at the blocksize of the file and use it. If the blocks are
relatively small (memoy-wise), using a multiple of the size can increase
performance a bit more. So if its row-based and you have plenty of memory to
spare, why not read blocks of (xsize, 128). Or if the blocksize is 128x128,
use blocks of 256x256 etc.

If the volume of data is really large, increasing GDAL's block cache can be
helpful. Although its best to avoid relying on the cache (if possible) by
specifying an appropriate blocksize.

Here are a few mini-benchmarks:
http://nbviewer.jupyter.org/gist/RutgerK/4faa43873ee6389873dd4de85de7c450

https://stackoverflow.com/questions/41742162/gdal-readasarray-for-vrt-extremely-slow/41853759#41853759

Regards,
Rutger




--
View this message in context: http://osgeo-org.1560.x6.nabble.com/gdal-dev-Raster-calculation-optimalisation-tp5324014p5324120.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev