8

I am working on a rockfall guide for QGIS and I'm currently stuck on one of the stages.

I currently at a point in the guide where I have two point layers one "rockfall release point"-layer and "braking point"-layer.

We are assuming that any rockfall will loosen from a release point, continue down the slope in the steepest direction, then begin to brake when it reaches a braking point, before stopping shortly thereafter.

Rockfall release and braking points with drainage

I was first planning on calculating drainage area for each braking point using watersheds like in ArcGIS Pro, where you can add a "pour point field", and then find which drainage field each release point is placed in (like in the figure above).

ArcGIS Watershed

The only similar thing I can find in QGIS is r.water.outlet where you have to type in one coordinate, do the calculation and start all over again - this is going to be hard considering I have 2.5m points.

A solution similar to the ArcGIS one would probably be the easiest. The other solution that I can come up with is doing a nearest neighbour analysis between the two layers, but then I also have to make sure that it only joins the point that is at the steepest slope.

0

1 Answer 1

4

Here you have a simple script that iterates over points and apply r.water.outlet to each one. You just need to add this code as a new script.

image1

It will ask for three things:

  • Points layer.
  • Flowdir layer (produced from dem file).
  • Output file.

The resulting rasters will contain the id of the point at the end of the filename.

##Rockfall=group
##RockfallAnalysis=name
##points=vector Point
##flowdir=raster
##outputfolder=folder

from PyQt4.QtCore import QFileInfo, QSettings, QVariant
from qgis.core import *
import qgis.utils
import os, glob, processing, string, time, shutil, os.path, subprocess

#crear capa vectorial a partir de la ruta anterior
layer = QgsVectorLayer(points, "shp", "ogr")
features = layer.getFeatures()

#raster extension for grass rwateroutlet
rasterlayer = QgsRasterLayer(flowdir,"rasterlayer")
rasterext = rasterlayer.extent()
xmin = rasterext.xMinimum()
xmax = rasterext.xMaximum()
ymin = rasterext.yMinimum()
ymax = rasterext.yMaximum()
extension = "%f,%f,%f,%f" %(xmin, xmax, ymin, ymax)

#iterate over points
count=int(1)
for f in features:
    geom = f.geometry()
    #create coordinate of point
    a = geom.asPoint()
    point = str(str(int(a[0]))+','+str(int(a[1])))
    id = str(count)
    #name of the raster file
    outputfile = str(outputfolder+"/"+"rockfall"+"_"+id+".tif")
    print ("Processing id:"+' '+id)
    #process
    processing.runalg("grass7:r.water.outlet",
    flowdir,
    point,
    extension,
    0,
    outputfile)
    count+=1

Here some picks:

image2

image3

image4


EDIT: SCRIPT FOR QGIS 3.8

Process resume: layer iteration --> r.water.outlet --> polygonize --> merge all

Here is the complete script ready to be add to your tools: https://github.com/cesarkero/QGIS_scripts/blob/master/QGIS_3.8/WaterOutletIterator.py

It have been hard as I tried to learn as the same time how to create an script for qgis 3.8 where different tools are used. I've set the keep temporal files by default because I needed to know how exactly the different steps where working, but, remember, if you want to execute it with your millions of points UNCHECK that boolean or you HDD/SDD will get out of the game.

Now it is open to you if you want to dissolve those polygons or join the attributes to the "losnepunkt" or "uphill"points...or whatever you are trying to do.

enter image description here

enter image description here

enter image description here

17
  • 2
    Could you share your layers to try something different ? Aug 27, 2019 at 12:31
  • 1
    Here is the data: bit.ly/2Zm6vVt Flowdir.tif, losnepunkt = release point and bremsepunkt = brake point. Aug 27, 2019 at 12:47
  • 1
    It's easy if you just need the final shp and not each one of the raster files. I hope I can try something in the afternoon. Aug 28, 2019 at 11:48
  • 1
    PD: don't be afraid of the code. Just create a new python script, copy it, save that .py and execute it. Sep 5, 2019 at 16:14
  • 1
    Thx a lot! I’ll check it out tomorrow morning 👌🏽 You have done an awesome job César 😊 Sep 5, 2019 at 18:07

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.