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.52k stars 2.99k forks source link

Issue on GDAL/GRASS raster calculator (3.22) #47656

Closed HelderSamuelBR closed 2 years ago

HelderSamuelBR commented 2 years ago

What is the bug or the crash?

I installed QGIS this week for 3.22, but I was using 3.16. I tried using the raster calculator GDAL and r.mapcalc.simple, but I had problems with both. The SAGA and native QGIS raster calculators had no problems. In this case, I'm calculating the albedo, which varies from 0 to 1. Could the cause be some missing library? I did the advanced installation of QGIS.

UPDATE: By some reason, GDAL/r.mapcalc.simple are multiplying the real results by 10.000.

erro

Steps to reproduce the issue

  1. Get Landsat 8: Bands 2, 3, 4, 5, 6 and 7.
  2. Plugin SCP: LS8 TOAR
  3. Cut a region of interest.
  4. Open GDAL Raster Calculator: A > LS8_B2 [...] F > LS8_B7 (TOAR) Equation: (A 0.300)+(B 0.277 )+(C 0.233)+(D 0.143)+(E 0.036)+(F 0.012)
  5. Open GRASS r.mapcalc.simple A > LS8_B2 [...] F > LS8_B7 Equation: (A 0.300)+(B 0.277)+(C 0.233)+(D 0.143)+(E 0.036)+(F 0.012)

The same result appears in both.

Versions

Versão do QGIS | 3.22.4-Białowieża | Código da versão do QGIS | ce8e65e9 -- | -- | -- | -- Qt version | 5.15.2 Versão do Python | 3.9.5 GDAL/OGR version | 3.4.1 PROJ version | 8.2.1 EPSG Registry database version | v10.041 (2021-12-03) GEOS version | 3.10.2-CAPI-1.16.0 SQLite version | 3.37.2 Versão PDAL | 2.3.0 PostgreSQL client version | 13.0 SpatiaLite version | 5.0.1 QWT version | 6.1.3 Versão QScintilla2 | 2.11.5 OS version | Windows 10 Version 2009   |   |   |   Active Python plugins autoSaver | 2.6 ee_plugin | 0.0.4 mapbiomascollection | 1.4 OSMDownloader | 1.0.3 partial_histogram | 1.2 pcraster_tools | 0.2.0 qfieldsync | v4.0.0-beta19 QuickOSM | 2.0.1 quick_map_services | 0.19.28 scpplugin | 2.0 SemiAutomaticClassificationPlugin | 7.10.5 temporalprofiletool | 2.0.3 db_manager | 0.1.20 grassprovider | 2.12.99 MetaSearch | 0.3.5 otbprovider | 2.12.99 processing | 2.12.99 sagaprovider | 2.12.99

Supported QGIS version

New profile

Additional context

No response

gioman commented 2 years ago

Get Landsat 8: Bands 2, 3, 4, 5, 6 and 7.

@HelderSamuelBR instead please attach/link the necessary sample dataset

Plugin SCP: LS8 TOAR

How the SCP plugin is involved in the steps to replicate the issue?

autoSaver 2.6
ee_plugin 0.0.4
mapbiomascollection 1.4
OSMDownloader 1.0.3
partial_histogram 1.2
pcraster_tools 0.2.0
qfieldsync v4.0.0-beta19
QuickOSM 2.0.1
quick_map_services 0.19.28
scpplugin 2.0
SemiAutomaticClassificationPlugin 7.10.5
temporalprofiletool 2.0.3

Despite having checked the checkbox in the issue template, you haven't really tested with a new QGIS profile, please do, and no 3rd party plugins installed.

HelderSamuelBR commented 2 years ago

Get Landsat 8: Bands 2, 3, 4, 5, 6 and 7.

@HelderSamuelBR instead please attach/link the necessary sample dataset

Plugin SCP: LS8 TOAR

How the SCP plugin is involved in the steps to replicate the issue?

autoSaver 2.6 ee_plugin 0.0.4 mapbiomascollection 1.4 OSMDownloader 1.0.3 partial_histogram 1.2 pcraster_tools 0.2.0 qfieldsync v4.0.0-beta19 QuickOSM 2.0.1 quick_map_services 0.19.28 scpplugin 2.0 SemiAutomaticClassificationPlugin 7.10.5 temporalprofiletool 2.0.3

