[GIS] How to calculate all paired shortest path using PyQGIS

pyqgisqgis

The code for shortest path analysis with PyQGIS specifies start and end nodes explicitly:

pStart = QgsPoint( -1.37144, 0.543836 )

pStop = QgsPoint( -1.1027, 0.699986 )

http://www.qgis.org/pyqgis-cookbook/network-analysis.html

Does anyone know whether it's possible to replace these statements with sets of all nodes, so that it's possible to get all-paired shortest path?

Thanks

Best Answer

OK, so this is what I came up with; it might be useful for novices like me out there. I know it's not perfect at all, so any edits are welcomed!

Given: I've a road network and a point layer which has x, y and TAZ( i.e. Traffic Analysis Zone id) with indexes of 191,192 and 1, respectively.

Desired: a distance matrix for all the points. for example, I'm interested in finding the distance from points with TAZ's of 125 to 135 to all the other points.

    x=[]
    y=[]
    taz=[]
    dict={}
    pntLayer=qgis.utils.iface.mapCanvas().currentLayer() //select point layer in TOC
    feat=QgsFeature()
    provider=pntLayer.dataProvider()
    provider.select([1,191,192],QgsRectangle(), True, False)
    while provider.nextFeature(feat): // fetch attributes, returns qgis objects
    ...     m=feat.attributeMap()
    ...     x.append(m[191])
    ...     y.append(m[192])
    ...     taz.append(m[1])
    ... 

    // I want a dictionary with points as keys and TAZ as values

    xp=[]
    yp=[]
    tazp=[]
    for i in taz:
    ...     tazp.append( str(i.toString()))

    for i in x:
    ...     a,v=i.toDouble()
    ...     xp.append(a)
    ... 
    for i in y:
    ...     a,v=i.toDouble()
    ...     yp.append(a)
    ... 
    qPoints=[]
    for i in xrange(900): // point layer has 900 entries
    ...     qPnt = QgsGeometry.fromPoint(QgsPoint(xp[i],yp[i]))
    ...     qPoints.append(qPnt)   // now I have reconstructed the points from x and y
    ... 
    net=qgis.utils.iface.activeLayer() // select network layer in TOC
    properter=QgsDistanceArcProperter()
    director = QgsLineVectorLayerDirector( net, -1, '', '', '', 3 )
    director.addProperter(properter)
    crs=QgsCoordinateReferenceSystem(26986, QgsCoordinateReferenceSystem.PostgisCrsId)
    builder=QgsGraphBuilder(crs)
    qp=[]
    for i in qPoints:
    ...     x=i.asPoint()
    ...     qp.append(x)
    ... 
    d={}
    for i in range(900):
    ...     d[qp[i]] = tazp[i] // d: dictionary with points as keys with corresponding TAZ
    tiedPoints=director.makeGraph(builder, d.keys()) 
    graph=builder.graph()

    //tiedPoints are now different points from d.keys() but the order is preserved.
    // so I'll create a new dictionary and connect new points to their corresponding TAZ values.
    dict={}
    i=0
    for j in d.keys():
    ...     dict[tiedPoints[i]] = d[j] 
    ...     i=i+1
    ... 
    t={}  // dictionary of graph id and TAZ
    for j in xrange(900):
        id=graph.findVertex(dict.keys()[j])
    ...     t[id] = dict[ dict.keys()[j]]


    costMat={} //cost matrix with graph id as keys 
    for id in t.keys():
    ...     (tree, costs) = QgsGraphAnalyzer.dijkstra(graph, id,0 )
    ...     c=[]
    ...     for i in t.keys():
    ...         c.append(costs[i])
    ...     costMat[id]=c   
    ...
Related Question