I'm working in QGIS (3.8) and PostGIS (2.5) and I have 2 layers: 1. set of points, 2. set of OSM roads.

enter image description here

I would like to generate a subset of the roads network which represents the minimum spanning tree that connects all the points, something like this:

enter image description here

How do I do this using QGIS-PostGIS?

The roads map can be downloaded from this link as a shapefile: https://download.geofabrik.de/

  • 1
    docs.pgrouting.org/dev/en/pgr_kruskal.html is what you are looking for. Or possibly, docs.pgrouting.org/dev/en/prim-family.html. – John Powell Jul 3 '19 at 15:20
  • It is likely you will have to write some SQL to match you points of interest to specific line segments to feed to the spanning tree algorithms. But, without any more information, it is hard to tell. – John Powell Jul 3 '19 at 15:41
  • Thank you John, I just edited the post to make it clearer. The pgr_kruskal function is very similar to what I’m looking for, but in this case, it would create the minimum spanning tree for all the intersections of the roads in the roadmap, but I just need to connect the points in the 'Points' layer. – Jony Zambrano Jul 7 '19 at 9:21

The easiest way that I found is using OSMnx library in Python to upload the OSM map as a network in NetworkX (You can also use geopandas), once there is easy to link each of the points in the 'Points' Layer to the closer edge of the OSM layer and then use the steiner_tree algorithm to get the optimal tree, this is an example of a simple implementation of this procedure:

import networkx as nx
import osmnx as ox
from networkx.algorithms import approximation as ap

#Import map
M = ox.graph_from_bbox(37.79, 37.78, -122.41, -122.43, network_type='drive')
M_projected = ox.project_graph(M)

#Create weighted graph from M
G = nx.Graph()
for u,v,data in M_projected.edges(data=True):
    w = data['weight'] if 'weight' in data else 1.0
    if G.has_edge(u,v):
        G[u][v]['weight'] += w
        G.add_edge(u, v, weight=w)

for v,data in M_projected.nodes(data=True):
        if G.has_node(v):
            G.nodes[v]['pos'] = (data['x'],data['y'])    

#Get tree
points = list(G.nodes)[1:20]
route2 = ap.steiner_tree(G, points)

nx.draw(G, pos = pos)
nx.draw(route2, pos = pos)

Sources: https://geoffboeing.com/2016/11/osmnx-python-street-networks/ https://networkx.github.io/documentation/latest/reference/algorithms/generated/networkx.algorithms.approximation.steinertree.steiner_tree.html#networkx.algorithms.approximation.steinertree.steiner_tree https://stackoverflow.com/questions/15590812/networkx-convert-multigraph-into-simple-graph-with-weighted-edges

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.