3

In the Upslope Area tool from SAGA it used to be possible to import the point file with multiple points and it would create several catchment/basins areas for those points. Now, it is only possible to MANUALLY enter coordinates for only 1 point. However, for a model we're building, we need to perform this analysis 24 times (as we have 24 different points) so we need to automate this process. I tried to use Run as Batch process + calculate by expression and get the XY coordinates from another layer, however, I run into the problem that I need to reference the XY coordinates that are stored in the other layer and looks like there isn't an easy way of doing so. Any suggestions on how to do that? Tried the geomwithin(targetLayer,targetField) from RefFunctions plugin but that produces empty result

so to resume, I have two questions:

  1. How to calculate upslope areas for multiple points?
  2. Is there a way to automatize it in the graphical modeler?
1

1 Answer 1

1

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})

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.