Wrong values in RasterMath?

[Agustin Lobo]

I’ve done the following operation in RasterMath:

(LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@5 - LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@4) / (LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@5 + LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@4)

Surprisingly, there are many ndvi values >>1 wherever the value should be <0 (rest of values are correct)

Log:

Python command:

processing.run('enmapbox:RasterMath', dict(code='# LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY := QgsRasterLayer("/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY.vrt")\n(LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@5 - LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@4) / (LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@5 + LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@4)', monolithic=False, outputRaster='/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/LC08_L2SP_202034_20211205_20211215_02_T1_SR_NDVI.TIF'))

Console command:

qgis_process run enmapbox:RasterMath --code="# LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY := QgsRasterLayer("/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY.vrt")\n(LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@5 - LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@4) / (LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@5 + LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@4)" monolithic=False outputRaster=/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/LC08_L2SP_202034_20211205_20211215_02_T1_SR_NDVI.TIF

Create Raster 7571x7711x1 -co INTERLEAVE=BAND COMPRESS=LZW TILED=YES BIGTIFF=YES /media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/LC08_L2SP_202034_20211205_20211215_02_T1_SR_NDVI.TIF
Execution completed in 9.11 seconds
Results:
{'outputRaster': '/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/LC08_L2SP_202034_20211205_20211215_02_T1_SR_NDVI.TIF'}
Execution completed in 13.97 seconds
Results:
{'outputRaster': '/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/LC08_L2SP_202034_20211205_20211215_02_T1_SR_NDVI.TIF'}

Loading resulting layers
Algorithm 'Raster math' finished

‌

[Andreas Janz]

I guess for those pixels, the input values are not strictly between 0 and 1 (or 0 and 10000). In those cases, a normalized difference can fall outside the 0 to 1 range.

‌

[Agustin Lobo]

The values must be in the -1, +1 range. In that pixel the actual correct value is

(7909-10182)/(7909+10182) = -0.1256426

Raster Calculator in QGIS calculates the index right.

[Andreas Janz]

Hi Agus,

what is the data type of your raster? Is it UINT16? This could explain the behaviour: 7909-10182 would result in an underflow when using Numpy arrays.

In this case, you need to cast the data before you calculate the NDVI:

array = array.astype(np.float32)

‌

[Agustin Lobo]

Yes, the data (original Landsat 8 images) are UINT16.

I presume this is going to be an endless source of ugly errors (the user is not aware), as there is a lot of RS imagery out there in UINT16 . Should you not be creating float32 by default? Then let the user cast the result if desired.

Anyway, by now I understand I could modify the code by adding ndvi.astype(np.float32) such as:

# find bands
red = {raster1}@4
nir = {raster1}@5

# calculate NDVI
ndvi = (nir - red) / (nir + red)

# mask no data region
noDataValue = 0
ndvi[~{raster1}Mask@4] = noDataValue
ndvi[~{raster1}Mask@5] = noDataValue

# set no data value and band name
ndvi.setNoDataValue(noDataValue)
ndvi.setBandName('NDVI', bandNo=1)
ndvi.astype(np.float32)

# clean up temp variables
del red, nir

‌

The problem is that I cannot even reproduce the code as I had it (without line 16) now that I use the dev version. I get:

Create Raster [7571x7711x1](Float64) -co INTERLEAVE=BAND COMPRESS=LZW TILED=YES BIGTIFF=YES /media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndvi.tif

Execution completed in 9.7 seconds

Results:
{'outputRaster': None, 'ndvi': '/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndvi.tif'}

Execution completed in 15.12 seconds

Results:

{'ndvi': '/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndvi.tif',
'outputRaster': None}

‌

Loading resulting layers

The following layers were not correctly generated.

/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndviv2.tif

You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

and in the Log Messages I just get Warnings related to this issue (it is not clear what the problem is):

2022-01-04T18:56:52     WARNING    warning:/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/pyqtgraph/functions.py:41: DeprecationWarning:

Flags not at the start of the expression '(?P[+-]?((((' (truncated)

