Shapely multilinestring: depth first search

Question

I have a Shapely multiline string that has multiple branches that I want removed. I've attached an image where the red is the main line that I want to keep and the black are to be removed.

I can't filter based on individual linestring lengths as some segments that need to be kept are shorter than the branches.
I was thinking a Depth-First Search would be appropriate but I'm unsure how to implement it when working with a shapely object. As well as calculate the full length of the main line.
Any tips would be appreciated.

Answer 1

You can use networkx:

Node the multiline to insert vertices at all line intersection, explode it into individual line segment.

Then iterate over each line and build a graph using start and end coordinates of the line as nodes, and with line length and line geometry as edge attributes.

Then find the longest shortest path between all end nodes.

from shapely.wkt import loads
import matplotlib.pyplot as plt
from shapely.plotting import plot_line
from shapely import union_all
import networkx as nx
from itertools import combinations

#Create a multiline from a wkt string and plot it
wkt_string = r"MULTILINESTRING((614330.676566 6826096.386227, 614343.019894 6826101.20304), (614333 6826179, 614350.427132 6826138.336691), (614341.294058 6826089.3079, 614343.019894 6826101.20304), (614343.019894 6826101.20304, 614363 6826109), (614361.167824 6826156.27977, 614346.41389 6826147.700924), (614363.073528 6826150.834902, 614348.455715 6826142.936665), (614363.345771 6826145.662278, 614350.427132 6826138.336691), (614350.427132 6826138.336691, 614363 6826109), (614363 6826109, 614433 6826089, 614471.289314 6826121.819412), (614434.673536 6826114.626534, 614462.993898 6826114.709055), (614434.945779 6826121.704861, 614471.289314 6826121.819412), (614466.670785 6826117.860673, 614474.965555 6826099.108661, 614476.054529 6826112.993073), (614471.289314 6826121.819412, 614471.289314 6826121.819412), (614471.289314 6826121.819412, 614503 6826149, 614526.419553 6826149), (614526.419553 6826149, 614525.602823 6826168.530722), (614526.419553 6826149, 614531.047691 6826149), (614531.047691 6826149, 614531.047691 6826129.599919), (614531.047691 6826149, 614537.581532 6826149), (614537.309288 6826166.080531, 614533.770124 6826169.619695), (614537.309288 6826166.080531, 614541.937426 6826169.347452), (614537.581532 6826149, 614537.309288 6826166.080531), (614537.581532 6826149, 614583 6826149, 614583 6826109))"
multiline = loads(wkt_string)
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(20,12))
_ = plot_line(line=multiline, ax=ax[0][0], add_points=False, color="green")
ax[0][0].set_title("Original multiline")

#Node the lines / insert vertices att all junctions
unioned = union_all(geometries=multiline, grid_size=1)
line_segments = list(unioned.geoms) #Explode the unioned multiline into individual line segments which are split at every line junction

#Build a graph
G = nx.Graph()
for line in line_segments: #For each split line
    coords = list(line.coords) #Extract its coordinates
    #Add an edge using the start and end coordinates as nodes, and append 
    #  line length and the line geometry as edge attributes
    G.add_edge(coords[0], coords[-1], length=line.length, geometry=line)
nx.draw(G, ax=ax[0][1])

#Find the longest "shortest path" between all pairs of line ends

#List all start or end nodes
end_nodes = [node for node, d in G.degree if d==1]
shortest_paths = [] #A list of all shortest paths
for n1, n2 in combinations(end_nodes, 2):
    shortest_path = nx.single_source_dijkstra(G, source=n1, target=n2, weight="length")
    shortest_paths.append(shortest_path)

#Sort the list by length so the longest paths is last [-1]
longest_shortest_path = sorted(shortest_paths, key=lambda x: x[0])[-1][1] #[1] is to extract only the path itself, without the length

#Iterate over each pair of nodes and extract the line geometries that are saved in the graph edges
lines = []
for n1, n2 in zip(longest_shortest_path, longest_shortest_path[1:]):
    lines.append(G.edges[n1, n2]["geometry"])

for line in lines:
    _ = plot_line(line=line, ax=ax[1][0], color="red", add_points=False)
ax[1][0].set_title("The longest shortest path between all end nodes")