Geospatial-Programming-Project/projectScript.py at main · s3763447/Geospatial-Programming-Project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#lOADING LAYERS AND REMOVING PREVIOUS ONES IF EXIST
QgsProject.instance().removeAllMapLayers()
dem=QgsRasterLayer('D:\OneDrive\MASTER\semester4\GeospatialProgramming\Module5\Module5\October14\data\data\dem.tif', "dem")
wdn=QgsVectorLayer('D:\OneDrive\MASTER\semester4\GeospatialProgramming\Module5\Module5\October14\data\data\wdn.shp', "WDN")

vald=''
#testing if layers are valid
if not dem.isValid():
        vald='DEM'
if not wdn.isValid():
        vald=vald+' WDN'

if len(vald)!=0:
    print('Layers are not valid')
else:
    #Loading layers
    QgsProject.instance().addMapLayer(dem)
    QgsProject.instance().addMapLayer(wdn)

    #1...............Create Slope Layer
    #generating the slope layer
    processing.runAndLoadResults("native:slope", {'INPUT': dem, 'Z_FACTOR':1, 'OUTPUT':'TEMPORARY_OUTPUT'})['OUTPUT']
    print('#1............... Slope Layer Generated succefully')
    #2...............Build Newtwork and junkions


    #water Network layer opening
    line_layer=wdn
    #MultiLine String to  LineString
    wdn=processing.run("native:multiparttosingleparts", {'INPUT': line_layer, 'OUTPUT': 'memory:water distribution network' })["OUTPUT"]
    #new layer of juctions
    point_layer = QgsVectorLayer("Point?crs=epsg:32760", "point_layer", 'memory')
    #Remove duplicated points
    junctions=processing.run("native:deleteduplicategeometries", {'INPUT': point_layer, 'OUTPUT': 'memory:junctions' })["OUTPUT"]

    pr = junctions.dataProvider()
    pr1 = wdn.dataProvider()
    junctions.startEditing()
    wdn.startEditing()

    pr.addAttributes([ QgsField("id", QVariant.Int), QgsField("x", QVariant.Double), QgsField("y", QVariant.Double)])
    pr1.addAttributes([ QgsField("start_id", QVariant.Int), QgsField("end_id", QVariant.Int)])
    i=0

    wdn.commitChanges()
    wdn.startEditing()

    junctions.commitChanges()
    junctions.startEditing()

    #loop over Network polylines and creating junctions

    for feature in wdn.getFeatures():
        #create a new feature of type type junction
        feat = QgsFeature(junctions.fields())
        geom = feature.geometry().asPolyline()
        #create new  start end points :junctions
        start_point = QgsPoint(geom[0])
        end_point = QgsPoint(geom[-1])

        feat.setGeometry(start_point)
        feat["id"]=i
        #retreiving the x,y coordinates of polyline start end points and a affect it to x,y columns of junction
        feat["x"]=start_point.x()
        feat["y"]=start_point.y()
        #adding end start feature
        pr.addFeatures([feat])
        #update the start_id of polyline with the id of the start junction
        feature["start_id"]=i
        wdn.updateFeature(feature)
        i+=1

        feat.setGeometry(end_point)
        feat["id"]=i
        feat["x"]=end_point.x()
        feat["y"]=end_point.y()
        #adding end point feature
        pr.addFeatures([feat])
        feature["end_id"]=i
        wdn.updateFeature(feature)
        i+=1

    junctions.commitChanges()
    wdn.commitChanges()
    #add junctions and network layers
    QgsProject.instance().addMapLayer(junctions)
    QgsProject.instance().addMapLayer(wdn)


    print('#2............... Network and Junctions Created  succefully')
    #3............... Fill Junctions with corresponding Altitudes And Slopes
    slope_layer = QgsProject.instance().mapLayersByName('Slope')[0]
    dem_layer = QgsProject.instance().mapLayersByName('dem')[0]
    junctions=QgsProject.instance().mapLayersByName('junctions')[0]
    #retreive altitude value from DEM file and affect it to junction points
    junctions_altitude=processing.run("native:rastersampling", {'INPUT': junctions,'RASTERCOPY': dem_layer,'COLUMN_PREFIX':'altitude','OUTPUT': 'TEMPORARY_OUTPUT'})["OUTPUT"]
    #retreive slope value from Slope layer and affect it to junction points
    junctions_altitude_slope=processing.run("native:rastersampling", {'INPUT': junctions_altitude,'RASTERCOPY': slope_layer,'COLUMN_PREFIX':'slope','OUTPUT': 'memory:junctions_altitude_slope'})["OUTPUT"]


    QgsProject.instance().addMapLayer(junctions_altitude_slope)
    print('#3............... Junctions are Filled succefully')