QGIS – Deleting Fixed Distance of a Line from a Point on the Line

linepointqgis

I have multiple points on lines in QGIS. Starting at every point I want to delete 50m of the line for each point. The direction does not matter.

If two points are too close the 50m have to delete just once.

My problem is besides how do it, what to do it without touching the lines near by. Therefore buffering is not possible. Also a problem is multiple points are close to each other on the same line.

I though about the unconventional way of creating a 50 m x 1 m polygon and put on the point (ca. 7800…)

How could I do this?

Best Answer

Use network tool: Menu processing > Toolbox > Service Area (from layer). Use the rails as network layer and the points as Vector layer with start points. Set a distance of 25 meters. The tool will find the nearest reailway line of the point and from there create a line along the network of all points within 25 m to both sides (thus 50 m in total).
Apply a very small buffer (like 0.1 meters) around the result.
Run Menu Vector > Geoprocessing tool > Difference with rail as Input layer and the buffer from step 2 as Overlay layer.

And here you are with the result:

Blue point: initial point; black lines: rails; yellow: service area; red: rails with deleted section:

Related Solutions

[GIS] Calculating distance from point to polyline

Script below designed to run from mxd. It assumes that you have empty a table (“nearLines”) to populate in mxd:

Where pointID and lineID are fields to store OIDs of input layers (type long), Distance field type double.

import arcpy
# parameters to re-type
maxDistance=100
mxd = arcpy.mapping.MapDocument("CURRENT")
# get lines
lines = arcpy.mapping.ListLayers(mxd,"lines")[0]
d=arcpy.Describe(lines)
fidLine = d.OIDFieldName
# get points
points = arcpy.mapping.ListLayers(mxd,"points")[0]
d=arcpy.Describe(points)
fidPoint = d.OIDFieldName
table = arcpy.mapping.ListTableViews(mxd,"nearLines")[0]
# process
curT=arcpy.da.InsertCursor(table,("POINTID","LINEID","DISTANCE"))
nodesDict={}
with arcpy.da.SearchCursor(points,(fidPoint,"Shape@")) as cursor:
    for fid,shp in cursor:nodesDict[fid]=shp.firstPoint
with arcpy.da.SearchCursor(lines,(fidLine,"Shape@")) as cursor:
    for fid,shp in cursor:
        for key, point in nodesDict.iteritems():
            dist=shp.distanceTo(point)
            if dist > maxDistance:continue
            curT.insertRow((key,fid,dist))

I hope there are enough comments in script to understand it’s logic...

It took 1 min 8 seconds to process 1000 points and 1200 lines on my rather solid machine. I summarised nearLines table, points below coloured by count of lines within 100 m:

Learn Python and you’ll be able to work around licensing limitations, e.g. this and sometimes absence of extension.

[GIS] Cutting line feature by the nearest point in QGIS

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)