traceback: File "", line 1, in
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 443, in startPlugin
              if not _startPlugin(packageName):
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 423, in _startPlugin
              plugins[packageName] = package.classFactory(iface)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/__init__.py", line 28, in classFactory
              from .profileplugin import ProfilePlugin
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/profileplugin.py", line 41, in
              from .tools.profiletool_core import ProfileToolCore
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/tools/profiletool_core.py", line 48, in
              from ..ui.ptdockwidget import PTDockWidget
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/ui/ptdockwidget.py", line 40, in
              from ..tools.plottingtool import *
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/tools/plottingtool.py", line 44, in
              from .. import pyqtgraph as pg
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1058, in _handle_fromlist
              File "", line 228, in _call_with_frames_removed
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/pyqtgraph/__init__.py", line 202, in
              from .graphicsItems.VTickGroup import *
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/pyqtgraph/graphicsItems/VTickGroup.py", line 7, in
              from .. import functions as fn
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1058, in _handle_fromlist
              File "", line 228, in _call_with_frames_removed
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/pyqtgraph/functions.py", line 41, in
              FLOAT_REGEX = re.compile(r'(?P[+-]?((((\d+(\.\d*)?)|(\d*\.\d+))([eE][+-]?\d+)?)|((?i)(nan)|(inf))))\s*((?P[u' + SI_PREFIXES + r']?)(?P\w.*))?$')
              File "/usr/lib/python3.9/re.py", line 252, in compile
              return _compile(pattern, flags)
              File "/usr/lib/python3.9/re.py", line 304, in _compile
              p = sre_compile.compile(pattern, flags)
              File "/usr/lib/python3.9/sre_compile.py", line 764, in compile
              p = sre_parse.parse(p, flags)
              File "/usr/lib/python3.9/sre_parse.py", line 948, in parse
              p = _parse_sub(source, state, flags & SRE_FLAG_VERBOSE, 0)
              File "/usr/lib/python3.9/sre_parse.py", line 443, in _parse_sub
              itemsappend(_parse(source, state, verbose, nested + 1,
              File "/usr/lib/python3.9/sre_parse.py", line 834, in _parse
              p = _parse_sub(source, state, sub_verbose, nested + 1)
              File "/usr/lib/python3.9/sre_parse.py", line 443, in _parse_sub
              itemsappend(_parse(source, state, verbose, nested + 1,
              File "/usr/lib/python3.9/sre_parse.py", line 834, in _parse
              p = _parse_sub(source, state, sub_verbose, nested + 1)
              File "/usr/lib/python3.9/sre_parse.py", line 443, in _parse_sub
              itemsappend(_parse(source, state, verbose, nested + 1,
              File "/usr/lib/python3.9/sre_parse.py", line 834, in _parse
              p = _parse_sub(source, state, sub_verbose, nested + 1)
              File "/usr/lib/python3.9/sre_parse.py", line 443, in _parse_sub
              itemsappend(_parse(source, state, verbose, nested + 1,
              File "/usr/lib/python3.9/sre_parse.py", line 809, in _parse
              warnings.warn(

2022-01-04T18:56:54     WARNING    warning:/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/geetimeseriesexplorer/site/pyqtgraph/graphicsItems/PlotItem/PlotItem.py:340: DeprecationWarning:

an integer is required (got type numpy.float64). Implicit conversion to integers using __int__ is deprecated, and may be removed in a future version of Python.

traceback: File "", line 1, in
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 448, in startPlugin
              plugins[packageName].initGui()
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/geetimeseriesexplorer/plugin.py", line 19, in initGui
              self.ui = DockWidget(iface=self.iface, parent=self.iface.parent())
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/geetimeseriesexplorer/dockwidget.py", line 130, in __init__
              self._initGui()
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/geetimeseriesexplorer/dockwidget.py", line 174, in _initGui
              self.mPlotWidget.showGrid(x=True, y=True, alpha=1.)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/geetimeseriesexplorer/site/pyqtgraph/graphicsItems/PlotItem/PlotItem.py", line 340, in showGrid
              self.ctrl.gridAlphaSlider.setValue(v)
              File "/usr/lib/python3.9/warnings.py", line 109, in _showwarnmsg
              sw(msg.message, msg.category, msg.filename, msg.lineno,

2022-01-04T18:56:55     WARNING    warning:/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/ee/data.py:230: ResourceWarning:

unclosed file

traceback: File "", line 1, in
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 443, in startPlugin
              if not _startPlugin(packageName):
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 423, in _startPlugin
              plugins[packageName] = package.classFactory(iface)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/qgis_gee_data_catalog/__init__.py", line 33, in classFactory
              from .qgis_gee_data_catalog import GeeDataCatalog
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/__init__.py", line 32, in __wrapping_ee_import__
              _module_ = __builtin_import__(name, *args, **kwargs)
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/qgis_gee_data_catalog/qgis_gee_data_catalog.py", line 26, in
              from .ee_interface import (add_ee_image_layer, download_ee_image_layer,
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/__init__.py", line 32, in __wrapping_ee_import__
              _module_ = __builtin_import__(name, *args, **kwargs)
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/qgis_gee_data_catalog/ee_interface.py", line 6, in
              import ee
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/__init__.py", line 36, in __wrapping_ee_import__
              _module_.Initialize()
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/ee/__init__.py", line 114, in Initialize
              credentials = data.get_persistent_credentials()
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/ee/data.py", line 230, in get_persistent_credentials
              tokens = json.load(open(oauth.get_credentials_path()))
              File "/usr/lib/python3.9/warnings.py", line 109, in _showwarnmsg
              sw(msg.message, msg.category, msg.filename, msg.lineno,

2022-01-04T18:56:55     WARNING    warning:/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/httplib2shim/__init__.py:45: DeprecationWarning:

Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3, and in 3.10 it will stop working

traceback: File "", line 1, in
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 443, in startPlugin
              if not _startPlugin(packageName):
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 423, in _startPlugin
              plugins[packageName] = package.classFactory(iface)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/qgis_gee_data_catalog/__init__.py", line 33, in classFactory
              from .qgis_gee_data_catalog import GeeDataCatalog
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/__init__.py", line 32, in __wrapping_ee_import__
              _module_ = __builtin_import__(name, *args, **kwargs)
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/qgis_gee_data_catalog/qgis_gee_data_catalog.py", line 26, in
              from .ee_interface import (add_ee_image_layer, download_ee_image_layer,
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/__init__.py", line 32, in __wrapping_ee_import__
              _module_ = __builtin_import__(name, *args, **kwargs)
              File "/usr/lib/python3/dist-packages/qgis/utils.py", line 888, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1007, in _find_and_load
              File "", line 986, in _find_and_load_unlocked
              File "", line 680, in _load_unlocked
              File "", line 855, in exec_module
              File "", line 228, in _call_with_frames_removed
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/qgis_gee_data_catalog/ee_interface.py", line 6, in
              import ee
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/__init__.py", line 36, in __wrapping_ee_import__
              _module_.Initialize()
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/ee/__init__.py", line 115, in Initialize
              data.initialize(
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/ee/data.py", line 210, in initialize
              _install_cloud_api_resource()
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/ee/data.py", line 280, in _install_cloud_api_resource
              _cloud_api_resource = _cloud_api_utils.build_cloud_resource(
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/ee/_cloud_api_utils.py", line 133, in build_cloud_resource
              http_transport = httplib2.Http(timeout=timeout)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/httplib2shim/__init__.py", line 117, in __init__
              pool = self._make_pool(proxy_info=proxy_info)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/ee_plugin/extlibs_linux/httplib2shim/__init__.py", line 45, in _default_make_pool
              if isinstance(proxy_info, collections.Callable):
              File "/usr/lib/python3.9/collections/__init__.py", line 62, in __getattr__
              warnings.warn("Using or importing the ABCs from 'collections' instead "
              File "/usr/lib/python3.9/warnings.py", line 109, in _showwarnmsg
              sw(msg.message, msg.category, msg.filename, msg.lineno,

2022-01-04T19:00:17     WARNING    warning::4: RuntimeWarning:

invalid value encountered in true_divide

traceback: File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/site-packages/typeguard/__init__.py", line 903, in wrapper
              retval = func(*args, **kwargs)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/enmapboxprocessing/algorithm/rastermathalgorithm/rastermathalgorithm.py", line 264, in processAlgorithm
              writers = self.makeWriter(code, filename, grid, readers, readers2, feedback)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/site-packages/typeguard/__init__.py", line 903, in wrapper
              retval = func(*args, **kwargs)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/enmapboxprocessing/algorithm/rastermathalgorithm/rastermathalgorithm.py", line 321, in makeWriter
              results = self.processBlock(code, block, readers, readers2, writers, 0, feedback, dryRun=True)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/site-packages/typeguard/__init__.py", line 903, in wrapper
              retval = func(*args, **kwargs)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/enmapboxprocessing/algorithm/rastermathalgorithm/rastermathalgorithm.py", line 513, in processBlock
              exec(code, namespace)
              File "", line 4, in
              File "/usr/lib/python3.9/warnings.py", line 109, in _showwarnmsg
              sw(msg.message, msg.category, msg.filename, msg.lineno,

2022-01-04T19:00:26     WARNING    warning::4: RuntimeWarning:

divide by zero encountered in true_divide

traceback: File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/site-packages/typeguard/__init__.py", line 903, in wrapper
              retval = func(*args, **kwargs)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/enmapboxprocessing/algorithm/rastermathalgorithm/rastermathalgorithm.py", line 284, in processAlgorithm
              results = self.processBlock(
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/site-packages/typeguard/__init__.py", line 903, in wrapper
              retval = func(*args, **kwargs)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/enmapboxprocessing/algorithm/rastermathalgorithm/rastermathalgorithm.py", line 513, in processBlock
              exec(code, namespace)
              File "", line 4, in
              File "/usr/lib/python3.9/warnings.py", line 109, in _showwarnmsg
              sw(msg.message, msg.category, msg.filename, msg.lineno,

2022-01-04T19:00:26     WARNING    warning::4: RuntimeWarning:

invalid value encountered in true_divide

traceback: File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/site-packages/typeguard/__init__.py", line 903, in wrapper
              retval = func(*args, **kwargs)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/enmapboxprocessing/algorithm/rastermathalgorithm/rastermathalgorithm.py", line 284, in processAlgorithm
              results = self.processBlock(
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/site-packages/typeguard/__init__.py", line 903, in wrapper
              retval = func(*args, **kwargs)
              File "/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/enmapboxplugin/enmapboxprocessing/algorithm/rastermathalgorithm/rastermathalgorithm.py", line 513, in processBlock
              exec(code, namespace)
              File "", line 4, in
              File "/usr/lib/python3.9/warnings.py", line 109, in _showwarnmsg
              sw(msg.message, msg.category, msg.filename, msg.lineno,

‌

[Andreas Janz]

The problem is that I cannot even reproduce the code as I had it (without line 16) now that I use the dev version. I get:

Hi Agus,

I’m a bit confused by this filename:
/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndviv2.tif

Please show a screenshot of the RasterMath dialog with all the inputs.

And also copy&paste the code snippet.

‌

[Andreas Janz]

All those warnings are from another plugin (profiletool ?):

WARNING    warning:/home/alobo/.local/share/QGIS/QGIS3/profiles/default/python/plugins/profiletool/pyqtgraph/functions.py:41:

[Andreas Janz]

Yes, the data (original Landsat 8 images) are UINT16.

I presume this is going to be an endless source of ugly errors (the user is not aware), as there is a lot of RS imagery out there in UINT16 . Should you not be creating float32 by default? Then let the user cast the result if desired.

Yes, let’s have a checkbox “Use 32-bit floating-point inputs” that is checked by default. See follow-up issue #908.

[Agustin Lobo]

I’m a bit confused by this filename:
/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndviv2.tif

Please show a screenshot of the RasterMath dialog with all the inputs.

And also copy&paste the code snippet.

The original code (without casting) and once the filenames have been automatically included:

# find bands
red = LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@4
nir = LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@5

# calculate NDVI
ndvi = (nir - red) / (nir + red)

# mask no data region
noDataValue = 0
ndvi[~LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAYMask@4] = noDataValue
ndvi[~LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAYMask@5] = noDataValue

# set no data value and band name
ndvi.setNoDataValue(noDataValue)
ndvi.setBandName('NDVI', bandNo=1)

# clean up temp variables
del red, nir

The Log:

QGIS version: 3.22.2-Białowieża
QGIS code revision: 1601ec46d0
Qt version: 5.15.2
Python version: 3.9.5
GDAL version: 3.2.2
GEOS version: 3.9.0-CAPI-1.16.2
PROJ version: Rel. 7.2.1, January 1st, 2021
PDAL version: 2.2.0 (git-version: Release)
Algorithm started at: 2022-01-05T10:38:09
Algorithm 'Raster math' starting…
Input parameters:
{ 'R1' : None, 'R10' : None, 'R2' : None, 'R3' : None, 'R4' : None, 'R5' : None, 'R6' : None, 'R7' : None, 'R8' : None, 'R9' : None, 'RS' : None, 'V1' : None, 'V10' : None, 'V2' : None, 'V3' : None, 'V4' : None, 'V5' : None, 'V6' : None, 'V7' : None, 'V8' : None, 'V9' : None, 'code' : '# LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY := QgsRasterLayer("/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY.vrt")\n# find bands\nred = LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@4\nnir = LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAY@5\n\n# calculate NDVI\nndvi = (nir - red) / (nir + red)\n\n# mask no data region\nnoDataValue = 0\nndvi[~LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAYMask@4] = noDataValue\nndvi[~LC08_L2SP_202034_20211205_20211215_02_T1_SR_DAYMask@5] = noDataValue\n\n# set no data value and band name\nndvi.setNoDataValue(noDataValue)\nndvi.setBandName('NDVI', bandNo=1)\n\n# clean up temp variables\ndel red, nir', 'grid' : None, 'monolithic' : False, 'outputRaster' : '/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndviv2.tif', 'overlap' : None }

Python command:

Console command:

Create Raster 7571x7711x1 -co INTERLEAVE=BAND COMPRESS=LZW TILED=YES BIGTIFF=YES /media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndvi.tif
Execution completed in 29.47 seconds
Results:
{'outputRaster': None, 'ndvi': '/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndvi.tif'}
Execution completed in 37.68 seconds
Results:
{'ndvi': '/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndvi.tif',
'outputRaster': None}

‌

[Andreas Janz]

And what is the problem now? Is the ndvi.tif not created?
/media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndvi.tif

I have no problems on the testdata:

‌

You can send me your Landsat raster and I’ll check it, maybe it is a data specific problem.

[Agustin Lobo]

Sorry, I missed to copy the last few lines of the Log

Loading resulting layers
The following layers were not correctly generated.
• /media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndviv2.tif
You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

Why should it write “ndvi.tif” if I selected “ndviv2.tif” as output?

‌

[Andreas Janz]

Ah ok, now I see the problem.

The logic is as follows: you can calculate as many outputs as you want in your code snippet. You do that by specifying new variables, e.g. red, green, ndvi:

All temp variables you don’t want to write, need to be deleted, e.g. red, nir:

‌

Now, in case you only have a single output (e.g. ndvi), you can specify where to store the result using the Output raster layer parameter (this is the default output). Here the file basename (without extension) need to match the variable name. In your case they don’t match: ndvi vs. ndviv2.

You get the error, because you haven’t specified a variable named ndvi2.

But in the same folder, you should find a ndvi.tif file for the “additional” ndvi output.

[Agustin Lobo]

Here the file basename (without extension) need to match the variable name.

ok, clarified.

But I think this is quite weird for the user: names of variables in the code should not condition output file names. Note input file names do not condition the names of the variables (you name “red” a variable which input is named “enmap_berlin@655nm” ), why putting restrictions to output filenames?

Also, while being able to clean up variables is good for memory efficiency, the default should be that all variables are cleaned after running, and result variables should be specifically saved, something like eg.

write (ndvi, file=“ /media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndviv2.bil”, format=”ENVIf”, size=”UINT16”, interleave=”BIL”)

Could this be contemplated for future versions?

[Andreas Janz]

But I think this is quite weird for the user: names of variables in the code should not condition output file names. Note input file names do not condition the names of the variables (you name “red” a variable which input is named “enmap_berlin@655nm” ), why putting restrictions to output filenames?

No, not really. When the user code is executed. The data is already mapped to a variable called “enmap_berlin”. The user assigns it to another variable called “red”.

Also, while being able to clean up variables is good for memory efficiency, the default should be that all variables are cleaned after running, and result variables should be specifically saved, something like eg.

write (ndvi, file=“ /media/alobo/LaCieNTFS2T/Ignacio2019/RioTinto/LC08_DAY/LC08_L2SP_202034_20211205_20211215_02_T1/ndviv2.bil”, format=”ENVIf”, size=”UINT16”, interleave=”BIL”)

Could this be contemplated for future versions?

We could change that. I will discuss this with the team.

‌