Cutting line feature by the nearest point in QGIS

Question

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:

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:

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:

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.

The output is:

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:

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:

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.

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, Commented Oct 20, 2016 at 3:12

Taras · Accepted Answer · 2021-12-07 08:04:23Z

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.

PyQGIS 2 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:

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

PyQGIS 3 code is:

registry = QgsProject.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().asMultiPolyline() for feat in feats_my_border ]

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

for line in my_line[0]:
    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.fromPolylineXY(first_segment))

prov.addFeatures([feat])

feat.setAttributes([1])

feat.setGeometry(QgsGeometry.fromPolylineXY(second_segment))

prov.addFeatures([feat])

QgsProject.instance().addMapLayer(mem_layer)

Editing Note:

When there are more than one cut point, next code (PyQGIS 2) 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:

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

With QGIS 3 following code also works:

registry = QgsProject.instance()

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

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().asMultiPolyline()[0] 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.fromPolylineXY(line))
    prov.addFeatures([feat])

QgsProject.instance().addMapLayer(mem_layer)

@Heinz You're welcome. I also modified my answer. It is included the possibility of more than one cut point. — xunilk, Commented Oct 19, 2016 at 20:28
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, Commented Mar 13, 2017 at 19:23
I tried this code but it doesn't seem to work now. I found out that the QgsMapRegistry has now been shifted to QgsProject. what changed needs to be done to his code to get the exact output ? — Ganesh Ambi Ramakrishnan, Commented Apr 14, 2020 at 3:58
@GaneshAmbiRamakrishnan I updated code for QGIS 3 (with one point) and it produces same output. — xunilk, Commented Apr 14, 2020 at 4:38

Knightshound · Accepted Answer · 2016-10-19 09:13:59Z

1

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

answered Oct 19, 2016 at 9:13

Knightshound

4,48019 silver badges39 bronze badges

Add a comment |

Babel · Accepted Answer · 2021-12-13 10:42:30Z

Based on the point layer, create auxiliary lines that connect each point with the closest point on the line layer, then extend this auxiliary lines so that they cross the line you want to split.

Use QGIS expressions with Geometry by expression and this expression (replace line on line 7 with the name of your line layer and adapt the value of 100 in the 2nd last line to your needs: the distance the auxiliary lines should be extended. Choose a small value to avoid auxiliary lines crossing several lines.:
```
 extend (
     make_line (
         $geometry, 
         closest_point (
             array_first (
                 overlay_nearest (
                     'line', 
                     $geometry
                 )
             ), 
             $geometry
         )
     ),
     0,
     100
 )
```
Use Menu Processing / Toolbox / Split with lines to cut the initial lines with the auxiliary lines created in step 1.

Initial line layer (black line), points (red dots) with the auxiliary lines in red. The colored lines show the new line features after the split operation:

Matt · Accepted Answer · 2022-12-20 13:25:28Z

I find another solution using PyQgis, from this post:

Split 3D polyline at point using QGIS

The code could be this:

def split_line(pointfeat, polyfeat):
    point = pointfeat.asPoint()
    poly = polyfeat.asPolyline()
    geo = QgsGeometry().fromPolylineXY(poly)
        
    (sqdist, projpoint, aftervertex, leftOf) = geo.closestSegmentWithContext(point)
    newline1 = poly[:aftervertex]  # from start to break point
    newline1.append(projpoint)  # append point
    newline2 = poly[aftervertex:]  # from break point to end
    newline2.insert(0, projpoint)  # insert at position 0
            
    newline1_Polyline=QgsGeometry.fromPolylineXY(newline1) # Convert list of Points in a Polyline
    new_line2_Polyline=QgsGeometry.fromPolylineXY(newline2)
        
    return new_line1_polyline, new_line2_polyline
    
    
# Load layer line to split
layer_line=QgsVectorLayer()
    
# Load layer rutometros for save partial Track
layer_points=QgsVectorLayer()
first=True
filter=sql expression to filter points
layer_points.selectByExpression(filter,QgsVectorLayer.SetSelection)
    
for feature in layer_points.selectedFeatures():
    geom_point=feature.geometry()
    if not first:
        new_line,geom_line=split_line(geom_previous,geom_line)
        ## add new line in a QgsVectorLayer
            
    geom_previous=geom_point
    first=False

Jean-Christophe Baudin · Accepted Answer · 2020-02-26 22:58:30Z

0

By now i made a new version of ClosestPoint plugin: version 4.0. for QGIS 3. It's directly available in Qgis plugin repository. This version is experimental so check that experimental plugins are aloud in your QGIS parameters.

This 4.0 version of ClosestPoint plugin does cut lines by nearest points.

It can deal with numerous points projected on a same vertex. This case wasn't delt with perfectly by previous code solutions in this post.

answered Feb 26, 2020 at 22:58

Jean-Christophe Baudin

572 bronze badges

Add a comment |

aldo_tapia · Accepted Answer · 2017-08-09 16:01:35Z

-1

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.

edited Aug 9, 2017 at 16:01

aldo_tapia

13.8k5 gold badges32 silver badges58 bronze badges

answered Aug 9, 2017 at 15:41

Jean-Christophe Baudin

572 bronze badges

Add a comment |

Stack Exchange Network

Cutting line feature by the nearest point in QGIS

6 Answers 6

PyQGIS 2 code is:

PyQGIS 3 code is:

Your Answer

Not the answer you're looking for? Browse other questions tagged
qgis
pyqgis
point
line
splitting
or ask your own question.

Linked

Hot Network Questions

Cutting line feature by the nearest point in QGIS

6 Answers 6

PyQGIS 2 code is:

PyQGIS 3 code is:

Your Answer

Sign up or log in

Post as a guest

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

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
qgis
pyqgis
point
line
splitting
or ask your own question.