Bring a list of unsorted linestrings in order so that they form a path from origin to destination

Question

I have a table of linestrings together with some attributes that, when displayed all at the same time, represent a connected path from an origin to a destination (each linestring has one predecessor, except for the first one, and one successor, except for the last one). However, the feature order in my table is randomized.

For another application, I need the features in the table to be ordered from origin to destination (or vice versa), so that for each row/linestring in the table, the preceding row represents the current row's predecessor linestring and the succeeding row represents the current row's successor linestring.

I am in control over the data, so writing them out in any format suitable for this operation should be possible (csv, geojson, shapefile, ...).

I have QGIS, ArcMap, PostGIS and Python easily available but any other solution will also be fine.

As an example that everybody can easily visualize, take the following GeoJSON:

{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "someAttribute": 2 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88175, 52.294 ], [ 6.88198, 52.29404 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 8 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88945, 52.29465 ], [ 6.88955, 52.29465 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 7 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88934, 52.29465 ], [ 6.88945, 52.29465 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 4 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88294, 52.29418 ], [ 6.8857, 52.29448 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 3 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88198, 52.29404 ], [ 6.88294, 52.29418 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 5 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.8857, 52.29448 ], [ 6.88836, 52.29462 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 9 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88955, 52.29465 ], [ 6.88963, 52.29466 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 6 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88836, 52.29462 ], [ 6.88934, 52.29465 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 11 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.89216, 52.29463 ], [ 6.89285, 52.2946 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 12 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.89285, 52.2946 ], [ 6.89528, 52.29444 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 1 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88108, 52.29389 ], [ 6.88175, 52.294 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 13 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.89528, 52.29444 ], [ 6.89794, 52.29414 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 10 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88963, 52.29466 ], [ 6.89216, 52.29463 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 14 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.89794, 52.29414 ], [ 6.90067, 52.2937 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 16 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.9034, 52.29311 ], [ 6.90369, 52.29304 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 15 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.90067, 52.2937 ], [ 6.9034, 52.29311 ] ] } }
]
}

For clarification, if you visualize this (e.g. with QGIS), you get the following picture:

The single linestrings form a connected path from 1 to 16, even though they are unordered in the data source. What I want is to sort my data source accordingly. The enumeration in someAttribute of course is only for clarification and not given in the real data, so simply sorting by attribute is not an option ;-)

Answer 1

First make sure your data is in an editable format. GeoJSON is not editable in QGIS. If the Editing button is disabled then this won't work - save the layer as a Shapefile or some other layer with enabled editing.

Then install the "Multiline Join" plugin.

Select the feature you want to merge, and run the plugin (from the "Vector"->"Join Multiple Lines" menu). Your lines will be merged, and your layer will be in edit mode. Click the pencil to commit the edits.

You should notice your sample data now has one feature, rather than 16!

Finding the sequence of lines to traverse that form the long line can be done using R and the igraph package.

library(igraph)
library(rjsonlite) # for reading data
j=fromJSON("locs.json", flatten=TRUE)
pts = j$features$geometry.coordinates

pts is now a list of 2x2 matrices:

> str(pts)
List of 16
 $ : num [1:2, 1:2] 6.88 6.88 52.29 52.29
 $ : num [1:2, 1:2] 6.89 6.89 52.29 52.29
 $ : num [1:2, 1:2] 6.89 6.89 52.29 52.29
 ... etc ...

Then this next bit of code turns the coordinates into an "edge list". You might be able to start here if you can output your data in this format. Up to this point all I've done is massage your JSON:

> EL = t(sapply(pts,function(m){c(paste(m[1,],collapse=","), paste(m[2,], collapse=","))}))
> head(EL)
     [,1]               [,2]              
[1,] "6.88175,52.294"   "6.88198,52.29404"
[2,] "6.88945,52.29465" "6.88955,52.29465"
[3,] "6.88934,52.29465" "6.88945,52.29465"
[4,] "6.88294,52.29418" "6.8857,52.29448" 
[5,] "6.88198,52.29404" "6.88294,52.29418"
[6,] "6.8857,52.29448"  "6.88836,52.29462"

Now each row has a text-encoded coordinate pair. We're going to use those to identify "nodes" (aka vertices) in a graph, so that each row connects two "nodes" and makes an "edge". This relies on the text representation being identical, so any rounding errors will screw up the process. EL[1,2] is identical to EL[5,1] and so that will produce a connection in the following graph:

g = graph.edgelist(EL, directed=FALSE)
plot(g)

The plot will show you it has made all the connections. The layout is arbitrary and not related to your coordinates.

Now we need to find the start and end points. In the jargon, these are nodes with degree 1, and all the other nodes should have degree 2 (one line in, one line out).

> degree(g)
  6.88175,52.294 6.88198,52.29404 6.88945,52.29465 6.88955,52.29465 
               2                2                2                2 
6.88934,52.29465 6.88294,52.29418  6.8857,52.29448 6.88836,52.29462 
               2                2                2                2 
6.88963,52.29466 6.89216,52.29463  6.89285,52.2946 6.89528,52.29444 
               2                2                2                2 
6.88108,52.29389 6.89794,52.29414  6.90067,52.2937  6.9034,52.29311 
               1                2                2                2 
6.90369,52.29304 
               1 
> which(degree(g)==1)
6.88108,52.29389 6.90369,52.29304 
              13               17 

That looks good. Now we just need to find the shortest path (the only path!) between node 13 and 17, and spit out the edges that make that path:

> get.shortest.paths(g, V(g)[13], V(g)[17], output="epath")
$vpath
NULL

$epath
$epath[[1]]
 [1] 11  1  5  4  6  8  3  2  7 13  9 10 12 14 16 15

So we start with edge 11, then 1, then 5, then 4 and so on. When I look these up in your original data attribute I see...

> j$features$properties.someAttribute[get.shortest.paths(g, V(g)[13], V(g)[17], output="epath")$epath[[1]]]
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16

Which I think is correct.

This process can go wrong in a few places, most likely the construction of the graph if the points don't exactly match up. If you have a graph where the node degrees are all 2 except for two 1-degrees, then you don't have a single line segment - you've either got disjoint segments or a tree.