[GIS] Create shapefile line from a list of point coordinates in a stand-alone python script

dijkstralinepyqgis

I've used pyqgis to calculate the shortest path through a network between two points. I want to turn this shortest path into a polyline.

I've managed to use the QGIS dijkstra function to calculate a shortest path between two points in a network.
If I print the resulting variable point_route as a string, it prints a list of coordinates.

Where do I go from here?

Is there a(nother) well-documented python library that I can use to turn this list of coordinates into a single polyline shapefile? I can't seem to find/understand the right instructions for pyqgis.

Relevant code
(copied from the shortest path chapter of 'pyqgis developer cookbook'):

(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, startId, 0)

if tree[endId] == -1:
    print "Path not found"
else:
    point_route = []        
    curPos = endId                                                              
    while curPos != startId:                                                    
        point_route.append(graph.vertex(graph.arc(tree[curPos]).inVertex()).point())  
        curPos = graph.arc(tree[curPos]).outVertex()                            
        print "appended " + str(curPos)                                         
    point_route.append(tStart)

Best Answer

The result of the PyQGIS script is a list of points (point_route)

In this example, the point_route result is (7 points between pStart = QgsPoint(110579.17281126765, 116585.26070464966), point 7, and pStop = QgsPoint(111819.28396552203, 116854.75854462344)), point 22)

[(111819,116855), (111576,116688), (111306,116504), (111145,116599), (111005,116681), (110724,116484), (110579,116585)]

If you want to transform these points in a LineString

line = QgsGeometry.fromPolyline( point_route)
print line.exportToWkt()
LineString (111819.28396552200138103 116854.75854462340066675, 111575.63678446662379429 116688.29510413877142128, 111305.62611302707227878 116503.81972814344044309, 111144.6106812154466752 116598.9506020281551173, 111004.94058255387062673 116681.47026111763261724, 110723.93371811544056982 116483.65621838794322684, 110579.17281126769375987 116585.26070464969961904)
print print line.exportToGeoJSON()
{"type": "LineString", "coordinates": [ [111819.28396552200138103, 116854.75854462340066675], [111575.63678446662379429, 116688.29510413877142128], [111305.62611302707227878, 116503.81972814344044309], [111144.6106812154466752, 116598.9506020281551173], [111004.94058255387062673, 116681.47026111763261724], [110723.93371811544056982, 116483.65621838794322684], [110579.17281126769375987, 116585.26070464969961904]]}

And you can use pyQGIS to save the result to a shapefile (in green).

But if you want a pure Python script you don't need QGIS (standalone script), simply use the module NetworkX (many examples in GIS SE) in combination with the module Shapely

import networkx as nx
G  = nx.read_shp("shapefile.shp")
# transform the directed graph to undirect
G = G.to_undirected()
# check that there is a path between the nodes/points 7 and 22
nx.has_path(G,(110579.17281126765, 116585.26070464966),(111819.28396552203, 116854.75854462344))
True 
# Dijkstra algorithm
nx.dijkstra_path(G,(110579.17281126765, 116585.26070464966),(111819.28396552203, 116854.75854462344))
[(111819,116855), (111576,116688), (111306,116504), (111145,116599), (111005,116681), (110724,116484), (110579,116585)]
# convert to polyline with Shapely
 result = LineString(nx.dijkstra_path(G,(110579.17281126765, 116585.26070464966),(111819.28396552203, 116854.75854462344)))
 print result.wkt
 LineString (111819.28396552200138103 116854.75854462340066675, 111575.63678446662379429 116688.29510413877142128, 111305.62611302707227878 116503.81972814344044309, 111144.6106812154466752 116598.9506020281551173, 111004.94058255387062673 116681.47026111763261724, 110723.93371811544056982 116483.65621838794322684, 110579.17281126769375987 116585.26070464969961904)
 # or specifying a weight
 nx.shortest_path(G,(110579.17281126765, 116585.26070464966),(111819.28396552203, 116854.75854462344),weight='distance')
 [(110579.17281126765, 116585.26070464966), (110723.93371811544, 116483.65621838794), (110996.30352863036, 116292.48585842559), (111305.62611302707, 116503.81972814344), (111575.63678446662, 116688.29510413877), (111819.28396552203, 116854.75854462344)]

