Efficient algorithms to split and join lines?

Question

Referring to the Figure below:

1. How to efficiently join segments into lines?
2. How to efficiently break lines into segments?

The above will be applied million times.

Answer 1

You can use effectively Shapely, and Fiona to read a shapefile for example:

import fiona
# open a line shapefile
file = fiona.open('lines.shp')
# first element of the shapefile
first = file.next
print first
{'geometry': {'type': 'LineString', 'coordinates': [(203317.23, 90448.75), (203679.62, 90105.68), (203882.57, 89902.74), (204143.49, 89641.81), (204394.75, 89385.72), (204563.87, 89235.93)]}, 'id': '0', 'properties': {u'id': "1"}}

Now import Shapely

from shapely.geometry import Point, LineString, shape
geom = shape(first['geometry'])
# now it is a Shapely geometry
print geom
LINESTRING (203317.23 90448.75, 203679.62 90105.68, 203882.57 89902.74, 204143.49 89641.81, 204394.75 89385.72, 204563.87 89235.93)

Splitting

We can split the line in segments, each pair of coordinates (Point) define a line segment:

def pair(list):
    '''Iterate over pairs in a list -> pair of points '''
    for i in range(1, len(list)):
        yield list[i-1], list[i]

for seg_start, seg_end in pair(geom.coords):
    line_start = Point(seg_start)
    line_end = Point(seg_end)
    segment = LineString([line_start.coords[0],line_end.coords[0]])
    print segment
LINESTRING (203317.23 90448.75, 203679.62 90105.68)
LINESTRING (203679.62 90105.68, 203882.57 89902.74)
LINESTRING (203882.57 89902.74, 204143.49 89641.81)
LINESTRING (204143.49 89641.81, 204394.75 89385.72)
LINESTRING (204394.75 89385.72, 204563.87 89235.93)

You could do it directly with Fiona:

for seg_start, seg_end in paires(first['geometry']['coordinates']):
    ....

Union

1. If the line is straight, simply use the first and the last point to make the LineString

   LINESTRING (203317.23 90448.75, 204563.87 89235.93)

2. In other cases, with shapely, you can use union, cascaded_union(geoms) or unary_union(geoms) and the result is a MULTILINESTRING;

With union:

line = LineString()
for element in line['geometry']['coordinates']:
   geom =shape(line['geometry'])
   profile = profile.union(geom)
print line
MULTILINESTRING ((203317.23 90448.75, 203679.62 90105.68), (203679.62 90105.68, 203882.57 89902.74), (203882.57 89902.74, 204143.49 89641.81), (204143.49 89641.81, 204394.756 89385.72), (204394.75 89385.72, 204563.87 89235.93))

or directly with shapely:

line = LineString()
for seg_start, seg_end in paires(geom.coords):
     line = line.union(LineString([Point(seg_start).coords[0],Point(seg_end).coords[0]]))

Intersection

The shapely function is geom1.intersection(geom2)

1. for two straight lines, it is simply

   line1.intersection(line2) # the result is a Point

2. for two lines with segments, it is

   (union of segments of line1).intersection(union of segments of line2) # the result is a Point or a MultiPoint (more than a single intersection)

And saving the results is easy with Fiona

Answer 2

It seems that Shapely can answer what you are expecting. See the manual. Create lines using native objects. For what you need, only loop ;) to get segments of a line. For getting intersection point, use operators like intersection

I can't say that is efficient but I suppose that because Shapely uses Geos behind, it's quite fast.

Answer 3

To join line segments you can assign a from-node# and a to-node# to each line segment, where the node# is a hash of x and y rounded off. For example, a point at 123.45,567.8 might be computed to be 123568. This implies a tolerance of 1. You could generalize this to round off to the closest multiple of the tolerance.

Once you have From and To node#'s, you can topologically sort the segments. This can be done by building a valence count for each node with a dictionary where the key is node# and the value is valence (the number of segments that share the node). Start at a node whose valence == 1 and walk until you reach another node whose valence==1, connecting the dots in other words. You can then use the topologically sorted list to build an ordered list of coordinates, i.e. a polyline.

Update: See whuber's similar question here.