search

qgis / QGIS

QGIS is a free, open source, cross platform (lin/win/mac) geographical information system (GIS)
https://qgis.org
GNU General Public License v2.0
10.54k stars 2.99k forks source link

GDAL Raster Calculator in the Processing Toolbox producing some incorrect results #42724

Open tombrennan06 opened 3 years ago

tombrennan06 commented 3 years ago

Describe the bug The GDAL Raster Calculator in the Processing Toolbox is producing incorrect results for certain calculations.

For example, when a raster with data from 1 to 255 has the formula A>128 applied via the GDAL Raster Calculator, you would expect data points greater than 128 to get the value 1 and those less than or equal to 128 to be 0. Instead, the output is much the same as the input, but with values from 1 to 254.

How to Reproduce

The input dataset is a simple 100 x 100 gradient raster with values from 255 (left) to 1 (right) raster_calc_issue.zip

  1. Open the attached project file
  2. In the Processing Toolbox, open the GDAL Raster Calculator tool
  3. Select Input Layer A: raster_calc_sample, Calculation: A>128, Output Raster Type: Byte, and then Run 2021_04_08_22_43_41_Raster_Calculator

Expected output: a 0/1 raster with the right half 0 and the left half 1 Actual output: a 0-254 raster

QGIS and OS versions

QGIS version 3.18.1-Zürich QGIS code revision 202f1bf7e5
Compiled against Qt 5.11.2 Running against Qt 5.11.2
Compiled against GDAL/OGR 3.1.4 Running against GDAL/OGR 3.1.4
Compiled against GEOS 3.8.1-CAPI-1.13.3 Running against GEOS 3.8.1-CAPI-1.13.3
Compiled against SQLite 3.29.0 Running against SQLite 3.29.0
PostgreSQL Client Version 11.5 SpatiaLite Version 4.3.0
QWT Version 6.1.3 QScintilla2 Version 2.10.8
Compiled against PROJ 6.3.2 Running against PROJ Rel. 6.3.2, May 1st, 2020
OS Version Windows 10 (10.0)
Active python plugins processing_wbt; qgis_resource_sharing; QuickOSM; db_manager; MetaSearch; processing

Additional context

gioman commented 3 years ago

Expected output: a 0/1 raster with the right half 0 and the left half 1 Actual output: a 0-254 raster

@tombrennan06 works fine here on QGIS master/linux

Screenshot_20210408_160908

nicogodet commented 3 years ago

I can confirm the bug on 3.18.1 with "A>128", output format does not matter. Works fine with greater(A, 128).

Works fine on 3.16.5

Gdal 3.2.2 on OSGeo4W Testing Win10

gioman commented 3 years ago

I can confirm the bug on 3.18.1

@nicogodet can you test master?

roya0045 commented 3 years ago

@gioman I tried with the mingw artefact, can't get the gdal file to perform, but the native calculator works fine. Were your tests on the gdal variant?

gioman commented 3 years ago

@roya0045 yes sure, I used the GDAL raster calc.

nicogodet commented 3 years ago

@gioman bug on master too.

gioman commented 3 years ago

@gioman bug on master too.

@nicogodet same correct result for me on 3.18.1 on Windows 10

Screenshot_20210408_202048

nicogodet commented 3 years ago

I installed 3.18.1 with the new standalone installer (.msi) to compare.

On OSGeo4W : 

