Joining lines when direction is not known

Question

I am struggling with this problem: I have a polyline made up of several segments. Each segment is made up of points, each point identified in 3D space by coordinates and elevation. Together, if plotted, the segments make up a more or less continuous line (there could be breaks between segments), but the segments themselves are not consecutive, nor do the points of all segments follow the same direction of travel.

The question is: how can I create, using Python preferably, a single line from these non-consecutive segments so that I can measure the distance between points and the total length of the line?

I don't even know which of the segments is the first one or the last one in the line, but somehow I have to put them in the right sequence and make sure they all point in the same direction so I can measure them.
I'll be happy to supply additional info, data, etc... I stress I am not asking for the actual python code (not that I would refuse it...), just the strategy.

Do segments in the polyline cross other segments? If so, is there a requirement to break the segments at these intersections? — Kirk Kuykendall, Commented Mar 3, 2011 at 15:42

Community · Accepted Answer · 2020-06-11 15:27:48Z

The segments can be used to form an abstract graph G in which they play the role of nodes. Consider a node that is the segment (arc) from point P to point Q, PQ. Let R be the closest endpoint among all the other segment endpoints to P and let S be the other endpoint of R's segment. G then contains an edge from node PQ to node RS and we will label this edge with the points P and R.

If we are to be successful, then G is either a linear graph or a single cycle. You can detect which is which by storing the degrees of the nodes as you create the edges. (The degree of a node counts the edges emanating from that node.) All nodes, except possibly two of them, must have degree 2. The other two must both have degree 2 (for a cycle) or degree 1: this marks them as the ends of the polyline. In the first case pick any node to start constructing the polyline; in the second case, pick one of those of degree 1. Any other combination of degrees is an error.

The polyline is initialized to the starting node (an arc). Look at one of the edges e incident on this node: its other node tells you which arc to process next and its label tells you which vertices of those arcs to join. (Join the vertices with a new line segment if they do not have identical coordinates.) Update the growing polyline in this fashion and, at the same time, remove edge e from the graph G. Continue until either an error occurs (and report that the edges do not form a non-branching connected polyline) or all edges are removed. If no error is raised, output the polyline you created.

Example

Sketch of the arcs

In the figure the arcs are AB, CD, EF, and FG. The nodes of the graph are therefore {AB, CD, EF, FG}. The edges are AB--CD (labeled with B and C), CD-EF (labeled with E and F), and EF--FG (labeled with F and F). The degrees of AB and FG are 1 whereas the degrees of CD and EF are 2. Here is a schematic of the abstract graph and its edge labels:

The graph

Let us arbitrarily start with FG, one of the degree-one nodes. Because it has degree 1, there is a unique edge EF--FG connected to it, labeled with F. Initialize the polyline with arc G-->F. Because the label designates a common endpoint of GF and EF, we don't have to make a new connection. Remove edge EF--FG from the graph and extend the polyline with EF via G-->F-->E.

This edge removal reduces the degree of EF from 2 to 1, leaving it with a single edge to arc CD labeled with E and D. This tells you to extend the polyline from E to D (with a new line segment there) and thence along arc CD: G-->F-->E-->D-->C. In the same manner, after removing edge ED--CD, you will extend the polyline further to its final form G-->F-->E-->D-->C-->B-->A. You stop because all edges have been removed from the graph: this indicates the procedure was successful.

