5

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.

enter image description here

  • 1
    What GIS software is involved? If none, can you say so and give more details of what you are looking for here, please? – PolyGeo May 22 '13 at 6:18
  • Are the lines and the points (where split has to take place) topologically connected? Do you know which points are present on a given line? – Devdatta Tengshe May 22 '13 at 6:29
  • @PolyGeo None. We will code in Python. The points above are intersection points with many other lines that are removed later. So our interest is to split each line to many segments based on the intersecting points (part left). Our another interest is reversing the act, that is having a series of aligned (parallel) segments to join them as one line if possible (part right). – Developer May 22 '13 at 7:09
  • @DevdattaTengshe Right, even after splitting the will have points in common. The gaps above (part left) are only for demonstration. We don't have a record of belonging points to a specific line. Indeed, the inputs are segments (or lines for question 2) and points. We drew them here together for demonstration. – Developer May 22 '13 at 7:13
  • Maybe wrong but why don't you use a spatial DB like PostGIS or Spatialite. More efficient if you use spatial index on your tables and low level. Built-in operations for your need FYI – ThomasG77 May 23 '13 at 6:23
8

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

  • 1
    Please note that we don't know which point belongs to which line. Line segments are not listed in any order. They are randomly located in an array. If you know which line is aligned (next) to a line, joining them is not an issue. But notice that we don't have any topology information of lines. No information between points and lines. – Developer May 23 '13 at 5:56
1

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.

  • We are not sure about that either. In our previous attempts for another issue we noticed however Shapely was 10 times slower than our mixed Python-Fortran implementation. – Developer May 23 '13 at 6:05
0

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.

  • If you wish to go further, it is the mathematical area of Planar Graphs (link) implemented in OpenJump (link – gene May 22 '13 at 18:56
  • Indeed, we have already implemented a similar idea in Python using defaultdict. This is not efficient as you need to update the dictionary entirely after joining every two segments. – Developer May 23 '13 at 6:03
  • @gene I don't think this approach requires the graph to be planar, does it? – Kirk Kuykendall May 23 '13 at 13:35
  • @Developer I wouldn't join two segments at a time, but rather just list the coordinates for each chain of edges after sorting. – Kirk Kuykendall May 23 '13 at 13:37

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.