1

I need an help to find in QGIS the MAX SLOPE of line from points and peaks along 8 cardinal direction. MAXSLOPE 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!"
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.

  • I don't need to find the TERRAIN MAXSLOPE around the rain gauge but the MAXSLOPE of the imaginary line from rain gauge to any peak in a buffer of 15km (line red in example picture). So your suggestion it's not good for my case, because the SLOPE raster is relative to the ground (grid to grid and not peak to rain gauge point). – Mr Prince Nov 16 '17 at 6:58

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.