1

I am running r.lake in the graphical modeler in QGIS as a component of a geoalgorithm I am working on.

At the moment the coordinates input for r.lake will only take a manual string input, I am wondering is there a way to create a python script that will take the coordinate values of a vector seed layer and input these as the coordinates.

I understand that I will have to iterate / loop over the vector layer if there are multiple coordinates involved.

My python knowledge is that of a beginner, so I am more comfortable with the graphical modeler interface for such tasks. If a script can be linked in the graphical modeler that would be great.


Script is now working well with the addition of several processing algorithms included to vectorize and run spatial query. However, where the spatial query should not return results, blank shapefiles are being produced.

See sample code below:

##Karst Surface Flooding=name
##DEM_input=raster
##Seed_vector_points=vector point
##Overspill_vector_boundary=vector polygon
##output=folder

# Import classes
from qgis.core import QgsRasterLayer
from qgis.utils import iface
from PyQt4.QtCore import QFileInfo
import processing

# Define layers
point_layer = processing.getObject(Seed_vector_points)
layer = processing.getObject(DEM_input)
boundary_layer = processing.getObject(Overspill_vector_boundary)
fileName = layer.source()
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(fileName, baseName)

# Get extent of the region
extent = iface.mapCanvas().extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

for feat in point_layer.getFeatures():
    geom = feat.geometry()
    attr = feat.attribute("Z_ADD_1")
    i = 0
    while i < 5: 
        output_1 = processing.runalg('grass:r.lake.coords', layer, "%f"%(attr),
"%f,%f"% (geom.asPoint().x(), geom.asPoint().y()) , False, "%f,%f,%f,%f"% (xmin, xmax, ymin, ymax), 0.00, None)
        output_2 = processing.runalg('saga:rastercalculator', output_1['lake'], '','a>0', False, 3, None)
        output_3 = processing.runalg('gdalogr:polygonize', output_2['RESULT'], "DN", None)
        output_4 = processing.runalg('qgis:dissolve',output_3['OUTPUT'], True, "DN", None)
        output_5 = processing.runalg('modeler:overspill_select', output_4['OUTPUT'], boundary_layer, output + '/' + baseName + "_%.6f_%.6f" %(geom.asPoint().x(), geom.asPoint().y()) + '_' + "%.2f"%(attr) + 'm' + '.shp')
        i += 1
        attr += 0.5
  • Please provide code as formatted text rather than pictures. – PolyGeo Aug 29 '17 at 4:21
2

You could use something like the following which defines the input layers, iterates through each point in your point layer and runs the r.lake tool. You could create a script from:

Processing Toolbox > Scripts > Tools > Create new script

The output rasters are saved in a directory you specify and is named by using its original name and the coordinates:

##Example=name
##input=raster
##points=vector point
##output=folder

# Import classes
from qgis.core import QgsRasterLayer
from qgis.utils import iface
from PyQt4.QtCore import QFileInfo

# Define layers
point_layer = processing.getObject(points)
layer = processing.getObject(input)
fileName = layer.source()
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(fileName, baseName)

# Get extent of the region
extent = rlayer.extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

for feat in point_layer.getFeatures():
    geom = feat.geometry()
    processing.runalg('grass7:r.lake.coords', layer, 1000, "%f,%f"% (geom.asPoint().x(), geom.asPoint().y()), False, "%f,%f,%f,%f"% (xmin, xmax, ymin, ymax), 0.00, output + '/' + baseName + "_%.6f_%.6f" %(geom.asPoint().x(), geom.asPoint().y()) + '.tif')

Save the script in your .qgis2/processing/scripts/ directory. Then add this script to your modeler which requires 3 input parameters as shown in a simple model:

Modeler

When you run the model, you can specify your raster layer, points layer and the output directory.


Tested on QGIS 2.18.2 for Win7 64-bit.

  • @Ian - Most welcome, hope it helped :) – Joseph Mar 7 '17 at 12:48
  • I have only resumed work on this in the past week or so due to other commitments. I have it working pretty well now. I have also added in processing algorithms in the script to vectorize the data. The last step I have in the script is an extract by location algorithm which I had hoped would allow me to produce only certain vectors that were within a larger polygon. The script does this successfully, however, it still outputs blank shapefiles where the spatial query didn't extract results. I would greatly appreciate your help. – Ian Aug 28 '17 at 9:08
  • @Ian - Can you include a screenshot of your modeler please? Did you set the Extract by location tool to have an output? – Joseph Aug 28 '17 at 9:39
  • I have included a screenshot in my edited question above. – Ian Aug 28 '17 at 10:09
  • @Ian - Have you tried calling the Extract by location tool directly from your script instead of calling your model which I assume uses this tool? Or does your model include other tools? – Joseph Aug 29 '17 at 9:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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