Thank you Whuber. I am going to need some time to evaluate your suggestion, but I already have some questions: you say "Let R be the closest endpoint". This I believe is the crux of the problem. There are three potential points of contact between segments, the x , the y and the z coordinates. Even if we eliminate the z coordinate because a line could be perfectly flat, that still leaves two choices. Remember the solution has to be programmed (in Python), no visual choice is possible. — user2117, Commented Mar 3, 2011 at 8:21
My current thinking is to use Matplotlib to plot the segments regardless of sequence, that would create the line but the points would need recreating. How? — user2117, Commented Mar 3, 2011 at 8:31
@Bob Finding the closest point is a relatively simple GIS operation (even in 3D). Each segment PQ has two endpoints P and Q. First find the closest point to P among the collection of endpoints of all segments (not including PQ itself). It is the endpoint, R, of some different segment RS. Create the edge PQ-->RS labeled by P and R. Repeat the search relative to the endpoint Q to obtain the other edge. One thing I did not note: you need a tolerance threshold; if the nearest point is greater than the tolerance, conclude there is no nearest point. — whuber, Commented Mar 3, 2011 at 15:24
The real difficulty in the entire procedure is that connecting disjoint segments can potentially create tiny self-intersections (depending on why and how the segment endpoints don't exactly match). This can be troublesome to fix in general. If such problems occur, consider cleaning up the polyline to eliminate such flaws. — whuber, Commented Mar 3, 2011 at 16:23

geographika · Accepted Answer · 2011-03-03 08:43:27Z

3

Following on from whuber's answer, if you are looking to do this in Python then have a look at the NetworkX library that has all sorts of functionality related to graphs, nodes, and edges. I think it is the Traversal functions that implement what whuber is describing. It also includes basic drawing functionality.

Once you have your line orders then you can easily construct geometries in Shapely to do further GIS type analysis, or display in MatPlotLib, export to GeoJSON etc.

answered Mar 3, 2011 at 8:43

geographika

14.4k4 gold badges55 silver badges79 bronze badges

Add a comment |

johanvdw · Accepted Answer · 2011-03-03 11:50:49Z

2

I would calculate the distance between each start segment and end segment, and then join the ones with the shortest distance between them?

answered Mar 3, 2011 at 11:50

johanvdw

6,24928 silver badges42 bronze badges

I thought of that too. I'm not sure if it would work as a general solution.
– George Silva
Commented Mar 3, 2011 at 13:46
@George You are correct. Problems occur in situations similar to the one I described in a comment to @nathanvda's post. You want to use endpoint-to-endpoint distances for any solution. @johanvdw's suggestion can be implemented if the distances are understood in this sense; it is analogous to the standard technique for constructing the Euclidean minimum spanning tree of a collection of points.
– whuber
Commented Mar 3, 2011 at 15:29
You are right. I meant start node to end node instead of segment.
– johanvdw
Commented Mar 3, 2011 at 16:08
That is indeed what i ment as well :)
– nathanvda
Commented Mar 3, 2011 at 16:20

Add a comment |

nathanvda · Accepted Answer · 2011-03-03 16:18:47Z

You have a lot of segments, no specific order or orientation. You do not know which actually touch or overlap, and which are merely close.

For each segment only the begin and endpoint is important. The goal is to make one big polyline, the orientation of that polyline is not important, i guess.

In that case i would make some kind of set/array of segments, start with the first one, which is entirely random.

In pseudocode do something like

all_segments = set of all segments

# take the first segment out of set
new_polyline = all_segments.pop

until all_segments.empty?
  start_segm = find_segment_closest_to(new_polyline.start_point)
  remove_from_all_segments(start_segm)
  expand_polyline_at_begin(new_polyline, start_segm) 
  end_segm = find_segment_closest_to(new_polyline.end_point)
  expand_polyline_at_end(new_polyline, end_segm)
  remove_from_all_segments(end_segm)
end

Something like that? That is a very high level. You will need to handle boundary cases. I am guessing you know or could the biggest possible gap/distance, because you will need to be able to somehow exclude points being found: if the closest possible point is at the other end of the polyline than it is not an option :) The easiest way to handle that is to define a maximum gap-distance. This will also limit the number of points you will have to look at for eacht segment.

[EDIT: detail the find_segment_closes_to]

To make it absolutely clear how i would find the closes segment, i will write first a very crude approach:

def find_segment_closes_to(point)
  closest_point = nothing
  closest_distance = MAX_GAP_RANGE
  all_segments.each do |segment|
    if distance(segment.end_point, point) < closest_distance
      closest_point = segment.end_point
      closest_segment = segment
      closest_distance = distance(segment.end_point, point)
   else if distance(segment.start_point, point) < closest_distance
      closest_point = segment.start_point
      closest_segment = segment
      closest_distance = distance(segment.start_point, point)
    end       
  end  
  # the closest segment
  return closest_segment
end

This is a very crude approach, that will iterate over all segments and check for each endpoint which is the nearest.

Ideally you would have some datastructure where you could ask for all the start- and endpoints that lie within range and only find the closest point among those.

Does that help? I hope it gets you started.

'find_segment_closest_to' will produce erroneous results in some cases. Imagine a polyline that almost closes on itself, but not quite. One vertex of this polyline might be extremely close to the middle of another of its segments, but it should not be connected to that segment. This is why a correct algorithm relies on point-to-point comparisons, not point-to-segment comparisons. — whuber, Commented Mar 3, 2011 at 15:26
Yes indeed: you should only be looking at begin- and endpoints of a segment. Maybe you said in your solution almost the same, but i found it very hard to read and understand. — nathanvda, Commented Mar 3, 2011 at 16:13
Knowing the language and basic techniques of graph theory is essential if you want to be able to read the literature on algorithms. — whuber, Commented Mar 3, 2011 at 16:25

Stack Exchange Network

Joining lines when direction is not known

4 Answers 4

Example

Your Answer

Hot Network Questions

Joining lines when direction is not known

4 Answers 4

Example

Your Answer

Sign up or log in

Post as a guest

Related

Hot Network Questions