[gdal-dev] Raster calculation optimalisation
Damian Dixon
damian.dixon at gmail.com
Tue Jun 13 06:36:11 PDT 2017
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 <bontepaarden at gmail.com> 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 <damian.dixon at gmail.com>:
>
>> 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 <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/5f7fb741/attachment-0001.html>
More information about the gdal-dev
mailing list