Hálózatelemzés

Válasszuk ki Magyarország gyorforgalmú útjait a Geofabrik által kiadott roads rétegből a következő attribútum szelekcióval:

"type" = 'primary_link' OR "type" = 'primary' OR "type" = 'motorway' OR "type" = 'motorway_link' OR "type" = 'trunk' OR "type" = 'trunk_link'

Majd mentsük csak a szelektált elemeket egyenesen EOV vetületbe. Tesztelés gyanánt vágjunk ki egy kisebb területet az egész országból. Például válasszuk Magyarország legnagyobb megyéjét:

1. ábra: Gyorsforgalmú utak Bács-Kiskun megyében.
1. ábra: Gyorsforgalmú utak Bács-Kiskun megyében.

Továbbá a megye városait és református templomait is felhasználjuk. Ezeket is vágjuk a megyére:

2. ábra: Városok, templomok és utak, Bács-Kiskun megye.
2. ábra: Városok, templomok és utak, Bács-Kiskun megye.

Töltsük be a rétegeket, valamint hozzunk létre egy új mezőt a városoknál:

from qgis.networkanalysis import *    
from PyQt4.QtCore import QVariant    
path = "d:/geodata/network/"    
filename1 = "road_net_clip"    
layer1 = QgsVectorLayer(path + filename1 + ".shp", filename1, "ogr")    
if not layer1.isValid():    
print "Layer1 failed to load!"    
filename2 = "town_city_clip"    
layer2 = QgsVectorLayer(path + filename2 + ".shp", filename2, "ogr")    
if not layer2.isValid():    
print "Layer2 failed to load!"    
if layer2.fieldNameIndex('tdist') == -1:    
layer2.dataProvider().addAttributes([QgsField("tdist", QVariant.Double)])    
layer2.updateFields()    
filename3 = "reformed_clip"    
layer3 = QgsVectorLayer(path + filename3 + ".shp", filename3, "ogr")    
if not layer3.isValid():    
print "Layer3 failed to load!"

Teszteljük a rétegek betöltését:

QgsMapLayerRegistry.instance().addMapLayer(layer1) QgsMapLayerRegistry.instance().addMapLayer(layer2) QgsMapLayerRegistry.instance().addMapLayer(layer3)`

Hozzunk létre egy irányítottságot (director), ami a későbbi gráf struktúra alapjaként fog szolgálni (pl. egyirányú vagy kétirányúak a gráfok), az elméleti sémája a következő:

director = QgsLineVectorLayerDirector(vl, directionFieldId, directDirectionValue, reverseDirectionValue, bothDirectionValue, defaultDirection)

  • directionFieldId — azon oszlop indexe, amely az irányultságot tartalmazza, ha -1, akkor nem használjuk. Egész szám
  • directDirectionValue — mezőérték normális irányultság esetén. Karakterlánc
  • reverseDirectionValue — mezőérték fordított. Karakterlánc bothDirectionValue — mezőrték kétirányú utakhoz. Karakterlánc
  • defaultDirection — alapértelmezett út irányultság. Ezt az értéket hazsnáljuk, amikor nincs előzetes irányultság információnk. Egész szám. 1 – direkt irány, 2 – fordított irány, 3 – mindkét irány.

Használjuk mindkét irányt:

director = QgsLineVectorLayerDirector(layer1, -1, '', '', '', 3) properter = QgsDistanceArcProperter() director.addProperter(properter)

Gyűjtsük fel a pontokat és a szükséges mezők indexeit:

points1 = [feature for feature in layer2.getFeatures()] points2 = [feature for feature in layer3.getFeatures()] idx1 = layer2.fieldNameIndex('name') idx2 = layer2.fieldNameIndex('tdist')

Kezdjünk a ciklus megírásába:

`for point1 in points1:
print point1.attributes()[idx1]
dlist = []
geom1 = point1.geometry()

print geom1.asPoint()

for point2 in points2:
tiedPoints = []
graph = []
builder = []
tree = []
builder = QgsGraphBuilder(layer1.crs(),False, 5.0,"GRS67")
geom2 = point2.geometry()
print geom2.asPoint()
tiedPoints = director.makeGraph(builder,
[QgsPoint(geom1.asPoint()),QgsPoint(geom2.asPoint())])
print tiedPoints
graph = builder.graph()
tStart = tiedPoints[0]
tStop = tiedPoints[1]
idStart = graph.findVertex(tStart)
idStop = graph.findVertex(tStop)
(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
if tree[idStop] != -1:
dlist.append(cost[idStop]/1000)
if not dlist:
dlist.append(-1)
print min(dlist)
attrs = { idx2 : min(dlist)}
layer2.dataProvider().changeAttributeValues({point1.id(): attrs})`

© Dr. Wirth Ervin