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

## 1 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:

• 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 Oct 8, 2017 at 18:25
• I added your suggested requests. Oct 9, 2017 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 Oct 9, 2017 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! Oct 12, 2017 at 22:03