And you can use NetworkX, Fiona or pyshp (as Clubdebambos propose) to save the result to a shapefile

Related Solutions

QGIS – Importing QGIS Processing in Standalone Python Script

Linux QGIS 2.18.1

With this code a got it, run Processing from standalone script:

#!/usr/bin/env python
import qgis
from qgis.core import *
import sys

app = QgsApplication([],True, None)
app.setPrefixPath("/usr", True)
app.initQgis()
sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

print Processing.getAlgorithm("qgis:creategrid")

[GIS] Shortest path between cities without duplicate calculation (many-to-many)

Since no one joined in, and if you don´t mind, I´ll share some thoughts without this being a satisfying answer (I share a work-related intereset in this).

I highly doubt you have much of an alternative to get a generalized network like the one you are looking for; this is due to the fact that the adjacency matrix of neighbouring cities you imply (to be able to, sort of, avoid unnecessary calculations) is based on a shortest path calculation as well, and for it to be the shortest path connecting each pair, well, by definition you will each have to traverse the whole graph to find it (or use a heuristic algorithm like A* maybe with its pros and cons).

In your example, as Evan said, there is no way to determine that an existing path {1, 2} is part of the shortest path of [1, 3] before walking all (dijkstra) or lots of (A*) other paths. Also, you shouldn't simply assume that if you reach a city with an existing path to the target of the running search, this is the shortest path and skip the search (imagine a triangle 'network': if you already have the path {B, C} and reach B in your search for [A, C], you shouldn`t skip and take {A, B, C} due to your assumption...).

With this said, I think you could at least try and speed up the 'all-to-all' search in some ways.

While the pgr_dijkstra(many-to-many) function is certainly a lot more sophisticated than what I could propose, writing your own function to do an 'all-to-all' search with pairwise one-to-one-dijkstra could give you the benefit of control:

For example, if you subsequentially merge all found paths to single linestrings (with accumulated cost and updated source/target of respective nodes) and recursively route on the updated network, you´ll get your generalized network in one turn without post-processing and the algorithm might get a speed boost with decreasing node/edge count.
(You´d have to clean your base network accordingly (i.e. get rid of dangling nodes/edges if there is no city assigned [you can use pgr_analyzeGraph to find those]) and at least keep track of visited cities to avoid reverse search.)
Additionally, you can try and carefully get only subsets of your data for each search, based on the bounding box of the city pair in question (fast) or by radius search (medium fast) or something similar. That could reduce the individual runtime a lot (or fetching those subsets for each iteration might defeat that purpose...not 100% sure, but I think it does help, with indexes and cleaned tables), but could also lead to unexpected errors (i.e. no path found due to cut-off network).

This should 'easily' be doable in PL/PgSQL (if you are familiar), maybe even by adjusting the pgr_dijkstra(many-to-many) function itself to your needs (although unlikely, due to external C dependencies or general inaccessibility...I haven´t checked).
You could as well switch to Python scripting with one of the many dijkstra implementations (I googled it and came up with heaps...I hope linking them here/now is not necessary), but you´d at least lose the benefits of spatial indexes (the whole control advantage is even more present then, though, since you can easily implement some break conditions in the main algorithm to skip the search whenever you want...if you really want).

As I said, I´m on it, for my own project. If I find something, I´ll update my answer. Maybe someone joins in with an ingenius (or obvious...) solution...

Best Answer

Related Solutions

QGIS – Importing QGIS Processing in Standalone Python Script

[GIS] Shortest path between cities without duplicate calculation (many-to-many)

Related Question