I wrote a python script that calculates upslope areas for multiple points. Your point shapefile needs to have x, y coordinates stored in the attribute table (you can add it via the field calculator). It is necessary to provide a physical path to the outputs, as 1. SAGA does not support memory layers and 2. in case you rerun the script it will add new features to the already existing memory file which will result in many features being saved in the memory file. This script works with QGIS 3.10.1
import processing
import os # This is is needed in the pyqgis console also
from qgis.core import (QgsVectorLayer)
#provide path to outflow points where x, y are stored in fields named 'x', 'y' (or rename geom.asPoint() to match your coordinates fields)
path_outflow_points = 'your_path/Ouflow_points_new.shp'
#provide path to filled and clipped DEM
path_DEM = 'yourpath/DEM_filled.sdat'
#provide path where to store temporary outputs from processing algorythims
path_to_temp_output1 = 'your_path/upslope_raster{}.sdat'
path_to_temp_output2 = 'your_path/polygonize{}.shp'
path_to_temp_output3 = 'your_path/extract{}.shp'
#provide path to final output
path_final = 'your_path/merged_outflow.shp'
#adds outflow point vector file to processing
outflow_points = QgsVectorLayer(path_outflow_points, "Outflow points layer", "ogr")
#creates an empty list to which extracted polygons will be appended and then this list will be used for merging all the polygons
output_3_list=[]
#iterates over all points in the outflow file. i is going to serve as polygone identifier
for i,f in enumerate(outflow_points.getFeatures()):
#specifies geometry
geom = f.geometry()
#creates a raster file with upslope area
output_1 = processing.run("saga:upslopearea", {'TARGET_PT_X': geom.asPoint().x(), 'TARGET_PT_Y': geom.asPoint().y(), 'ELEVATION': path_DEM, 'METHOD': 0, 'AREA': path_to_temp_output1.format(str(i))})
#converts raster to polygon
output_2 = processing.run("gdal:polygonize", {'INPUT': output_1['AREA'], 'BAND': '1', 'OUTPUT': path_to_temp_output2.format(str(i))})
#extracts only polygons with value 100
output_3 = processing.run("qgis:extractbyattribute", {'INPUT': output_2['OUTPUT'], 'FIELD': 'DN', 'OPERATOR': '0', 'VALUE': '100', 'OUTPUT': path_to_temp_output3.format(str(i))})
#adds a file name to the list that will be used for merging
output_3_list.append(output_3['OUTPUT'])
#merges all extracted polygons into one shapefile
output_4 = processing.runAndLoadResults('qgis:mergevectorlayers', {'LAYERS' : output_3_list, 'OUTPUT': path_final})