How to use SAGA Raster processing in processing scripts (QGIS3)

Question

I want to make a processing script which uses the saga:inversedistanceweightedinterpolation algorithm (Interpolating a point layer - magnetometry measurements - to a raster. The SAGA runs so far but the result could not be loaded by QGIS. I have the following script:

from PyQt5.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
 QgsProcessingAlgorithm,
 QgsProcessingParameterFeatureSource,
 QgsProcessingParameterRasterDestination)
import processing
class ExampleProcessingAlgorithm(QgsProcessingAlgorithm):
 INPUT = 'INPUT'
 OUTPUT = 'OUTPUT'
 def tr(self, string):
 return QCoreApplication.translate('Processing', string)
 def createInstance(self):
 return ExampleProcessingAlgorithm()
 def name(self):
 return 'test_IDW'
 def displayName(self):
 return self.tr('Test IDW')
 def group(self):
 return self.tr('Example scripts')
 def groupId(self):
 return 'examplescripts'
 def shortHelpString(self):
 return self.tr("Example algorithm short description")
 def initAlgorithm(self, config=None):
 self.addParameter(
 QgsProcessingParameterFeatureSource(
 self.INPUT,
 self.tr('Input layer'),
 [QgsProcessing.TypeVectorPoint]))
 self.addParameter(
 QgsProcessingParameterRasterDestination(
 self.OUTPUT,
 self.tr("Output Raster"),
 None, False))
 def processAlgorithm(self, parameters, context, feedback):
 output = processing.run("saga:inversedistanceweightedinterpolation", {
 'DW_BANDWIDTH' : 1, 
 'DW_IDW_OFFSET' : False, 
 'DW_IDW_POWER' : 1, 
 'DW_WEIGHTING' : 1, 
 'FIELD' : 'READING', 
 'SEARCH_DIRECTION' : 0, 
 'SEARCH_POINTS_ALL' : 0, 
 'SEARCH_POINTS_MAX' : 20, 
 'SEARCH_POINTS_MIN' : -1, 
 'SEARCH_RADIUS' : 1, 
 'SEARCH_RANGE' : 0, 
 'SHAPES' : parameters['INPUT'], 
 'TARGET_DEFINITION' : 0, 
 'TARGET_OUT_GRID' : parameters['OUTPUT'], #'memory:', 
 'TARGET_TEMPLATE' : None, 
 'TARGET_USER_FITS' : 0, 
 'TARGET_USER_SIZE' : 0.25, 
 'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX' : None 
 }, context=context, feedback=feedback)['TARGET_OUT_GRID']
 return {self.OUTPUT: output}

I get the following error after execution:
Loading resulting layers
The following layers were not correctly generated.

• C:/Users/Stefan/AppData/Local/Temp/processing_03bc47ec83e54650bd96292f861633b5/11bb8145ca214423b41e7c1c33c3965c/OUTPUT.tif

You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.
the Log Message Panel shows: 2019年01月30日T18:06:45 INFO Task complete : Executing "Test IDW"

Answer 1

SAGA algorithms cannot natively write to TIFF formats -- they only support SAGA's native sdat format. Try using a file path with a .sdat extension.

• Paulo Raposo: True. You can then convert back to GeoTiff using GDAL, either gdal_translate on the command line (see gdal.org/gdal_translate.html), or by loading the .sdat/.sgrd in a Python script using the gdal library and saving back out to a .tif (see pcjericks.github.io/py-gdalogr-cookbook/index.html). SAGA has native import/export too, but I've not had good luck with those tools.

Answer 2

You are right - SAGA runs fine, but the output is not what QGIS processing expects.

Processing insists on TIFF for raster output, so parameters['OUTPUT'] will be a TIFF file reference. Saga does not mind getting a TIFF file reference for its output - it just changes it to an sdat file reference by switching the extension in the reference string.

So the SAGA algorithm will run fine, but produce sdat format output (you will probably find that you can open C:/Users/Stefan/AppData/Local/Temp/processing_03bc47ec83e54650bd96292f861633b5/11bb8145ca214423b41e7c1c33c3965c/OUTPUT.sdat). But since SAGA has changed the file reference, there will be no C:/Users/Stefan/AppData/Local/Temp/processing_03bc47ec83e54650bd96292f861633b5/11bb8145ca214423b41e7c1c33c3965c/OUTPUT.tif - resulting in the error message you observe.

If you would like the result to be loaded once the script has been run, you can use GDAL-translate on the result of the SAGA algorithm:

processing.run('gdal:translate', {'INPUT': output['TARGET_OUT_GRID'],
 'OUTPUT': parameters['OUTPUT']})

But first you will have to fix the value that you provide for 'TARGET_OUT_GRID' in the SAGA algorithm:

newtarget = parameters['OUTPUT'].replace(".tif", ".sdat")

and then use newtarget instead of parameters['OUTPUT'] in the SAGA algorithm.

It is a bit of a hack, but it should work.