I got a line shapefile and a point shapefile (or a csv file containing coordinates of the point), and I want to clip this line feature by those point in the 'nearest' way.

So far, I have always used edit function in QGIS GUI to clip line feature manually:

enter image description here enter image description here

I am wondering are there QGIS built-in functions could handle this problem (for example, input line feature and point feature, output cut line feature (new file is even better).)? or could I cut line feature in this way using PYQGIS command or user-written python function (even better)?


Update#1

I used PyQGIS code from @xunilk to split line feature showed in figure above, and named my input line and point shapefile as 'my_border' and 'my_point' respectively, but the results looked like:

enter image description here

The result is not what I expect.

Besides, I followed the way in the answer provided by @Micha in splitting a lines layer using points, trying to split line by GRASS.

Because my point is not on the line, I have to run the following module to get 'error point':

v.distance
v.patch
v.clean

I think there are no problem use the first 2 module, but after I used v.clean:

v.clean --overwrite input=tcall@nl output=tcerrpoint tool=break 

I just got a line feature, not error points:

enter image description here

Were there something I missed or did wrong in the procedure?


Update#2

I used the updated code provided by @xunilk, which trying to handle multiple cut points. I still used my line and point shapefile as input layers, and the following figures are their attribute tables.

enter image description here enter image description here

The output is:

enter image description here

Still not what I expected, I am wondering that maybe it outputted strange line feature because the input line has too many segments.

If the PyQGIS script being modified to capable of handling selected feature of input line shapefile, and before run the script we select the nearest line segment to one point, maybe the output can be right?


Update#3

This time I edited input line shapefile 'my_border' to a line with only one attribute:

enter image description here

and used PyQGIS code in the edited answer of @xunilk:

registry = QgsMapLayerRegistry.instance()

my_border = registry.mapLayersByName('my_border')
my_point = registry.mapLayersByName('my_point')

feats_my_border = [ feat for feat in my_border[0].getFeatures() ]
feats_my_points = [ feat for feat in my_point[0].getFeatures() ]

points = [ feat.geometry().asPoint() for feat in feats_my_points ]

my_tuplas = []

for point in points:
    tupla = feats_my_border[0].geometry().closestSegmentWithContext(point)
    my_tuplas.append(tupla)

my_line = [ feat.geometry().asPolyline() for feat in feats_my_border ]

new_line = []

k = 0

try:
    for i, line in enumerate(my_line):
        for j, point in enumerate(line):
            new_line.append(point)
            if j == my_tuplas[k][2] - 1:
                new_line.append(my_tuplas[k][1])
                k += 1

except IndexError:
    pass

k = 0

lines = [ [] for i in range(len(points) + 1) ]

epsg = my_border[0].crs().postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'split_line',
                           'memory')

prov = mem_layer.dataProvider()

try:
    for point in new_line:
        lines[k].append(point)
        if point == my_tuplas[k][1]:
            k += 1
            lines[k].append(point)

except IndexError:
    pass

for i, line in enumerate(lines):

    feat = QgsFeature()
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolyline(line))
    prov.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

This time the attribute of output 'split_line' looked rational because I inputted 8 point to cut line:

enter image description here

But no matter what I tried including set layer CRS, zoom to layer, show in overview, etc., I could not see this output file. Besides, I saved this output file as a new ESRI shapefile and checked its size, which was 208 kb (input 'my_border.shp' is 21 kb), all seemed to be rational, but I just can't see the output.

  • 1
    I just got a message from my friend that he wrote a Arcpy script which could cut line by point, but I want to do this with open-sources... – Heinz Oct 20 '16 at 3:12

You can use 'closestSegmentWithContext' from QgsGeometry. The argument is the cut nearest QgsPoint and, it produces a tuple with the squared cartesian distance from closest segment to polyline, the nearest point on the segment and, finally, the index of the vertex (at polyline) after the closest segment.

To illustrate the approach, I used next shapefiles for only cutting in two features the polyline layer.

