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:
