Bug report #16646

Small errors in SAGA Raster Calculator output

Added by Pedro Venâncio over 3 years ago. Updated over 3 years ago.

Status:Closed
Priority:Normal
Assignee:Victor Olaya
Category:Processing/SAGA
Affected QGIS version:2.18.9 Regression?:No
Operating System: Easy fix?:No
Pull Request or Patch supplied:No Resolution:
Crashes QGIS or corrupts data:No Copied to github as #:24546

Description

This question was rised in the users mailing list.

We have 3 rasters:

a: Z_top_res.tif
b: t_top_res.tif
c1: topo.tif

The idea is to apply 2 formulas. The first one:

ifelse(c1>0,-(0.0065*c1-16),(0.0045*c1+16))

witch the result gives the raster "c".

And then the second formula:

ifelse(a<500,0,((b-c)/a)*1000)

For instance, taking a pixel of the images with these values (spreadsheet attached):

Z_top_res            a     1202.03510276
t_top_res            b       53.6123428345
topo                 c1    -182.584899902
SAGA Result 1        c       15.1783676147
SAGA Result 2        d       31.9437332153
Expected Result 1            15.1783679504
Expected Result 2            31.9740872222
Difference 1                  0.0303540069
Difference 2                  0.000000335741

Using SAGA GUI with same sample data, there is no error.

Could this be because of some roundings in the input / output (io_gdal)?

topo.tif - Raster "c1" (1.27 KB) Pedro Venâncio, 2017-06-01 12:30 AM

t_top_res.tif - Raster "b" (1.27 KB) Pedro Venâncio, 2017-06-01 12:30 AM

Z_top_res.tif - Raster "a" (2.14 KB) Pedro Venâncio, 2017-06-01 12:30 AM

teste_saga_ifelse.model - Processing Model (7.72 KB) Pedro Venâncio, 2017-06-01 12:33 AM

Errors.ods - Spreadsheet with errors (13.4 KB) Pedro Venâncio, 2017-06-01 12:36 AM

Associated revisions

Revision b49f53ba
Added by Pedro Venâncio over 3 years ago

Changes SAGA io_gdal RESAMPLING method to B-Spline Interpolation, as SAGA default, and add the Resampling Method parameter to SAGA Raster Calculator, as explained in https://issues.qgis.org/issues/16646

Revision 00ca8ccb
Added by Pedro Venâncio over 3 years ago

Changes SAGA io_gdal RESAMPLING method to B-Spline Interpolation, as SAGA default, and add the Resampling Method parameter to SAGA Raster Calculator, as explained in https://issues.qgis.org/issues/16646 - Backport to 2.18

Revision a511ecdd
Added by Alexander Bruy over 3 years ago

Merge pull request #4743 from PedroVenancio/master

