Dealing with vectors in processing script

Question

So I am writing a processing script with inputs:

• multiple raster layers
• a text file

the output is vector file (multipolygon)

I initialize a vector file called mergedLayer

mergedLayer = QgsVectorLayer('MultiPolygon','temp','memory')
print(mergedLayer)

gives <QgsVectorLayer: 'temp' (memory)>

the steps of the algorithm are in summary

for each raster layer

• raster calculator with a formula that populates the raster with n when informed and 0 elsewhere
• Polygonize (raster to vector)
• merge the polygonized output with mergedLayer
• save in mergedLayer the output of the raster to polygon algorithm before following iteration

I am blocked at the merge step

alg_params = {
 'LAYERS': [outputs['PolygonizeRasterToVector']['OUTPUT'],mergedLayer],
 'CRS':parameters['crs'], 
 'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
 }
outputs['MergeVectorLayers'] = processing.run('native:mergevectorlayers', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
mergedLayer=outputs['mergevectorlayers']['OUTPUT']
 

the scripts runs with no error but there is no output

when I print the output of the Polygonize algorithm, I get a path to a geopackage

print(outputs['PolygonizeRasterToVector']['OUTPUT'])

/tmp/processing_SMBizu/59eb9bda47b04b318c9e0ff885ba6ac6/OUTPUT.gpkg

but when I look at the mergedLayer or output of the merge

print(mergedLayer)

I get a string with a temporary layer Merged_1785d2cf_af36_466c_8ea6_12058b411d54

I am guessing that my inputs are not compatible and that is why the merge is not working

How do I get to use both vectors in the same format?

---

Edit after Ben W answer 8/02/2025

I have tried Ben W's code and this one modified and with both I get the same error. When putting a shapefile as ouptput:

Feature could not be written to /media/geotyr/geotyr_data/PROJETS/SOT/GIS/PC_2022/RASTER/t9.shp: Feature creation error (OGR error: Coordinates with non-finite values are not allowed)
Could not write feature into OUTPUT

When I put try a geopackage file, the file is created, but none of the features are visible. When I try to compute areas, they are all "nan".

So I am getting close, but I do not understand what is happening with the output file

So here is the adapted script modified with most of @BenW’s suggestions.

from qgis.core import (QgsProcessing,QgsProcessingAlgorithm,QgsProcessingMultiStepFeedback,QgsProcessingParameterRasterLayer,QgsProcessingParameterVectorLayer,QgsProcessingParameterCrs,QgsProcessingParameterRasterDestination,QgsProcessingParameterBoolean,QgsProcessingParameterMultipleLayers,QgsCoordinateReferenceSystem,QgsProcessingParameterVectorDestination,QgsVectorLayer,QgsRasterLayer)
from PyQt5.QtCore import QFileInfo
import processing
class VectoritzarAProfunditatsvstack(QgsProcessingAlgorithm):
 
 OUTPUT = 'OUTPUT'
 def initAlgorithm(self, config=None):
 
 self.addParameter(QgsProcessingParameterMultipleLayers('threshold_layers', 'threshold layers', layerType=QgsProcessing.TypeRaster, defaultValue=None))
 self.addParameter(QgsProcessingParameterCrs('crs', 'CRS', defaultValue='EPSG:25831'))
 self.addParameter(QgsProcessingParameterVectorDestination(self.OUTPUT,'Capa vectorial final',QgsProcessing.TypeVectorPolygon))
 
 def processAlgorithm(self, parameters, context, model_feedback):
 
 results = {}
 outputs = {}
 # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
 # overall progress through the model
 steps = (len(parameters['threshold_layers']))*2+1
 feedback = QgsProcessingMultiStepFeedback(steps, model_feedback)
 step = 1
 
 
 
 # Raster counter
 raster_n=1
 #####################################################################################
 # Loop on raster files
 #####################################################################################
 
 for rasterFile in parameters['threshold_layers']:
 # QFileInfo for all the info on the file
 fileInfo = QFileInfo(rasterFile)
 # EExtract the name
 rlayerName = fileInfo.baseName()
 # create the raster layer
 rlayer=QgsRasterLayer(rasterFile,rlayerName)
 #########################################
 # Raster calculator
 #########################################
 
 # create the formula
 # for defined pixel we assign the raster_n value
 
 formula = "(A>-9999999)*{}".format(raster_n)
 r_calc_params = {
 'BAND_A': 1,
 'FORMULA': formula,
 'INPUT_A': rlayerName,
 'RTYPE': 5,
 'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
 }
 
 feedback.setCurrentStep(step)
 step+=1
 if feedback.isCanceled():
 return {}
 
 
 # results of calculator
 outputs['RasterCalculator'] = processing.run('gdal:rastercalculator', r_calc_params, context=context, feedback=feedback, is_child_algorithm=True)
 
 #########################################
 # Polygonize (raster to vector)
 #########################################
 vectorize_params = {
 'BAND': 1,
 'EIGHT_CONNECTEDNESS': False,
 'FIELD': 'Nivell',
 'INPUT': outputs['RasterCalculator']['OUTPUT'],
 'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
 }
 # resultats of vectorize
 feedback.setCurrentStep(step)
 step+=1
 outputs[f'vectorized_{raster_n}']=processing.run("gdal:polygonize", vectorize_params, feedback=feedback, context=context, is_child_algorithm=True)
 results[f'vectorized_{raster_n}']=outputs[f'vectorized_{raster_n}']['OUTPUT']
 raster_n+=1
 
 
 ########################################
 # Merge 
 ########################################
 
 # Build a list of the vectorized outputs in the results dictionary
 v_lyrs = [v for k, v in results.items() if 'vectorized' in k]
 
 # Parametres del Merge
 merge_params = {
 'LAYERS':v_lyrs,
 'CRS':parameters['crs'],
 'OUTPUT':parameters[self.OUTPUT]
 }
 
 feedback.setCurrentStep(step)
 step+=1
 outputs['merged']=processing.run("native:mergevectorlayers", merge_params, feedback=feedback, context=context, is_child_algorithm=True)
 results['merged']=outputs['merged']['OUTPUT']
 
 return results
 def name(self):
 return 'Vectoritzar a profunditats vstack'
 def displayName(self):
 return 'Vectoritzar a profunditats vstack'
 def group(self):
 return 'Sot'
 def groupId(self):
 return 'Sot'
 def createInstance(self):
 return VectoritzarAProfunditatsvstack()

• Ben W: Ok, this seems more a problem with your data or one which arises because of your data. I tested my answer (in QGIS 3.34) and it works fine with with temporary, geopackage and shapefile outputs. For testing purposes, my input rasters were 2 srtm DEMs. Does your workflow work when run manually via the gui? A tip for getting params & syntax right is to run an alg via gui then go to Processing History (clock icon) and click on the entry to see the processing.run call which was generated. I am not sure about your raster calc params or formula- what version of QGIS are you using?
• Ben W: Because, if I understand what you are trying to do, the formula would be something like: '({}@1>-9999999)*{}@1'.format(rlayerName, rlayerName). For further assistance would you be able to share a couple of your rasters for testing? And regarding the ogr error, see: gis.stackexchange.com/questions/426998/…
• geotyr: so, I posted a simplified script because this port is out of the scope of the question, what I actually do is extract from the raster name 1, 2, 3 last characters depending on the length of the name and store them in a variable called level which is the one used in the formula. Regarding the Ogr error...I did see that thread, but did no read all post...I now see it is because of the repeated fid of the geopackage temporary file. I did some tests...but I still I am not able to get the output vector layer stored as a string...I will pause it for now
• geotyr: here is a link with a couple of rasters, my current script (sorry, all comments are in catalan) and a text file that I also need in the script owncloud.dedikam.com/index.php/s/fk3ssZ5z6SDeE4Q (I was not able to find a way to officially add data to my question)
• geotyr: so I managed to edit fid in the new version of the script, but I still get the same Ogr error. I am guessing it is because I am doing it after the merge, so I will try doing it before and see if it solves anything. Thanks again for all the inputs Ben W.

Answer 1

Instead of doing an iterative, interim merge at every step, I would suggest a slightly different approach.

Outside the loop, at the completion of the raster calculation & polygonizing steps, retrieve the outputs of all the "gdal:polygonize" calls and use them all as inputs to a single merge operation.

If you also use a QgsProcessingParameterVectorDestination as the output of your parent algorithm, this has the added advantage that you can feed this to the 'OUTPUT' param of the "native:mergevectorlayers" child algorithm so that the output layer is properly handled. Otherwise you are left with the memory layer output of the final merge call which you would then need to handle yourself- perhaps via context.addLayerToLoadOnCompletion() or add another child algorithm e.g. "native:savefeatures".

Here is a basic example:

from qgis.PyQt.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsProcessingParameterMultipleLayers,
 QgsProcessing, QgsProcessingAlgorithm,
 QgsProcessingParameterVectorDestination,
 QgsProcessingMultiStepFeedback)
