8

I have a set of lines representing bus lines. Some of the lines are overlapping and take the same roads.

enter image description here

I am able to extract the nodes. enter image description here

However I am interested in extracting only crossings like this: enter image description here

How can I do this? I am looking for ways with QGIS or Python.

I tried the intersection method from GDAL Python but this basically returns me only the vertices.

The Line Intersections method from QGIS returns me the crossings if two lines cross. However in the case that two bus lines go far part of their route on the same road, it doesn't give me they point where they merge.

  • Have you tried the line intersection tool in QGIS: Vector Analysis Tool > Line Intersections... It won't give you the end and starting nodes of a line, but all intersections. – Jakob Mar 6 '15 at 14:17
  • Yeah I wrote about this in the question. – ustroetz Mar 6 '15 at 14:52
  • I'm unclear on what you're asking, in part because all lines are symbolized the same way in your images - I can't tell different routes apart to understand what nodes you're looking at or why there are so many in the second image. Are the routes coincident on roads? Are they all two-point line segments, or continuous polylines? I note the behavior you describe is the same as ArcGIS' Intersect tool - lines / lines with output of lines gives you overlap, but lines / lines with point output only gives crossings. – Chris W Mar 6 '15 at 22:04
  • Based on that, to get what I think you want, you might have to use both methods. Get the crossings (line/line=point) and then get the overlaps (line/line=line) and extract the start/end nodes for those overlap lines. Those should be all the points/nodes you're looking for. – Chris W Mar 6 '15 at 22:12
16

The nodes:

You want two things, the end points of the polylines (without intermediate nodes) and the intersection points. There are an additional problem, some polylines end points are also intersection points:

enter image description here

A solution is to use Python and the modules Shapely and Fiona

1) Read the shapefile:

from shapely.geometry import Point, shape
import fiona
lines = [shape(line['geometry']) for line in fiona.open("your_shapefile.shp")]

2) Find the end Points of the lines (how would one get the end points of a polyline?):

endpts = [(Point(list(line.coords)[0]), Point(list(line.coords)[-1])) for line  in lines]
# flatten the resulting list to a simple list of points
endpts= [pt for sublist in endpts  for pt in sublist] 

enter image description here

3) Compute the intersections (iterating through pairs of geometries in the layer with the itertools module). The result of some intersections are MultiPoints and we want a list of points:

import itertools
inters = []
for line1,line2 in  itertools.combinations(lines, 2):
  if  line1.intersects(line2):
    inter = line1.intersection(line2)
    if "Point" == inter.type:
        inters.append(inter)
    elif "MultiPoint" == inter.type:
        inters.extend([pt for pt in inter])
    elif "MultiLineString" == inter.type:
        multiLine = [line for line in inter]
        first_coords = multiLine[0].coords[0]
        last_coords = multiLine[len(multiLine)-1].coords[1]
        inters.append(Point(first_coords[0], first_coords[1]))
        inters.append(Point(last_coords[0], last_coords[1]))
    elif "GeometryCollection" == inter.type:
        for geom in inter:
            if "Point" == geom.type:
                inters.append(geom)
            elif "MultiPoint" == geom.type:
                inters.extend([pt for pt in geom])
            elif "MultiLineString" == geom.type:
                multiLine = [line for line in geom]
                first_coords = multiLine[0].coords[0]
                last_coords = multiLine[len(multiLine)-1].coords[1]
                inters.append(Point(first_coords[0], first_coords[1]))
                inters.append(Point(last_coords[0], last_coords[1]))

enter image description here

4) Eliminate duplicates between end points and intersection points (as you can see in the figures)

result = endpts.extend([pt for pt in inters  if pt not in endpts])

5) Save the resulting shapefile

from shapely.geometry import mapping
# schema of the shapefile
schema = {'geometry': 'Point','properties': {'test': 'int'}}
# creation of the shapefile
with fiona.open('result.shp','w','ESRI Shapefile', schema) as output:
    for i, pt in enumerate(result):
        output.write({'geometry':mapping(pt), 'properties':{'test':i}})

Final result:

enter image description here

The line segments

If you want also the segments between the nodes, you need to "planarize" (Planar Graph, no edges cross each other) your shapefile. This can be done by the unary_union function of Shapely

from shapely.ops import unary_union
graph = unary_union(lines)

enter image description here

  • Thanks @gene for the detailed answer. I edited the part where it goes over the different geometrie types. In my case the intersection also returns lines, geometrycollections, etc. But this depends on the input data. I was not clear enough in my question. – ustroetz Mar 9 '15 at 12:56
  • Great answer. Might I add that it is not necessary to do the following: result = endpts.extend([pt for pt in inters if pt not in endpts]) since it appears that the .extend method modifies endpt. In my case result = None after that operation. It is endpts which ends up containing the result sett – user32882 Dec 13 '18 at 9:29

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.