QGIS version: 3.19.0-Master
QGIS code revision: ca7c033a4
Qt version: 5.15.2
GDAL version: 3.2.2
GEOS version: 3.9.1-CAPI-1.14.2
PROJ version: Rel. 8.0.0, March 1st, 2021
PDAL version: 2.2.0 (git-version: 24dd45)
Processing algorithm…
Algorithm 'Raster calculator' starting…
Input parameters:
{ 'BAND_A' : 1, 'BAND_B' : None, 'BAND_C' : None, 'BAND_D' : None, 'BAND_E' : None, 'BAND_F' : None, 'EXTRA' : '', 'FORMULA' : 'A>128', 'INPUT_A' : 'C:/Users/GODET/AppData/Local/Temp/processing_bbedrI/372c8ab275394b2eae455f6bf51e85a2/OUTPUT.tif', 'INPUT_B' : None, 'INPUT_C' : None, 'INPUT_D' : None, 'INPUT_E' : None, 'INPUT_F' : None, 'NO_DATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'RTYPE' : 5 }

GDAL command:
gdal_calc.bat --overwrite --calc "A>128" --format GTiff --type Float32 -A C:/Users/GODET/AppData/Local/Temp/processing_bbedrI/372c8ab275394b2eae455f6bf51e85a2/OUTPUT.tif --A_band 1 --outfile C:/Users/GODET/AppData/Local/Temp/processing_bbedrI/0d0f61c73b8e43d8abc3e5c5c4872d7d/OUTPUT.tif
GDAL command output:
Process completed successfully
Execution completed in 0.68 seconds
Results:
{'OUTPUT': 'C:/Users/GODET/AppData/Local/Temp/processing_bbedrI/0d0f61c73b8e43d8abc3e5c5c4872d7d/OUTPUT.tif'}

Loading resulting layers
Algorithm 'Raster calculator' finished

See GDAL command output is empty !

On standalone install:

QGIS version: 3.18.1-Zürich
QGIS code revision: 202f1bf7e5
Qt version: 5.15.2
GDAL version: 3.2.2
GEOS version: 3.9.1-CAPI-1.14.2
PROJ version: Rel. 8.0.0, March 1st, 2021
PDAL version: 2.2.0 (git-version: 24dd45)
Processing algorithm…
Algorithm 'Raster calculator' starting…
Input parameters:
{ 'BAND_A' : 1, 'BAND_B' : None, 'BAND_C' : None, 'BAND_D' : None, 'BAND_E' : None, 'BAND_F' : None, 'EXTRA' : '', 'FORMULA' : 'A>128', 'INPUT_A' : 'D:/raster_calc_sample.tif', 'INPUT_B' : None, 'INPUT_C' : None, 'INPUT_D' : None, 'INPUT_E' : None, 'INPUT_F' : None, 'NO_DATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'RTYPE' : 5 }

GDAL command:
gdal_calc.bat --calc "A>128" --format GTiff --type Float32 -A D:/raster_calc_sample.tif --A_band 1 --outfile C:/Users/GODET/AppData/Local/Temp/processing_IOytHx/2d2464a2133c4634b852133238777765/OUTPUT.tif
GDAL command output:
Traceback (most recent call last):
File "C:/PROGRA~1/QGIS31~1.1/apps/qgis/./python/plugins\processing\algs\gdal\GdalUtils.py", line 130, in on_stderr
val = ba.data().decode('UTF-8')
UnicodeDecodeError: 'utf-8' codec can't decode byte 0x8a in position 3: invalid start byte

The above exception was the direct cause of the following exception:

SystemError: <class 'PyQt5.QtCore.QByteArray'> returned a result with an error set

The above exception was the direct cause of the following exception:

SystemError: <class 'PyQt5.QtCore.QByteArray'> returned a result with an error set

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
File "C:/PROGRA~1/QGIS31~1.1/apps/qgis/./python/plugins\processing\algs\gdal\GdalAlgorithm.py", line 135, in processAlgorithm
GdalUtils.runGdal(commands, feedback)
File "C:/PROGRA~1/QGIS31~1.1/apps/qgis/./python/plugins\processing\algs\gdal\GdalUtils.py", line 146, in runGdal
res = proc.run(feedback)
SystemError: <built-in method run of QgsBlockingProcess object at 0x0000022FC1E7E160> returned a result with an error set

Execution failed after 0.08 seconds

Loading resulting layers
Algorithm 'Raster calculator' finished

With greater(A, 128) on both install:

QGIS version: 3.19.0-Master
QGIS code revision: ca7c033a4
Qt version: 5.15.2
GDAL version: 3.2.2
GEOS version: 3.9.1-CAPI-1.14.2
PROJ version: Rel. 8.0.0, March 1st, 2021
PDAL version: 2.2.0 (git-version: 24dd45)
Processing algorithm…
Algorithm 'Raster calculator' starting…
Input parameters:
{ 'BAND_A' : 1, 'BAND_B' : None, 'BAND_C' : None, 'BAND_D' : None, 'BAND_E' : None, 'BAND_F' : None, 'EXTRA' : '', 'FORMULA' : 'greater(A,128)', 'INPUT_A' : 'C:/Users/GODET/AppData/Local/Temp/processing_bbedrI/372c8ab275394b2eae455f6bf51e85a2/OUTPUT.tif', 'INPUT_B' : None, 'INPUT_C' : None, 'INPUT_D' : None, 'INPUT_E' : None, 'INPUT_F' : None, 'NO_DATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'RTYPE' : 5 }

GDAL command:
gdal_calc.bat --overwrite --calc "greater(A,128)" --format GTiff --type Float32 -A C:/Users/GODET/AppData/Local/Temp/processing_bbedrI/372c8ab275394b2eae455f6bf51e85a2/OUTPUT.tif --A_band 1 --outfile C:/Users/GODET/AppData/Local/Temp/processing_bbedrI/504796b75b60402cbba14d68b3352bb7/OUTPUT.tif
GDAL command output:
0.. 20.. 40.. 60.. 80.. 100 - Done
Process completed successfully
Execution completed in 0.58 seconds
Results:
{'OUTPUT': 'C:/Users/GODET/AppData/Local/Temp/processing_bbedrI/504796b75b60402cbba14d68b3352bb7/OUTPUT.tif'}

Loading resulting layers
Algorithm 'Raster calculator' finished

If I print the value of ba and ba.decode() from standalone install:

b'Acc\x8as refus\x82.\r\n'
b'Acc\x8as refus\x82.\r\n'

If I print the value of ba and ba.decode() from OSgeo4W install:

b''
b''

Locale issue

ping @nyalldawson as you implemented this in https://github.com/qgis/QGIS/pull/40897

@gioman could you set you win VM in French or any accented language to confirm ?

nicogodet commented 3 years ago

https://github.com/qgis/QGIS/blob/3ee49a60a382b7fc8e97e5cf7796620efb950ef6/python/plugins/processing/algs/gdal/GdalUtils.py#L141 It seems to break the command. Output is: ['gdal_calc.bat', '--overwrite', '--calc', 'A>128', '--format', 'GTiff', '--type', 'Float32', '-A', 'D:/raster_calc_sample.tif', '--A_band', '1', '--outfile', 'C:/Users/GODET/AppData/Local/Temp/processing_hZjELV/1c1ad78c8ada42848f48772fbb57310f/OUTPUT.tif']

If I run gdal_calc.bat --overwrite --calc A>128 --format GTiff --type Float32 -A D:/raster_calc_sample.tif --A_band 1 --outfile D:/raster_calc_sample_truc.tif in a shell, It fails and it results in a file named 128

C:\OSGeo4W64_new>.\apps\Python39\Scripts\gdal_calc.bat --overwrite --calc A>128 --format GTiff --type Float32 -A D:/raster_calc_sample.tif --A_band 1 --outfile D:/raster_calc_sample_truc.tif

C:\OSGeo4W64_new>dir
 Le volume dans le lecteur C s’appelle OS
 Le numéro de série du volume est CABC-B79C

 Répertoire de C:\OSGeo4W64_new

04/09/2021  09:03 AM    <DIR>          .
04/09/2021  09:03 AM    <DIR>          ..
04/09/2021  09:03 AM                36 128
04/08/2021  05:35 PM    <DIR>          apps
04/08/2021  05:43 PM    <DIR>          bin
03/12/2021  05:32 PM    <DIR>          doc
04/08/2021  05:33 PM    <DIR>          etc
03/12/2021  05:32 PM    <DIR>          lib
02/18/2021  02:33 PM               321 OSGeo4W.bat
10/09/2020  11:59 PM             5,430 OSGeo4W.ico
04/08/2021  05:33 PM    <DIR>          share
03/12/2021  05:37 PM    <DIR>          var
               3 fichier(s)            5,787 octets
               9 Rép(s)  80,020,582,400 octets libres

C:\OSGeo4W64_new>
gioman commented 3 years ago

I installed 3.18.1 with the new standalone installer (.msi) to compare.

@nicogodet so if you use the official Windows installer (the ones linked on QGIS's we site) no issue?

jef-n commented 3 years ago

If I run gdal_calc.bat --overwrite --calc A>128 --format GTiff --type Float32 -A D:/raster_calc_sample.tif --A_band 1 --outfile D:/raster_calc_sample_truc.tif in a shell, It fails and it results in a file named 128

That's normal . You have to quote A>128, otherwise cmd will redirect stdout to 128 and gdal_calc.bat will only see A:

gdal_calc.bat --overwrite --calc "A>128" --format GTiff --type Float32 -A D:/raster_calc_sample.tif --A_band 1 --outfile D:/raster_calc_sample_truc.tif

But what that outputs probably has some content that is better handled by python 3.9 than by 3.7 (or PyQt or Qt)

nicogodet commented 3 years ago

@jef-n The unquoted A>128 is the output of https://github.com/qgis/QGIS/blob/3ee49a60a382b7fc8e97e5cf7796620efb950ef6/python/plugins/processing/algs/gdal/GdalUtils.py#L141

        print("fused_command")
        print(fused_command)
        print("QgsRunProcess.splitCommand(fused_command)")
        print(QgsRunProcess.splitCommand(fused_command))
fused_command
gdal_calc.bat --calc "A>128" --format GTiff --type Float32 -A D:\raster_calc_sample.tif --A_band 1 --outfile C:/Users/GODET/AppData/Local/Temp/processing_VZTszq/016ad59a3ab2455da63429d6ef3f23ac/OUTPUT.tif
QgsRunProcess.splitCommand(fused_command)
['gdal_calc.bat', '--calc', 'A>128', '--format', 'GTiff', '--type', 'Float32', '-A', 'D:\\raster_calc_sample.tif', '--A_band', '1', '--outfile', 'C:/Users/GODET/AppData/Local/Temp/processing_VZTszq/016ad59a3ab2455da63429d6ef3f23ac/OUTPUT.tif']

@gioman Same issue on Official Windows installer But the refused access is caused by the console trying to create 128 in a folder that require an admin right.

gioman commented 3 years ago

@gioman Same issue on Official Windows installer

@nicogodet so it not clear to me why I can't replicate on my machine.

jef-n commented 3 years ago

gdal_calc.bat --calc "A>128" --format GTiff --type Float32 -A D:\raster_calc_sample.tif --A_band 1 --outfile C:/Users/GODET/AppData/Local/Temp/processing_VZTszq/016ad59a3ab2455da63429d6ef3f23ac/OUTPUT.tif

That output has the quotes. But I was referring to what you entered into the shell and that was missing the quotes. Running that with quotes might show what in the output might causes the hickup.

nicogodet commented 3 years ago

@jef-n Look at the output of QgsRunProcess.splitCommand(fused_command) and you will see that quotes are missing in the list of arguments. This list of arguments is then throw in https://github.com/qgis/QGIS/blob/3ee49a60a382b7fc8e97e5cf7796620efb950ef6/python/plugins/processing/algs/gdal/GdalUtils.py#L142

jef-n commented 3 years ago

@jef-n Look at the output of QgsRunProcess.splitCommand(fused_command) and you will see that quotes are missing in the list of arguments.

I'm not sure what your point is. I commented on what you ran manually. I'm interested in it's output, but with the quotes. AFAICT the command is run fine, but the "old" standalone can't parse it's output, while the "new" one can.

nicogodet commented 3 years ago

The output with quotes is not the actual command executed.

Of course if I execute manually the command with quotes everything runs correctly.

But. The command with quotes is sent to QgsRunProcess.splitCommand() then QgsBlockingProcess(). The problem is that splitCommand() removes quotes. That's the purpose of my previous comments.

jef-n commented 3 years ago

Of course if I execute manually the command with quotes everything runs correctly.

And it outputs?

nicogodet commented 3 years ago

Outputs are corrects.

nicogodet commented 3 years ago

If I run avec A > 128 output is correct

QGIS version: 3.18.1-Zürich
QGIS code revision: 202f1bf7e5
Qt version: 5.15.2
GDAL version: 3.2.2
GEOS version: 3.9.1-CAPI-1.14.2
PROJ version: Rel. 8.0.0, March 1st, 2021
PDAL version: 2.2.0 (git-version: 24dd45)
Processing algorithm…
Algorithm 'Raster calculator' starting…
Input parameters:
{ 'BAND_A' : 1, 'BAND_B' : None, 'BAND_C' : None, 'BAND_D' : None, 'BAND_E' : None, 'BAND_F' : None, 'EXTRA' : '', 'FORMULA' : 'A > 128', 'INPUT_A' : 'D:/raster_calc_sample.tif', 'INPUT_B' : None, 'INPUT_C' : None, 'INPUT_D' : None, 'INPUT_E' : None, 'INPUT_F' : None, 'NO_DATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'RTYPE' : 5 }

GDAL command:
gdal_calc.bat --calc "A > 128" --format GTiff --type Float32 -A D:\raster_calc_sample.tif --A_band 1 --outfile C:/Users/GODET/AppData/Local/Temp/processing_EuzHpE/3f8ccc7a6c2d41708226f60791fd86fc/OUTPUT.tif
GDAL command output:
0.. 20.. 40.. 60.. 80.. 100 - Done
Process completed successfully
Execution completed in 1.21 seconds
Results:
{'OUTPUT': 'C:/Users/GODET/AppData/Local/Temp/processing_EuzHpE/3f8ccc7a6c2d41708226f60791fd86fc/OUTPUT.tif'}

Loading resulting layers
Algorithm 'Raster calculator' finished

image

A>128 incorrect

nicogodet commented 3 years ago

@tombrennan06 could you test with A > 128 (space around >) instead of A>128 ?

gioman commented 3 years ago

@tombrennan06 could you test with A > 128 (space around >) instead of A>128 ?

@nicogodet for what is worth... correct results on 3.18.1 on both Ubuntu and Windows (using the official installer).

nicogodet commented 3 years ago

@gioman Even without spaces around > ?

gioman commented 3 years ago

@gioman Even without spaces around > ?

@nicogodet yes.

gioman commented 3 years ago

A>128 incorrect

@nicogodet why you say that the image is incorrect: if you change the project background to something different from white you'll se that the left part of the output is white (1) and the right part is black (0).

nicogodet commented 3 years ago

@gioman No. With A>128 without spaces the command executed is gdal_calc.bat --overwrite --calc A>128 --format GTiff --type Float32 -A D:/raster_calc_sample.tif --A_band 1 --outfile D:/raster_calc_sample_truc.tif so in a shell, gdal_calc.bat --overwrite --calc A redirected to 128.

image Same with spaces image

I can not give more info, it's out of my competences... But I think all my previous messages show where is the error (https://github.com/qgis/QGIS/issues/42724#issuecomment-816524879) but I can not explain why the output is correct when you run A>128...

tombrennan06 commented 3 years ago

@nicogodet - same result as you.

With A>128 (no spaces), I get the same error as I originally described.

With A > 128 (spaces), I get the expected (correct) output (half white/half black).

FWIW, I am using the official installer for 3.18.1

github-actions[bot] commented 3 years ago

The QGIS project highly values your report and would love to see it addressed. However, this issue has been left in feedback mode for the last 14 days and is being automatically marked as "stale". If you would like to continue with this issue, please provide any missing information or answer any open questions. If you could resolve the issue yourself meanwhile, please leave a note for future readers with the same problem and close the issue. In case you should have any uncertainty, please leave a comment and we will be happy to help you proceed with this issue. If there is no further activity on this issue, it will be closed in a week.

nicogodet commented 3 years ago

@gioman I think feedback label could be removed.

gioman commented 3 years ago

@gioman I think feedback label could be removed.

@nicogodet done. It does not affect QGIS on Linux.

nicogodet commented 3 years ago

@gioman Indeed, it could be a Windows issue linked to QProcess and its Windows specificities regarding command arguments.

Windows: The arguments are quoted and joined into a command line that is compatible with the CommandLineToArgvW() Windows function. For programs that have different command line quoting requirements, you need to use setNativeArguments(). One notable program that does not follow the CommandLineToArgvW() rules is cmd.exe and, by consequence, all batch scripts.

https://doc.qt.io/qt-5/qprocess.html#start