/
mca_tools.py
executable file
·461 lines (380 loc) · 15.1 KB
/
mca_tools.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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
# Build-in dependencies
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.networkanalysis import *
from qgis.utils import *
import math
import sys
import inspect
# Custom cost builder
from arc_properter import customProperter
# Loading shapely and SciPy
try:
<<<<<<< Updated upstream
from shapely.ops import cascaded_union, polygonize
from shapely.geometry import MultiPoint
from scipy.spatial import Delaunay
ex_dep_loaded = True
=======
from shapely.ops import cascaded_union, polygonize
from shapely.geometry import MultiPoint, MultiLineString
from scipy.spatial import Delaunay
ex_dep_loaded = True
>>>>>>> Stashed changes
except ImportError,e:
ex_dep_loaded = False
def graph_builder(network_lines, origin_points, origins_column, tolerance,custom_cost, cost_column):
# Settings
crs = network_lines.crs()
epsg = crs.authid()
otf = False
default_value = 0
network_fields = network_lines.pendingFields()
cost_index = network_fields.indexFromName(cost_column)
# Setting up graph build director
director = QgsLineVectorLayerDirector(network_lines, -1, '', '', '', 3)
# Determining cost calculation
if custom_cost == True:
properter = customProperter(cost_index,default_value)
else:
properter = QgsDistanceArcProperter()
# Building graph
director.addProperter(properter)
builder = QgsGraphBuilder(crs, otf, tolerance, epsg)
# Reading origins and making list of coordinates
origins = []
origins_name = {}
for i,f in enumerate(origin_points.getFeatures()):
geom = f.geometry().asPoint()
origins.append(geom)
if origins_column:
origins_name[i] = f[origins_column]
# Connect origin points to the director and build graph
tied_origins = director.makeGraph(builder, origins)
graph = builder.graph()
print origins_name
return graph, tied_origins, origins_name
def alpha_shape(points, alpha):
<<<<<<< Updated upstream
# Empty triangle list
pl_lines = []
# Transform points into Shapely's MultiPoint format
multi_points = MultiPoint(points)
# Create delaunay triangulation
triangles = Delaunay(multi_points)
# Assess triangles
for a, b, c in triangles.vertices:
coord_a = points[a]
coord_b = points[b]
coord_c = points[c]
# Calculating length of triangle sides
a = math.sqrt((coord_a[0] - coord_b[0]) ** 2 + (coord_a[1] - coord_b[1]) ** 2)
b = math.sqrt((coord_a[0] - coord_c[0]) ** 2 + (coord_a[1] - coord_c[1]) ** 2)
c = math.sqrt((coord_c[0] - coord_b[0]) ** 2 + (coord_c[1] - coord_b[1]) ** 2)
# Test triangle area
if ((a + b + c) * (b + c - a) * (c + a - b) * (a + b - c)) <=0:
pass
else :
# Calculating circumcircle radius
circum_rad = (a * b * c) / math.sqrt((a + b + c) * (b + c - a) * (c + a - b) * (a + b - c))
# Circumcircle radius filter
if circum_rad < alpha:
pl_lines.append((coord_a, coord_b))
pl_lines.append((coord_a, coord_c))
pl_lines.append((coord_c, coord_b))
# Circumcircle radius filter
pl_triangles = list(polygonize(pl_lines))
pl_polygon = cascaded_union(pl_triangles)
return pl_polygon
=======
# Empty triangle list
pl_lines = []
# Transform points into Shapely's MultiPoint format
multi_points = MultiPoint(points)
# Create delaunay triangulation
triangles = Delaunay(multi_points)
# Assess triangles
for a, b, c in triangles.vertices:
coord_a = points[a]
coord_b = points[b]
coord_c = points[c]
# Calculating length of triangle sides
a = math.sqrt((coord_a[0] - coord_b[0]) ** 2.0 + (coord_a[1] - coord_b[1]) ** 2.0)
b = math.sqrt((coord_a[0] - coord_c[0]) ** 2.0 + (coord_a[1] - coord_c[1]) ** 2.0)
c = math.sqrt((coord_c[0] - coord_b[0]) ** 2.0 + (coord_c[1] - coord_b[1]) ** 2.0)
# Semi-perimeter of the triangle
s = (a + b + c)/2.0
# Area of
tri_area = math.sqrt(s * (s-a) * (s-b)) * (s-c)
# Calculating circumcircle radius and area
circum_rad = a * b * c / (4.0 * tri_area)
circum_area = 3.14159 * circum_rad
#(a * b * c) / math.sqrt((a + b + c) * (b + c - a) * (c + a - b) * (a + b - c))
#print circum_rad
# Circumcircle radius filter
if circum_area < tri_area * 5:
pl_lines.append((coord_a, coord_b))
pl_lines.append((coord_a, coord_c))
pl_lines.append((coord_c, coord_b))
# Writing the polygon
pl_triangles = list(polygonize(pl_lines))
pl_polygon = cascaded_union(pl_triangles)
return pl_polygon
>>>>>>> Stashed changes
def mca_network_writer(output_network, mca_network):
output_network.dataProvider().addAttributes([QgsField("min_dist", QVariant.Int)])
output_network.updateFields()
arc_geom_length = []
i = 0
for arc in mca_network:
arc_geom = QgsGeometry.fromPolyline(arc['arcGeom'])
if arc_geom.length() in arc_geom_length: # removes duplicates
pass
elif not arc['arcCost']: # removes unconnected lines or lines out of reach
pass
else:
f = QgsFeature(output_network.pendingFields())
f.setAttribute("id", i)
geom = QgsGeometry.fromPolyline(arc['arcGeom'])
f.setGeometry(geom)
cost_list = []
for cost_dict in arc['arcCost']:
for origin_index,origin in enumerate(cost_dict):
cost_list.append(cost_dict[origin])
f.setAttribute(origin_index + 1, int(cost_dict[origin]))
if len(cost_list) > 0:
f.setAttribute('min_dist', int(min(cost_list)))
output_network.dataProvider().addFeatures([f])
arc_geom_length.append(arc_geom.length())
i += 1
def mca_catchment_writer(output_catchment, mca_catchments, origins_name, alpha):
for i,j in enumerate(mca_catchments):
origin = j.keys()[0]
points = j[origin]
p = QgsFeature(output_catchment.pendingFields())
if origins_name:
p.setAttribute("origin", "%s" % origins_name.get(i))
else:
p.setAttribute("origin", "origin_%s" % i)
p_geom = QgsGeometry.fromWkt((alpha_shape(points, alpha)).wkt)
p.setGeometry(p_geom)
output_catchment.dataProvider().addFeatures([p])
def mca_vector_writer(layer, path, crs):
shp_writer = QgsVectorFileWriter.writeAsVectorFormat(
layer,
r"%s" % path,
"utf-8",
crs,
"ESRI Shapefile")
def mca(graph,
<<<<<<< Updated upstream
tied_origins,
origins_name,
output_network,
output_catchment,
alpha,
radius):
output_network.dataProvider().addAttributes([QgsField("id", QVariant.Int)])
output_network.updateFields()
# dictionary with id's, list of lines and their respective costs
mca_network = []
# list with dictionaries of polygon points
mca_catchments = []
n = 0
# distance calculation for all the arcs
while n < graph.arcCount():
# setup the network properties dictionary
arc_prop = {'arcId': 0, 'arcGeom': [], 'arcCost': [],}
# give id
arc_prop['arcId'] = n;
# get the vertex id's
inVertexId = graph.arc(n).inVertex()
outVertexId = graph.arc(n).outVertex()
# find the vertex geometries
inVertexGeom = graph.vertex(inVertexId).point()
outVertexGeom = graph.vertex(outVertexId).point()
# add the geometries to the network properties dictionary
arc_prop['arcGeom'] = [inVertexGeom, outVertexGeom];
mca_network.append(arc_prop)
# next arc
n += 1
# iteration through all tied origin points
i = 0
for o in tied_origins:
mca_catchment_points = {i: []}
origin_vertex_id = graph.findVertex(tied_origins[i])
# update output network fields
if origins_name:
origin_field_name = str(origins_name.get(i))
else:
origin_field_name = "origin_%s" % (i + 1)
output_network.dataProvider().addAttributes([QgsField("%s" % (origin_field_name), QVariant.Int)])
output_network.updateFields()
# run the graph builder for origin
(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, origin_vertex_id, 0)
# Analysing the costs of the tree
x = 0
while x < graph.arcCount():
inVertexId = graph.arc(x).inVertex()
# origin lines
if inVertexId == origin_vertex_id:
arcCost = {i: 0}
mca_network[x]['arcCost'].append(arcCost)
# lines within the radius
if cost[inVertexId] < radius and tree[inVertexId] != -1:
outVertexId = graph.arc(x).outVertex()
if cost[outVertexId] < radius:
arc_cost = cost[outVertexId]
arcCost = {i: arc_cost}
mca_network[x]['arcCost'].append(arcCost)
mca_catchment_points[i].append(graph.vertex(inVertexId).point())
# lines at the edge of the radius
elif cost[inVertexId] > radius and tree[inVertexId] != -1:
outVertexId = graph.arc(x).outVertex()
if cost[outVertexId] < radius:
# constructing cut down edge lines
edge_line_length = radius - cost[outVertexId]
edge_line_azimuth = graph.vertex(outVertexId).point().azimuth(
graph.vertex(inVertexId).point()) # degrees from north
new_point_adjacent = math.sin(math.radians(edge_line_azimuth)) * edge_line_length
new_point_opposite = math.cos(math.radians(edge_line_azimuth)) * edge_line_length
new_point_x = graph.vertex(outVertexId).point()[0] + new_point_adjacent
new_point_y = graph.vertex(outVertexId).point()[1] + new_point_opposite
# add edge lines
l = QgsFeature(output_network.pendingFields())
if origins_name:
old_cost = l[origin_field_name]
# check if there is already a cost
if old_cost > int(cost[outVertexId]):
l.setAttribute(origin_field_name, int(cost[outVertexId]))
else:
l.setAttribute(i + 1, int(cost[outVertexId]))
l.setGeometry(QgsGeometry.fromPolyline(
[graph.vertex(outVertexId).point(), QgsPoint(new_point_x, new_point_y)]))
output_network.dataProvider().addFeatures([l])
# add edge points
mca_catchment_points[i].append((new_point_x, new_point_y))
x += 1
mca_catchments.append(mca_catchment_points)
i += 1
# running writers
mca_network_writer(output_network, mca_network)
mca_catchment_writer(output_catchment, mca_catchments, origins_name, alpha)
=======
tied_origins,
origins_name,
output_network,
output_catchment,
alpha,
radius):
output_network.dataProvider().addAttributes([QgsField("id", QVariant.Int)])
output_network.updateFields()
# dictionary with id's, list of lines and their respective costs
mca_network = []
# list with dictionaries of polygon points
mca_catchments = []
n = 0
while n < graph.arcCount():
arc_prop = {'arcId': 0, 'arcGeom': [], 'arcCost': [],}
arc_prop['arcId'] = n;
inVertexId = graph.arc(n).inVertex()
outVertexId = graph.arc(n).outVertex()
inVertexGeom = graph.vertex(inVertexId).point()
outVertexGeom = graph.vertex(outVertexId).point()
arc_prop['arcGeom'] = [inVertexGeom, outVertexGeom];
mca_network.append(arc_prop)
n += 1
# iteration through all tied origin points
i = 0
for o in tied_origins:
mca_catchment_points = {i: []}
origin_vertex_id = graph.findVertex(tied_origins[i])
if origins_name:
origin_field_name = str(origins_name.get(i))
else: # No name in record
origin_field_name = "origin_%s" % (i + 1)
output_network.dataProvider().addAttributes([QgsField("%s" % (origin_field_name), QVariant.Int)])
output_network.updateFields()
(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, origin_vertex_id, 0)
# Analysing the costs of the tree
x = 0
while x < graph.arcCount():
inVertexId = graph.arc(x).inVertex()
# origin lines
if inVertexId == origin_vertex_id:
arcCost = {i: 0}
mca_network[x]['arcCost'].append(arcCost)
# lines within the radius
if cost[inVertexId] < radius and tree[inVertexId] != -1: # lines within the radius
outVertexId = graph.arc(x).outVertex()
if cost[outVertexId] < radius:
arc_cost = cost[outVertexId]
arcCost = {i: arc_cost}
mca_network[x]['arcCost'].append(arcCost)
# mca_catchment_points[i].append(graph.vertex(inVertexId).point())
# lines at the edge of the radius
elif cost[inVertexId] > radius and tree[inVertexId] != -1:
outVertexId = graph.arc(x).outVertex()
if cost[outVertexId] < radius:
# constructing cut down edge lines
edge_line_length = radius - cost[outVertexId]
edge_line_azimuth = graph.vertex(outVertexId).point().azimuth(
graph.vertex(inVertexId).point()) # degrees from north
new_point_adjacent = math.sin(math.radians(edge_line_azimuth)) * edge_line_length
new_point_opposite = math.cos(math.radians(edge_line_azimuth)) * edge_line_length
new_point_x = graph.vertex(outVertexId).point()[0] + new_point_adjacent
new_point_y = graph.vertex(outVertexId).point()[1] + new_point_opposite
# add edge lines
l = QgsFeature(output_network.pendingFields())
l.setAttribute(i + 1, int(cost[outVertexId]))
l.setGeometry(QgsGeometry.fromPolyline(
[graph.vertex(outVertexId).point(), QgsPoint(new_point_x, new_point_y)]))
output_network.dataProvider().addFeatures([l])
# add edge points
mca_catchment_points[i].append((new_point_x, new_point_y))
x += 1
mca_catchments.append(mca_catchment_points)
i += 1
# running writers
mca_network_writer(output_network, mca_network)
mca_catchment_writer(output_catchment, mca_catchments, origins_name, alpha)
>>>>>>> Stashed changes
def mca_network_renderer(output_network, radius):
# settings for 10 color ranges depending on the radius
color_ranges = (
(0, (0.1 * radius), '#ff0000'),
((0.1 * radius), (0.2 * radius), '#ff5100'),
((0.2 * radius), (0.3 * radius), '#ff9900'),
((0.3 * radius), (0.4 * radius), '#ffc800'),
((0.4 * radius), (0.5 * radius), '#ffee00'),
((0.5 * radius), (0.6 * radius), '#a2ff00'),
((0.6 * radius), (0.7 * radius), '#00ff91'),
((0.7 * radius), (0.8 * radius), '#00f3ff'),
((0.8 * radius), (0.9 * radius), '#0099ff'),
((0.9 * radius), (1 * radius), '#0033ff'))
# list with all color ranges
ranges = []
# for each range create a symbol with its respective color
for lower, upper, color in color_ranges:
symbol = QgsSymbolV2.defaultSymbol(output_network.geometryType())
symbol.setColor(QColor(color))
symbol.setWidth(0.5)
range = QgsRendererRangeV2(lower, upper, symbol, '')
ranges.append(range)
# create renderer based on ranges and apply to network
renderer = QgsGraduatedSymbolRendererV2('min_dist', ranges)
output_network.setRendererV2(renderer)
# add network to the canvas
QgsMapLayerRegistry.instance().addMapLayer(output_network)
def mca_catchment_renderer(output_catchment):
# create a black dotted outline symbol layer
symbol_layer = QgsMarkerLineSymbolLayerV2()
symbol_layer.setColor(QColor('black'))
symbol_layer.setWidth(1)
# create renderer and change the symbol layer in its symbol
renderer = output_catchment.rendererV2()
renderer.symbols()[0].changeSymbolLayer(0, symbol_layer)
output_catchment.setRendererV2(renderer)
# add catchment to the canvas
QgsMapLayerRegistry.instance().addMapLayer(output_catchment)