GDAL raster calculator RuntimeWarnings

[ar-siddiqui]

The runoff calculations here cause runtime errors.

https://github.com/Dewberry/QNSPECT/blob/fbff1a01a10e5e40d58d71b3c8c8302b032b6cbe/algorithms/run_analysis/Runoff_Volume.py#L138

<string>:1: RuntimeWarning: overflow encountered in square

This might be because of #29

• Checked that it is not caused by misaligned cells

[ar-siddiqui]

After checking and fixing #29, this is not because of it.

It is because of raster not able to store the resultant value which is large enough.

Possible solutions:

1. Change dtype to float64 (but there maybe cases when data is bigger than float 64, this will also mean that it would be unnecessarily float64 in some cases when it does not need to be)
2. Make dtype a user input (advanced parameter)
3. Change units of calculations and use bigger units.

More details:
https://stackoverflow.com/a/7559950/11428410

[ar-siddiqui]

@DaveEslinger @NSPECT do you remember what we settled for as a solution for this one?

[DaveEslinger]

I believe we were going to test doing a conversion before accumulation on these variables as well. I think it will be similar to what you did to correct the units on erosion, but we need to go from mg to kg, to 10^6 instead of 10^3.

Also, if there is already a conversion at the end, we'll need to drop that.

Note: My comments are about the pollutants. I suspect your previous unit correction will have fixed the erosion issue.

[ar-siddiqui]

After investigating, it turns out that this is not an overflow issue but GDAL/NUMPY calculation issue. #61 is also caused by GDAL/NUMPY calculation.

The raster values are nowhere close to overflow limit for float32. The error can be replicated outside of QGIS as well.

A temporary fix has been put in place in 97ae3f1 that may or may not work for a different test area.

A ticket has been filed with GDAL OSGeo/gdal#5609

Alternative tools were considered, like the QGIS processing raster calculator, but the idea was dropped because several bugs exist in the QGIS repo for the QGIS processing raster calculator.

It is important to note that these messages are warnings, and doesn't seem to impact results

[ar-siddiqui]

Tested results of GDAL calculator with and without warnings and they are exactly the same. Thus changing the label from bug to warnings.

[ar-siddiqui]

@DaveEslinger when you get a chance can you please confirm if you are consistently getting this warning message at the same place in erosion modules for different test datasets.

[DaveEslinger]

With a small 99x100 cell (cell size 30 m) artificial data set, I do NOT get the error at all. I do get an error about no NODATA values to fill. Note that I think I'm not really using good R-Factor values here since I'm getting 9,745 kg/year of soil from one 30x30m cell. That seems a lot.

With real data for Hawaii and Puerto Rico I DO get the error, consistently at the same spot. I'm getting new data to test later.

With HI data, get error here:
Generating local sediments raster ...
GDAL command:
gdal_calc.bat --overwrite --calc "A * B * 907.18474" --format GTiff --type Float32 --NoDataValue -999999.0 -A C:/Users/davee/AppData/Local/Temp/processing_MbaJAK/35d078bbe5f1406f9317da5d9429f351/OUTPUT.tif --A_band 1 -B C:/Users/davee/AppData/Local/Temp/processing_MbaJAK/bfc9ac4f95fa4b55a07db626327684af/OUTPUT.tif --B_band 1 --overwrite --outfile "E:\GIS\QNSPECT\Testing\base05142022\Sediment Local.tif"
GDAL command output:
:1: RuntimeWarning: overflow encountered in multiply
C:\OSGeo4W\apps\Python39\lib\site-packages\GDAL-3.4.2-py3.9-win-amd64.egg\osgeo_utils\gdal_calc.py:522: RuntimeWarning: invalid value encountered in multiply
myResult = ((1 * (myNDVs == 0)) * myResult) + (myOutNDV * myNDVs)

FOR PUERTO RICO DATA error occurs same spot:
Generating local sediments raster ...
GDAL command:
gdal_calc.bat --overwrite --calc "A * B * 907.18474" --format GTiff --type Float32 --NoDataValue -999999.0 -A C:/Users/davee/AppData/Local/Temp/processing_MbaJAK/0e9fa30bbaa44f1eb2869fa2ca003009/OUTPUT.tif --A_band 1 -B C:/Users/davee/AppData/Local/Temp/processing_MbaJAK/ae84be7b31cb43ebb684934f488db742/OUTPUT.tif --B_band 1 --overwrite --outfile "E:\GIS\QNSPECT\Testing\PR_clipped_base\Sediment Local.tif"
GDAL command output:
:1: RuntimeWarning: overflow encountered in multiply
C:\OSGeo4W\apps\Python39\lib\site-packages\GDAL-3.4.2-py3.9-win-amd64.egg\osgeo_utils\gdal_calc.py:522: RuntimeWarning: invalid value encountered in multiply
myResult = ((1 * (myNDVs == 0)) * myResult) + (myOutNDV * myNDVs)

WITH ARTIFICAL DATA, NO ERROR:
Generating local sediments raster ...
GDAL command:
gdal_calc.bat --overwrite --calc "A * B * 907.18474" --format GTiff --type Float32 --NoDataValue -999999.0 -A C:/Users/davee/AppData/Local/Temp/processing_MbaJAK/d0aba6b299b34237b1d81a7376282039/OUTPUT.tif --A_band 1 -B C:/Users/davee/AppData/Local/Temp/processing_MbaJAK/e161a91c2dc8432cb16aac10c0b29a76/OUTPUT.tif --B_band 1 --overwrite --outfile "E:\GIS\QNSPECT\Testing\100x99x30_base\Sediment Local.tif"
GDAL command output:
0.. 20.. 40.. 60.. 80.. 100 - Done
Process completed successfully

[ar-siddiqui]

GDAL