qgis - How to calculate difference in elevation between a POINT and the surrounding cardinal directions?

I have a POINT shapefile (Rain gauges) and a 100m DTM in the same coordinate system.

Is there a method using QGis to calculate the difference in elevation between the elevation of each rain gauge point and the surrounding maximum height along the 8 cardinal directions (maximum distance from Rain Gauge 30km)?

Software: QGIS Input: Rain gauge POINT Shapefile, DTM 100m grid, Output required: LINEAR Shapefile (with 8 lines for each Raing Gauge and Attribute table with the elevation difference between each Raing Gauge and hightest point on each cardinal directions and its cell distance)

Answer

It can be done with following PyQGIS code:

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)

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""&index=yes"

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

prov = mem_layer.dataProvider()

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


dif_elev = []
distances = []

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

    dif_elev.append(elev)
    distances.append(distance)
    feat.setAttributes([i,float(elev), float(distance)])
    feat.setGeometry(geom)

prov.addFeatures(feats)


QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

print "Done!"

It uses a numpy method for generating points corresponding to eight sides polygon around each original points (Rain Gauge) and grouped them by pairs to produce each cardinal directions as line memory layer.

Afterward, each cardinal direction is interpolated by cell distance and each point is used to get raster values along that direction for selecting maximum value and corresponding index.

These values are used to compute elevation differences between each original points (Rain Gauge) and highest point on each cardinal directions and its cell distance and introducing them in attributes table of line vectorial layer.

I tried above code out with layers of following image:

Value Tool QGIS plugin was used to corroborate that code produces expected results. At attributes table of line memory layer, selecting first feature, it can be observed that original point is highest point. Then, elevation difference and its cell distance are both zero.

Editing Note:

New code based in your commentary:

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!"

and result after running it at Python Console of QGIS: