Calculate Max Slope of line from points and peaks

Question

I need an help to find in QGIS the MAX SLOPE of line from points and peaks along 8 cardinal direction. The procedure is the same of this Script by @xunilk but in this case I need to find the MAX SLOPE and not the MAX ELEVATION DIFFERENCE.

It's really important, in MAX_SLOPE calculus, to ignore ELEVATION DIFFERENCE < 200meters.

I Think a solution could be:

• FOR peak with dif_elev from rain gauge >200meters calculate distance and dif_elev from rain gauge
• search SLOPE MAX with someting like: np.max(dif_elev/distance)

EDIT1: The script should analyze ALL the imaginary line from rain gauge to peaks along the cardinal direction and retain only the steeper one in 15km bufferlenght.

But how can do this with python?

Software: QGIS (python console) Input: Rain gauge POINT Shapefile, DTM 100m grid, Output required: LINEAR Shapefile (with 8 lines for each Rain Gauge and Attribute table with the MAX_SLOPE on each cardinal directions and its cell distance)

MAX_DIF_ELEV SCRIPT:

import numpy as np

print "Processing..."

bufferLength = 15000
polygonSides = 8

registry = QgsMapLayerRegistry.instance()

layer = registry.mapLayersByName('point')
raster = registry.mapLayersByName('utah_demUTM2')

xsize = raster[0].rasterUnitsPerPixelX()

prov_raster = raster[0].dataProvider()

feat_points = [ feat for feat in layer[0].getFeatures() ] 
points = [ feat.geometry().asPoint() for feat in layer[0].getFeatures() ]

epsg = layer[0].crs().postgisSrid()

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'points',
                           'memory')

prov = mem_layer.dataProvider()

group_points = [] 

for i, point in enumerate(points):

    int_points = [ QgsPoint(point[0] + np.sin(angle)*bufferLength, point[1] + np.cos(angle)*bufferLength)
                        for angle in np.linspace(0, 2*np.pi, polygonSides, endpoint = False) ]

    group_points.append(int_points)

lines = []
idx_lines = []

for i, group in enumerate(group_points):
    for point in group:
        lines.append([points[i], point])
        idx_lines.append(i)

RainGaugeCode = {0:1001, 1:1002}

for group in group_points:

    feats = [ QgsFeature() for i in range(len(group)) ]

    for i, feat in enumerate(feats):
        feat.setAttributes([i])
        feat.setGeometry(QgsGeometry.fromPoint(group[i]))

    prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

uri = "LineString?crs=epsg:" + str(epsg) + \
      "&field=id:integer&field=dif_elev:double&field=distance:double&field=tan(ratio):double&field=rain_gauge:integer&field=direction:string""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'Output_line',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(lines)) ]

dif_elev = []
distances = []

directions = ['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW']

for i, feat in enumerate(feats):
    geom = QgsGeometry.fromPolyline(lines[i])
    int_points = []

    for distance in range(0, int(geom.length()), int(xsize)):
        values = []
        point = geom.interpolate(distance)
        pt = point.asPoint()
        int_points.append(pt)
        for p in int_points:
            value = prov_raster.identify(p,
                                         QgsRaster.IdentifyFormatValue).results()[1]
            values.append(value)

    init_value = prov_raster.identify(points[idx_lines[i]],
                                      QgsRaster.IdentifyFormatValue).results()[1]

    elev = np.max(values) - init_value
    distance = values.index(np.max(values))*xsize
    if distance != 0:
        coc_el_dist = elev/distance
    else:
        coc_el_dist = -999

    if i % 8 == 0:
        k = 0

    dif_elev.append(elev)
    distances.append(distance)
    feat.setAttributes([i,float(elev), float(distance), float(coc_el_dist), RainGaugeCode[idx_lines[i]], directions[k]])
    feat.setGeometry(geom)
    k += 1

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

print "Done!"

Answer 1

Use the DEM you have to run a slope aspect in QGIS, and then use that as a second raster input in your script. Where it is already looping through the DEM raster, have it also grab the values of the new slope raster:

for distance in range(0, int(geom.length()), int(xsize)):
    values = []
    point = geom.interpolate(distance)
    pt = point.asPoint()
    int_points.append(pt)
    for p in int_points:
        value = prov_raster.identify(p, QgsRaster.IdentifyFormatValue).results()[1]
        values.append(value)

        init_value = prov_raster.identify(points[idx_lines[i]], QgsRaster.IdentifyFormatValue).results()[1]

        slopevalue = slope_raster.identify(p, QgsRaster.IdentifyFormatValue).results()[1]
        slopevalues.append(slopevalue)

Then just grab the max slope like it is already doing with the elevation:

maxslope = np.max(slopevalues)

You already have the elevation difference, so use it to check if it is less than 200 meters, in which case ignore that slope value or whatever you need to do with it.