Despite having checked the checkbox in the issue template, you haven't really tested with a new QGIS profile, please do, and no 3rd party plugins installed.

  1. Dataset:
    • LC08_L1TP_218072_20211002_20211013_02_T1 ( https://earthexplorer.usgs.gov ).
    • Reprojection to SRC: Sirgas 2000:31983 UTM 23S
    • Crop raster by extension: 600 km² (Center: 16°41'60.0"S, 43°52'53.4"W)
  2. The layers used in the calculations are relative to the Top of Atmosphere Reflectance (TOAR) which I got from the SCP plugin after processing the dataset.
  3. Yes, I tested it with a new profile. The same results using the original raster (without SCP TOAR correction). Extremely high and irreais values in the resulting rasters. In this case, i used to calculate NDVI (B5 - B4) / (B5 + B4).
gioman commented 2 years ago

2. The layers used in the calculations are relative to the Top of Atmosphere Reflectance (TOAR) which I got from the SCP plugin after processing the dataset.

@HelderSamuelBR then please attach/link the necessary input. We must test with this already processed tif. Do not expect one of us will do it.

LC08_L1TP_218072_20211002_20211013_02_T1 ( https://earthexplorer.usgs.gov ).

link/attach it. You are supposing we may have an account or know of to use earthexplorer, which may not be the case (and and very least is time consuming).

HelderSamuelBR commented 2 years ago
  1. The layers used in the calculations are relative to the Top of Atmosphere Reflectance (TOAR) which I got from the SCP plugin after processing the dataset.

@HelderSamuelBR then please attach/link the necessary input. We must test with this already processed tif. Do not expect one of us will do it.

LC08_L1TP_218072_20211002_20211013_02_T1 ( https://earthexplorer.usgs.gov ).

link/attach it. You are supposing we may have an account or know of to use earthexplorer, which may not be the case (and and very least is time consuming).

Alright, here it is. Link: https://1drv.ms/u/s!AoVzv1X2xEFtnQGYyZoLE7_PmzN4?e=S9iKQR

Equation: (A 0.300)+(B 0.277)+(C 0.233)+(D 0.143)+(E 0.036)+(F 0.012) RTLC08..._B2 > A, RTLC08..._B3 > B [...] RTLC08..._B7 > F The result must be something like 0.01 until > 0.6

HelderSamuelBR commented 2 years ago

A member of QGIS Brazil Community did help me. Only me have this strange results. I've already reinstalled QGIS about 3 times and nothing changes. I'm going to do a downgrade to 3.16 to try to solve.

HelderSamuelBR commented 2 years ago

Results with downgrade and new user.

image

Versão do QGIS | 3.16.16-Hannover | Código da versão do QGIS | f5778a89df -- | -- | -- | -- Compilado sobre Qt | 5.11.2 | Rodando sobre Qt | 5.11.2 Compilado sobre GDAL/OGR | 3.1.4 | Rodando sobre GDAL/OGR | 3.1.4 Compilado sobre GEOS | 3.8.1-CAPI-1.13.3 | Rodando sobre GEOS | 3.8.1-CAPI-1.13.3 Compilado no SQLite | 3.29.0 | Executando contra SQLite | 3.29.0 Versão do cliente PostgreSQL | 11.5 | Versão SpatiaLite | 4.3.0 Versão do QWT | 6.1.3 | Versão QScintilla2 | 2.10.8 Compilado com PROJ | 6.3.2 | Em execução com PROJ | Rel. 6.3.2, May 1st, 2020 Versão SO | Windows 10 (10.0) Ativar complementos python | db_manager; MetaSearch; processing
gioman commented 2 years ago

@HelderSamuelBR I'm investigating.

gioman commented 2 years ago

UPDATE: By some reason, GDAL/r.mapcalc.simple are multiplying the real results by 10.000.

@HelderSamuelBR no, from what I'm seeing is the exact opposite. GDAL and GRASS are doing it right. The (big) issue seems to be in QGIS. The legend values for the rasters you have provided is wrong, they are divided by 10000, and then the QGIS raster calc uses this wrong values for its computation.

Screenshot 2022-03-07 074941

See the stats ( STATISTICS_MAXIMUM and STATISTICS_MINIMUM) at the bottom of the GDALINFO commands

gdalinfo -stats d:\RT_LC08_L1TP_218072_20211002_20211002_01_RT_2021-10-02_B2.tif
Driver: GTiff/GeoTIFF
Files: d:\RT_LC08_L1TP_218072_20211002_20211002_01_RT_2021-10-02_B2.tif
       d:\RT_LC08_L1TP_218072_20211002_20211002_01_RT_2021-10-02_B2.tif.aux.xml
Size is 765, 730
Coordinate System is:
PROJCRS["SIRGAS 2000 / UTM zone 23S",
    BASEGEOGCRS["SIRGAS 2000",
        DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4674]],
    CONVERSION["UTM zone 23S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-45,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Brazil - between 48┬░W and 42┬░W, northern and southern hemispheres, onshore and offshore."],
        BBOX[-33.5,-48,5.13,-42]],
    ID["EPSG",31983]]
Data axis to CRS axis mapping: 1,2
Origin = (610365.000000040046871,8161885.000057072378695)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  610365.000, 8161885.000) ( 43d57'54.90"W, 16d37'22.14"S)
Lower Left  (  610365.000, 8139985.000) ( 43d57'51.06"W, 16d49'14.72"S)
Upper Right (  633315.000, 8161885.000) ( 43d45' 0.41"W, 16d37'17.88"S)
Lower Right (  633315.000, 8139985.000) ( 43d44'55.76"W, 16d49'10.41"S)
Center      (  621840.000, 8150935.000) ( 43d51'25.53"W, 16d43'16.39"S)
Band 1 Block=765x5 Type=UInt16, ColorInterp=Gray
  Min=80.000 Max=2678.000
  Minimum=80.000, Maximum=2678.000, Mean=317.152, StdDev=126.433
  NoData Value=65535
  Offset: 0,   Scale:0.0001
  Metadata:
    STATISTICS_MAXIMUM=2678
    STATISTICS_MEAN=317.15193481959
    STATISTICS_MINIMUM=80
    STATISTICS_STDDEV=126.43304365887
    STATISTICS_VALID_PERCENT=100

GDAL RASTER CALC

gdal_calc.bat --overwrite --calc "A*2" --format GTiff --type Float32 -A D:/RT_LC08_L1TP_218072_20211002_20211002_01_RT_2021-10-02_B2.tif --A_band 1 --outfile D:/gdal_raster_calc.tif

gdalinfo -stats d:\gdal_raster_calc.tif
Driver: GTiff/GeoTIFF
Files: d:\gdal_raster_calc.tif
Size is 765, 730
Coordinate System is:
PROJCRS["SIRGAS 2000 / UTM zone 23S",
    BASEGEOGCRS["SIRGAS 2000",
        DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4674]],
    CONVERSION["UTM zone 23S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-45,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Brazil - between 48┬░W and 42┬░W, northern and southern hemispheres, onshore and offshore."],
        BBOX[-33.5,-48,5.13,-42]],
    ID["EPSG",31983]]
Data axis to CRS axis mapping: 1,2
Origin = (610365.000000040046871,8161885.000057072378695)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  610365.000, 8161885.000) ( 43d57'54.90"W, 16d37'22.14"S)
Lower Left  (  610365.000, 8139985.000) ( 43d57'51.06"W, 16d49'14.72"S)
Upper Right (  633315.000, 8161885.000) ( 43d45' 0.41"W, 16d37'17.88"S)
Lower Right (  633315.000, 8139985.000) ( 43d44'55.76"W, 16d49'10.41"S)
Center      (  621840.000, 8150935.000) ( 43d51'25.53"W, 16d43'16.39"S)
Band 1 Block=765x2 Type=Float32, ColorInterp=Gray
  Minimum=160.000, Maximum=5356.000, Mean=634.304, StdDev=252.866
  NoData Value=3.4028234663852886e+38
  Metadata:
    STATISTICS_MAXIMUM=5356
    STATISTICS_MEAN=634.30386963916
    STATISTICS_MINIMUM=160
    STATISTICS_STDDEV=252.86608731775
    STATISTICS_VALID_PERCENT=100

QGIS RASTER CALC

{ 'CELLSIZE' : 0, 'CRS' : QgsCoordinateReferenceSystem('EPSG:31983'), 'EXPRESSION' : '"RT_LC08_L1TP_218072_20211002_20211002_01_RT_2021-10-02_B2@1" * 2', 'EXTENT' : '610365.000000000,633315.000000000,8139985.000100000,8161885.000100000 [EPSG:31983]', 'LAYERS' : ['C:/Users/giovanni/Desktop/RT_LC08_L1TP_218072_20211002_20211002_01_RT_2021-10-02_B2.tif'], 'OUTPUT' : 'D:/qgis_raster_calc.tif' }

gdalinfo -stats d:\qgis_raster_calc.tif
Driver: GTiff/GeoTIFF
Files: d:\qgis_raster_calc.tif
Size is 765, 730
Coordinate System is:
PROJCRS["SIRGAS 2000 / UTM zone 23S",
    BASEGEOGCRS["SIRGAS 2000",
        DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4674]],
    CONVERSION["UTM zone 23S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-45,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Brazil - between 48┬░W and 42┬░W, northern and southern hemispheres, onshore and offshore."],
        BBOX[-33.5,-48,5.13,-42]],
    ID["EPSG",31983]]
Data axis to CRS axis mapping: 1,2
Origin = (610365.000000000000000,8161885.000099999830127)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  610365.000, 8161885.000) ( 43d57'54.90"W, 16d37'22.14"S)
Lower Left  (  610365.000, 8139985.000) ( 43d57'51.06"W, 16d49'14.72"S)
Upper Right (  633315.000, 8161885.000) ( 43d45' 0.41"W, 16d37'17.88"S)
Lower Right (  633315.000, 8139985.000) ( 43d44'55.76"W, 16d49'10.41"S)
Center      (  621840.000, 8150935.000) ( 43d51'25.53"W, 16d43'16.39"S)
Band 1 Block=765x2 Type=Float32, ColorInterp=Gray
  Min=0.018 Max=0.395
  Minimum=0.016, Maximum=0.536, Mean=0.063, StdDev=0.025
  NoData Value=-3.4028234663852886e+38
  Metadata:
    STATISTICS_MAXIMUM=0.53560000658035
    STATISTICS_MEAN=0.063430387086918
    STATISTICS_MINIMUM=0.016000000759959
    STATISTICS_STDDEV=0.025286608767605
    STATISTICS_VALID_PERCENT=100
agiudiceandrea commented 2 years ago

@gioman see the "scale" metadata displayed in the gdalinfo output of the input raster layers:

Band 1 Block=765x5 Type=UInt16, ColorInterp=Gray
  Min=80.000 Max=2678.000
  Minimum=80.000, Maximum=2678.000, Mean=317.152, StdDev=126.433
  NoData Value=65535
  Offset: 0,   Scale:0.0001

Those raster layers have cells with integer (UInt16) value which needs to be scaled by the 0.0001 scale factor (i.e multiply by 0.0001 or divide by 10000) in order to obtain the actual floating point value.

See the https://gdal.org/programs/gdal_edit.html#cmdoption-scale and https://gdal.org/programs/gdal_translate.html#cmdoption-gdal_translate-unscale

It seems to me the same behaviour occurs in QGIS 3.16.16, 3.22.4 and 3.24.0.

gioman commented 2 years ago

It seems to me the same behaviour occurs in QGIS 3.16.16, 3.22.4 and 3.24.0.

@agiudiceandrea yes.

Offset: 0, Scale:0.0001

@agiudiceandrea I missed that. So as far as I understand:

1) QGIS does the right thing at showing values in the legend and handling the data in further processing

2) QGIS should show somewhere that there is this scaling factor, I don't see it anywhere in QGIS gui or layer properties

3) external tools like GDAL and GRASS need to handle scaling beforehand (i.e. gdal_translate), but still users in QGIS are left with no warning about scaling

Is the above correct?

agiudiceandrea commented 2 years ago

Is the above correct?

Yes, I think so.

For point 2, QGIS displays the raster data type as "Float32" (instead of the raw "UInt16") because it detects the presence of the scale (and offset) and applies it (since QGIS 2.3: https://github.com/qgis/QGIS/pull/1252): https://github.com/qgis/QGIS/blob/667886b3b93fab4ed3ad071ae46681bcd46e9ba6/src/core/providers/gdal/qgsgdalprovider.cpp#L1558-L1602

Maybe it could be useful for the user that QGIS would display also the raw data type and the scale/offset applied.

agiudiceandrea commented 2 years ago

For point 3:

instead,

gioman commented 2 years ago

For point 3:

* when a GRASS processing algorithm (as "r.mapcalc.simple") is used, then the GRASS module `r.in.gdal` is used beforehand in order to convert the input raster layer in a GRASS raster: `r.in.gdal` [doesn't apply](https://trac.osgeo.org/grass/ticket/1640) any scale or offset to the raster

* when the GDAL processing algorithm "Raster calculator" is used, it seems the corresponding `gdal_calc` Python script doesn't take into consideration the scale / offset metadata of the input raster layer

instead,

* when a SAGA processing algorithm (as "Raster calculator") is used then the `io_gdal` SAGA module is used beforehand in order to convert the input raster layer in a SAGA raster: `io_gdal` [do apply](https://github.com/saga-gis/saga-gis/blob/bc940fdbeb34ca5a9403481913a1235b5f9ff199/saga-gis/src/tools/io/io_gdal/gdal_driver.cpp#L965-L971) the scale and offset to the raster

@agiudiceandrea yes, that was also my observation. One more reason to have QGIS really show when a raster has scaling applied.

gioman commented 2 years ago

Closing if favor of https://github.com/qgis/QGIS/issues/47673