import processing
class ExAlgo(QgsProcessingAlgorithm):
 INPUT_RASTERS = 'INPUT_RASTERS'
 OUTPUT = 'OUTPUT'
 
 def __init__(self):
 super().__init__()
 
 def name(self):
 return "calc-vectorize-merge"
 
 def tr(self, text):
 return QCoreApplication.translate("calc-vectorize-merge", text)
 
 def displayName(self):
 return self.tr("Calc Vectorize Merge")
 
 def group(self):
 return self.tr("Examples")
 
 def groupId(self):
 return "examples"
 
 def shortHelpString(self):
 return self.tr("Raster calculation, vectorize and merge pipeline")
 
 def helpUrl(self):
 return "https://qgis.org"
 
 def createInstance(self):
 return type(self)()
 
 def initAlgorithm(self, config=None):
 self.addParameter(QgsProcessingParameterMultipleLayers(
 self.INPUT_RASTERS,
 self.tr("Input rasters"),
 QgsProcessing.TypeRaster))
 
 self.addParameter(QgsProcessingParameterVectorDestination(
 self.OUTPUT,
 self.tr("Output layer"),
 QgsProcessing.TypeVectorPolygon))
 
 def processAlgorithm(self, parameters, context, model_feedback):
 outputs = {}
 results = {}
 input_rasters = self.parameterAsLayerList(parameters, self.INPUT_RASTERS, context)
 
 steps = (len(input_rasters)*2)+1
 feedback = QgsProcessingMultiStepFeedback(steps, model_feedback)
 step = 1
 
 for i, r_lyr in enumerate(input_rasters):
 r_calc_params = {'LAYERS':[r_lyr.source()],
 'EXPRESSION':f'"{r_lyr.name()}@1">100',
 'EXTENT':None,
 'CELL_SIZE':0.008,
 'CRS':'EPSG:4326',
 'OUTPUT':'TEMPORARY_OUTPUT'}
 feedback.setCurrentStep(step)
 step+=1
 outputs[f'r_calc_{i}']=processing.run("native:rastercalc", r_calc_params, feedback=feedback, context=context, is_child_algorithm=True)
 results[f'r_calc_{i}']=outputs[f'r_calc_{i}']['OUTPUT']
 
 vectorize_params = {'INPUT':results[f'r_calc_{i}'],
 'BAND':1,
 'FIELD':'DN',
 'EIGHT_CONNECTEDNESS':False,
 'EXTRA':'',
 'OUTPUT':'TEMPORARY_OUTPUT'}
 feedback.setCurrentStep(step)
 step+=1
 outputs[f'vectorized_{i}']=processing.run("gdal:polygonize", vectorize_params, feedback=feedback, context=context, is_child_algorithm=True)
 results[f'vectorized_{i}']=outputs[f'vectorized_{i}']['OUTPUT']
 
 # Build a list of the vectorized outputs in the results dictionary
 v_lyrs = [v for k, v in results.items() if 'vectorized' in k]
 
 merge_params = {'LAYERS':v_lyrs,
 'CRS':None,
 'OUTPUT':parameters[self.OUTPUT]}
 feedback.setCurrentStep(step)
 step+=1
 outputs['merged']=processing.run("native:mergevectorlayers", merge_params, feedback=feedback, context=context, is_child_algorithm=True)
 results['merged']=outputs['merged']['OUTPUT']
 
 return results

P.s. as a kind of bonus answer to your original question- generally, if you have a layer id string which you need to convert to a map layer object you can use QgsProcessingContext.getMapLayer(). It calls QgsProcessingUtils.mapLayerFromString() which you could also use.

• geotyr: thank you so much for your detailed answer, I could not have dreamed of a better one...although it is still not working you gave me material to highly improve the script, especially the merge part. There are plenty of details I still don not understand, but it is only my second script so I hope I will get there at some point. Regarding your PS note, I understand the instance works for strings with the name of the internal memory layer, and not the full path of a layer and weirdly, when applying the polygonize step, the output is a path to a temporary geopackage file..