Rockfall analysis (Nearest neighbor downhill) in QGIS

Question

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.

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

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.

Answer 1

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.

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:

---

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.