-
Notifications
You must be signed in to change notification settings - Fork 0
/
wrapper.py
140 lines (121 loc) · 5.13 KB
/
wrapper.py
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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
"""
Wrapper script to get path from a point with an algorithm
Author : Ismail Sunni
Email : imajimatika@gmail.com
Date : Jun 2019
"""
from osgeo import ogr, osr
import os
import networkx as nx
from qgis.core import QgsApplication, QgsVectorLayer, QgsPointXY, QgsGeometry, QgsFeature
from utils import get_points, get_spatial_reference
from qgis_utils import get_nearest_feature
from a_star import shortest_path_a_star as a_star_algorithm
from a_star_landmark import shortest_path_a_star as a_star_landmark_algorithm
def algorithm_wrapper(start_point, end_point, node_layer, input_data_path, output_file, algorithm, node_id_attribute='nodeID'):
"""
`start_point` : the starting point as QgsPointXY
`end_point` : the end point as QgsPointXY
`node_layer` : a point vector layer that contains the node. Naturally it's located in input_data_path
`input_data_path` : a path to directory with nodes and edges layer.
`output_file` : a path to the output shape file.
"""
# Get the nearest node from start and end
start_node = get_nearest_feature(node_layer, start_point)
print('start nodeID:', start_node['nodeID'])
end_node = get_nearest_feature(node_layer, end_point)
print('end nodeID:', end_node['nodeID'])
# generate path with the `algorithm`
start_node_id = (node_id_attribute, start_node['nodeID'])
end_node_id = (node_id_attribute, end_node['nodeID'])
path = algorithm(start_node_id, end_node_id, input_data_path, output_file)
# Sort the path
coordinates = [start_point]
G = nx.Graph(nx.read_shp(input_data_path, strict=False, geom_attrs=True)) # Read and convert to Graph
for i in range(len(path) - 1):
node1 = path[i]
node2 = path[i+1]
edge_key = [node1, node2]
# print('edge key: ', edge_key)
points = get_points(G, edge_key)
# print(points)
if i == 0:
for p in points:
coordinates.append(QgsPointXY(p[0], p[1]))
else:
start_point_path = QgsPointXY(points[0][0], points[0][1])
if start_point_path == coordinates[-1]:
# Same end point, add a usual
for p in points[1:]:
coordinates.append(QgsPointXY(p[0], p[1]))
else:
# Different endpoint, reverse
for p in points[::-1][1:]:
coordinates.append(QgsPointXY(p[0], p[1]))
# Append end point
coordinates.append(end_point)
# Spatial reference
spatial_reference = get_spatial_reference(input_data_path)
# Write result to a shapefile (TODO: put it in a function)
# Create geometry for the whole line
line_geometry = QgsGeometry.fromPolylineXY(coordinates)
# Convert to Wkt for ogr
line_wkt = line_geometry.asWkt()
# set up the shapefile driver
driver = ogr.GetDriverByName("ESRI Shapefile")
# create the data source
data_source = driver.CreateDataSource(output_file)
# create the layer
layer = data_source.CreateLayer("A Star Shortest Path", spatial_reference, ogr.wkbLineString)
feature = ogr.Feature(layer.GetLayerDefn())
geom = ogr.CreateGeometryFromWkt(line_wkt)
# Set the feature geometry using the geom
feature.SetGeometry(geom)
# Create the feature in the layer (shapefile)
layer.CreateFeature(feature)
# Dereference the feature
feature = None
data_source = None
print(output_file)
if __name__ == "__main__":
from qgis.core import QgsApplication, QgsVectorLayer, QgsPointXY
QgsApplication.setPrefixPath('/usr', True)
qgs = QgsApplication([], False)
qgs.initQgis()
routes = [
[
'A',
QgsPointXY(405812, 5756851),
QgsPointXY(404984, 5757671)
],
[
'B',
QgsPointXY(406423, 5757091),
QgsPointXY(405600, 5757755)
],
[
'C',
QgsPointXY(405814, 5758009),
QgsPointXY(404769, 5757892)
]
]
algorithms = [
a_star_algorithm,
a_star_landmark_algorithm
]
# Test
# start_point = QgsPointXY(404346.48, 5757913.23)
# end_point = QgsPointXY(404447.56, 5757876.82)
node_path = '/home/ismailsunni/Documents/GeoTech/Routing/topic_data/nodes_single.shp'
# node_path = '/home/ismailsunni/Documents/GeoTech/Routing/processed/small_data/smaller_nodes_single.shp'
node_layer = QgsVectorLayer(node_path, 'nodes', 'ogr')
if not node_layer.isValid():
print('Node layer is not valid.')
input_data_path = '/home/ismailsunni/Documents/GeoTech/Routing/topic_data'
# input_data_path = '/home/ismailsunni/Documents/GeoTech/Routing/processed/small_data'
base_output_file = '/home/ismailsunni/dev/python/routing/test/output/route_landmark_survey'
output_file = os.path.join(base_output_file, 'landmark_wrapper_A.shp')
for route in routes:
for algorithm in algorithms:
output_file = os.path.join(base_output_file, algorithm.name + '_' + route[0] + '.shp')
algorithm_wrapper(route[1], route[2], node_layer, input_data_path, output_file, algorithm)