[processing] change resampling methods to be like SAGA default (fix #16646)

Revision f98ec626
Added by Alexander Bruy over 3 years ago

Merge pull request #4744 from PedroVenancio/saga_resampling_method

[processing] change resampling methods to be like SAGA default (fix #16646)

History

#1 Updated by Alexander Bruy over 3 years ago

Probably this caused by inputs conversion, as Processing itself does nothing to input rasters, all operations performed by SAGA.

#2 Updated by Pedro Venâncio over 3 years ago

Hi Alexander,

I was investigating and you're right.

On the one hand, SAGA Import Raster Module is being used in Processing with Resampling method 0 (Nearest Neighbour), when the SAGA default for LTR is method 3 (B-Spline Interpolation):

http://www.saga-gis.org/saga_tool_doc/2.3.0/io_gdal_0.html

I think this should be changed, to use the default method: line 64 of SagaAlgorithm230.py

return 'io_gdal 0 -TRANSFORM 1 -RESAMPLING 3 -GRIDS "' + destFilename + '" -FILES "' + source + '"'

With this, the result of the first formula is the expected.

On the other hand, Processing Raster Calculator always considers "Additional Layers" as Grids from different Systems (XGRIDS), which leads to an interpolation, more precisely, a B-Spline Interpolation, as it is the default.

Adding the parameter

ParameterSelection|RESAMPLING|Resampling Method|[0] Nearest Neighbour;[1] Bilinear Interpolation;[2] Bicubic Spline Interpolation;[3] B-Spline Interpolation|3

to the GridCalculator.txt descritpion file, it is possible to choose the interpolation method.

Running Raster Calculator with Nearest Neighbour method, leads to a result for the pixel I've checked before, more close to the expected result: 31.9740869429.

Is there anything that can be done to avoid the interpolation, when it is not necessary? SAGA uses "gx" and "hx" to distinguish between them:

http://www.saga-gis.org/saga_tool_doc/2.3.0/grid_calculus_1.html

But I suspect that Victor have changed that notation to a,b,c,... and to consider all additional layers as grids from different Systems, for some reason.

#3 Updated by Pedro Venâncio over 3 years ago

Hi,

I've been going deeper into this.

In this image, we can see how SAGA GUI deals with rasters in "different systems".

So GRIDS are selected from rasters within the group of
- Resolution: 25
- Size: 11360 x 23200 pixels
- Lower left corner coordinate: -119987.5, -301987.5

XGRIDS come from the group with:
- Resolution: 80
- Size: 3551 x 7251 pixels
- Lower left corner coordinate: -120000, -302000

The RESAMPLING option only show up when there are XGRIDS selected.

I think that what could be done in Processing Raster Calculator is to have two MultipleSelection fields, one for GRIDS and another for XGRIDS.

I've changed Processing Raster Calculator, transforming the GRIDS ParameterRaster, to a ParameterMultipleInput.

GridCalculator.txt

Raster calculator|Grid Calculator
grid_calculus
AllowUnmatching
ParameterMultipleInput|GRIDS|Main input layers (GRIDS)|3|False
ParameterMultipleInput|XGRIDS|Additional layers (XGRIDS)|3|True
ParameterString|FORMULA|Formula|
ParameterBoolean|USE_NODATA|Use NoData|False
ParameterSelection|TYPE|Output Data Type|[0] bit;[1] unsigned 1 byte integer;[2] signed 1 byte integer;[3] unsigned 2 byte integer;[4] signed 2 byte integer;[5] unsigned 4 byte integer;[6] signed 4 byte integer;[7] 4 byte floating point number;[8] 8 byte floating point number|7
OutputRaster|RESULT|Calculated

I made a test, using the three rasters of the sample data (Z_top_res, t_top_res and topo_resultado) in GRIDS.

And I made another, with Z_top_res as GRIDS, and t_top_res and topo_resultado as XGRIDS.

Then I've done the same tests, using saga_cmd and saga_gui.

These are the results for the pixel I've checked before:

So, implementing the GRIDS / XGRIDS logic in Processing, gives the exact same results as SAGA native! As well as using current logic, one GRID and the remaining rasters as XGRIDS.

In conclusion, I think it is clearly confirmed that the differences are related only with RESAMPLING.

So, I'm going to make a pull request with the changes I proposed above:

- Change line 64 of SagaAlgorithm230.py to

return 'io_gdal 0 -TRANSFORM 1 -RESAMPLING 3 -GRIDS "' + destFilename + '" -FILES "' + source + '"'

- Add the parameter to GridCalculator.txt descritpion file

ParameterSelection|RESAMPLING|Resampling Method|[0] Nearest Neighbour;[1] Bilinear Interpolation;[2] Bicubic Spline Interpolation;[3] B-Spline Interpolation|3

because this way, anyone can choose the best interpolator that fits the data used.

I think it is the best thing to do quickly.

But to make things even better, the logic that could be implemented in Processing Raster Calculator was:

1) Select the Main Input Layer (as it is at this moment);
2) Add a new ParameterMultipleInput and stay with two MultipleSelection lists:
    - One ParameterMultipleInput for GRIDS;
    - Another ParameterMultipleInput for XGRIDS;
3) Perform a check and compare the "system" of the Main Input Layer with the remaining rasters loaded in the QGIS, and:
    - For those with the same "system" of Main Input Layer, be listed in the GRIDS MultipleSelection;
    - For the others, be listed on XGRIDS MultipleSelection;
4) Rasters would be referenced in FORMULA as in SAGA ("g" for GRIDS and "h" for XGRIDS);
5) When there was XGRIDS in the formula, RESAMPLING was applied. When there was only GRIDS, there was no RESAMPLING.

Victor, Alexander, Giovanni, do you think this could make sense? Should I open a feature request for this?

Thanks!

#4 Updated by Alexander Bruy over 3 years ago

  • % Done changed from 0 to 100
  • Status changed from Open to Closed

#5 Updated by Giovanni Manghi over 3 years ago

Victor, Alexander, Giovanni, do you think this could make sense? Should I open a feature request for this?

Hi Pedro, yes,make a feature request for this please! cheers!

Also available in: Atom PDF