1

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)

Example

3

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:

enter image description here

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:

enter image description here

  • Thanks @xunilk, the script works but is it possible to add to "Output_line attribute table" the following columns: <1> the Rain Gauges references for every set of 8 lines (from Rain Gauge attribute table "Field_1") <2> A column with nearest cardinal direction to each line <3> a column with the difference "dif_elev/distance"? EXAMPLE IMG: i.imgur.com/MQA55Z4.png – Mr Prince Oct 8 '17 at 18:25
  • I added your suggested requests. – xunilk Oct 9 '17 at 2:02
  • thanks for your precious time and knowledge!! Solutions for point <2> & <3> works!! But point <1> about Rain Gauge references give me Keyerror:2 i.imgur.com/as2ikMK.png when I have more than 2 Rain Gauges; and the value should be taken from original Rain Gauge.shp Attribute table (same value of Field_1) EXAMPLE IMG: i.imgur.com/AB13wwv.png – Mr Prince Oct 9 '17 at 8:14
  • Is there some suggestion about how could I relate ID FIELD from source Attribute table to the new generated by the script? (i.imgur.com/AB13wwv.png) Please I need only this thing to can really use the script for my project! – Mr Prince Oct 12 '17 at 22:03

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.