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