enter image description here

PyQGIS code is:

registry = QgsMapLayerRegistry.instance()

my_border = registry.mapLayersByName('my_border')
my_point = registry.mapLayersByName('my_point')

feats_my_border = [ feat for feat in my_border[0].getFeatures() ]

feats_my_point = [ feat for feat in my_point[0].getFeatures() ]

my_tupla = feats_my_border[0].geometry().closestSegmentWithContext(feats_my_point[0].geometry().asPoint())

my_line = [ feat.geometry().asPolyline() for feat in feats_my_border ]

first_segment = []
second_segment = [my_tupla[1]]

for line in my_line:
    for i, point in enumerate(line):
        if i < my_tupla[2]:
            first_segment.append(point)
        else:
            second_segment.append(point)

first_segment.append(my_tupla[1])

epsg = my_border[0].crs().postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'split_line',
                           'memory')

prov = mem_layer.dataProvider()

feat = QgsFeature()

feat.setAttributes([0])

feat.setGeometry(QgsGeometry.fromPolyline(first_segment))

prov.addFeatures([feat])

feat.setAttributes([1])

feat.setGeometry(QgsGeometry.fromPolyline(second_segment))

prov.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

After running it, at the Python Console of QGIS, I got:

enter image description here

It can be observed the first selected feature at the created memory layer.

Editing Note:

When there are more than one cut point, next code works well:

registry = QgsMapLayerRegistry.instance()

my_border = registry.mapLayersByName('my_border')
my_point = registry.mapLayersByName('my_point')

feats_my_border = [ feat for feat in my_border[0].getFeatures() ]
feats_my_points = [ feat for feat in my_point[0].getFeatures() ]

points = [ feat.geometry().asPoint() for feat in feats_my_points ]

my_tuplas = []

for point in points:
    tupla = feats_my_border[0].geometry().closestSegmentWithContext(point)
    my_tuplas.append(tupla)

my_line = [ feat.geometry().asPolyline() for feat in feats_my_border ]

new_line = []

k = 0

try:
    for i, line in enumerate(my_line):
        for j, point in enumerate(line):
            new_line.append(point)
            if j == my_tuplas[k][2] - 1:
                new_line.append(my_tuplas[k][1])
                k += 1

except IndexError:
    pass

k = 0

lines = [ [] for i in range(len(points) + 1) ]

epsg = my_border[0].crs().postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'split_line',
                           'memory')

prov = mem_layer.dataProvider()

try:
    for point in new_line:
        lines[k].append(point)
        if point == my_tuplas[k][1]:
            k += 1
            lines[k].append(point)

except IndexError:
    pass

for i, line in enumerate(lines):

    feat = QgsFeature()
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolyline(line))
    prov.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

After running above code at the Python Console of QGIS, I got:

enter image description here

At the image, it was selected only the second feature of polyline memory layer.

  • @Heinz You're welcome. I also modified my answer. It is included the possibility of more than one cut point. – xunilk Oct 19 '16 at 20:28
  • 1
    I edited my post again! – Heinz Oct 20 '16 at 1:27
  • 1
    I tried using this code to split a polyline using a point file containing multiple (156) points. It created the right number of features, but only split the line up into a 4 segments. The attribute table has 157 features, but all but the first 4 features have no geometry. Do you have any idea what is going on, or how I could modify the code so it continues to split the original line at the points? Thanks! – Brian Fisher Mar 13 '17 at 19:23

I know this wasn't what you asked for but its possible in postGIS. There is ST_LineLocatePoint:

If you have an existing point table that you want to reference to a linear network, use the ST_LineLocatePoint function, which takes a line and point, and returns the proportion along the line that the point can be found.

The use ST_Split

I tried to make such a tool but without the elegant use of line.closestSegmentWithContext(point).

You can find it here: Plugin repository on github

Please, see this plugin attempt to be used with a QGIS 2.18.x version, no lesser one.

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.