[GIS] Selecting/tracing part of line feature using QGIS

lineqgisselectshapefiletrace

I have a shapefile which comprises a number of lines, each representing a bus route. I'm wanting to create a new bus route which is the same as an existing line to a certain point, but then branches off and follows a different route (see image below).

What's the best way to do this?

I'll be doing it numerous times. I'm thinking the best way is to draw a polygon around the parts I want and then clip using the same polygon, and dissolving the result. I don't really want to draw new lines by snapping the vertices to existing routes as each route has thousands of vertices.

I am using QGIS 3.2.3.

Best Answer

You can use the Snapping toolbar and click on Enable tracing to trace the existing polylines this much more easier to draw a new line using existing lines'

Existing line'

Adding new line (purple color)

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

ArcGIS Desktop – Extending Lines to Points Without Altering Original Geometry

This is what snap does:

As one can see it is not what @Stella wants

The following workflow based on this structure of original pipes table:

arcpy.FeatureVerticesToPoints_management("ORIGINAL","D:/Scratch/ENDS.shp", "BOTH_ENDS")
arcpy.AddField_management("ENDS","WHAT", "SHORT")
arcpy.CalculateField_management("ENDS", "WHAT", "!FID! % 2", expression_type="PYTHON_9.3")
arcpy.AddField_management("ENDS", "UNIQUE","TEXT", field_length="10")
arcpy.CalculateField_management("ENDS", "UNIQUE", expression="'%s %s' %( !LINE_ID!, !WHAT!)", expression_type="PYTHON_9.3")
arcpy.Near_analysis("ENDS","NODE", search_radius="", location="LOCATION")
arcpy.TableToTable_conversion(in_rows="ENDS", out_path="D:/Scratch", out_name="xy.dbf")

ADD XY data using NEAR_X and NEAR_Y fields from XY.dbf and convert them to DOUBLES.shp (no clue what ESRI did to this tool)

arcpy.Merge_management("DOUBLES;ENDS","D:/Scratch/MERGED.shp")
arcpy.PointsToLine_management("MERGED", "D:/Scratch/connections.shp", Line_Field="UNIQUE")

RESULTS:

Because connections table contains parent pipe name in field UNIQUE:

it can be derived into a new field:

Afterwards connections can be merged with PIPES and result dissolved using LID. This is no brainer to bring attributes, e.g. DIAMETER, back from PIPES table.

Dissolving however will modify geometry, because without it you task cannot be accomplished.

Best Answer

Related Solutions

QGIS Python – How to Intersect Lines to Get Crossings

The nodes:

The line segments

ArcGIS Desktop – Extending Lines to Points Without Altering Original Geometry

Related Question