3

I have a Polyline and would like to add new vertex before start point or end point or at any part of it by using python in QGIS.

How can I do that?

I have used approach of following link:

How to add vertex to an existing polyline programatically?

However, when I did its procedure I got a polyline like this:

enter image description here

But I don't want the line to be snapped; just only to insert vertex.

2

This is a difficult question. When we want to introduce a new vertex (as point) inside a polyline we need to be very careful to avoid changing adequate order of each vertex in final polyline. You got above result because mentioned order was lost.

Following code was adapted of another situation where I wanted to split a polyline in several parts (equals or not) by points (it was practically not written in PyQGIS but it can be adapted to it). Points are inside code but is very easy to get its coordinates by snapping them to polyline. Algorithm to find adequate order is based in azimuths and it looks as follow:

import fiona
from shapely.geometry import shape, mapping, Point, LineString
import math
from fiona.crs import from_epsg

path = '/home/zeito/pyqgis_data/new_line.shp'

data = fiona.open(path)

crs = data.crs

line = data

geom = [ feat["geometry"] for feat in line ]

points = [ feat["geometry"]["coordinates"] for feat in line ]

line = shape(geom[0])

coords = points[0]

length = line.length

new_points = [ (400488.72071200114, 4447409.7246439075), 
               (401049.78326624096, 4447112.397201004),
               (401301.203395, 4447127.5991) ]

azimuths = []

for i in range(len(coords)-1):
    point1 = coords[i]
    point2 = coords[i+1]
    azimuths.append(math.atan2(point2[0] - point1[0], point2[1] - point1[1]))

idxs = []

k = 0

for i in range(len(coords)-1):
    for point in new_points:
        new_az = math.atan2(point[0] - coords[i][0], point[1] - coords[i][1])
        for j, item in enumerate(azimuths):
            if math.fabs(item - new_az) < 1e-6:
                idxs.append([j, new_points[k]])
                k +=1

values = []

for item in idxs:
    if item[0] not in values:
        values.append(item[0])

list = [ 1 for i in range(len(values)) ]

for i, item in enumerate(idxs):
    try:
        if idxs[i][0]== idxs[i+1][0] and values[0] == 0:
            list[idxs[i][0]] += 1
        if idxs[i][0]== idxs[i+1][0] and values[0] != 0:
            list[idxs[i][0]-1] += 1

    except IndexError:
        pass

new_idxs = [ [] for i in range(len(values)) ] 

k = 0

for i, item in enumerate(list):
    new_idxs[i].append(idxs[k][0])

    for j in range(item):
        new_idxs[i].append(idxs[k][1])
        k +=1

new_coords = []

for i in range(len(coords)-1):
    new_coords.append([coords[i],coords[i+1]])

complete_points = []

for i, element in enumerate(new_coords):
    complete_points.append(new_coords[i][0])
    for item in new_idxs:
        if item[0] == i:
            tmp = item
            for j, element in enumerate(tmp):
                if j != 0:
                    complete_points.append(element)

complete_points.append(new_coords[-1][1])

line = LineString(complete_points) 

# creation of the resulting shapefile
schema1 = {'geometry': 'Point','properties': {'id': 'int'},}
with fiona.open('/home/zeito/pyqgis_data/result.shp', 'w', 'ESRI Shapefile', schema1, crs)  as output:

    for i, point in enumerate(new_points):
        output.write({'geometry':mapping(Point(point[0], point[1])),'properties': {'id':i}})

with fiona.open('/home/zeito/pyqgis_data/result3.shp', 'w', 'ESRI Shapefile', schema1, crs)  as output:

    for i, point in enumerate(complete_points):
        output.write({'geometry':mapping(Point(point[0], point[1])),'properties': {'id':i}})


#creation of the resulting shapefile
schema2 = {'geometry': 'LineString','properties': {'id': 'int'},}

with fiona.open('/home/zeito/pyqgis_data/result2.shp', 'w', 'ESRI Shapefile', schema2, crs)  as output:

    output.write({'geometry':mapping(line),'properties': {'id':0}})


result2 = QgsVectorLayer('/home/zeito/pyqgis_data/result2.shp',
                         'line_with_added_vertexes',
                         'ogr')

QgsMapLayerRegistry.instance().addMapLayer(result2)

At following image it can be observed points to add inside polyline (in edition mode for corroborating that these points are not belong to it).

enter image description here

After running it at Python Console of QGIS I got:

enter image description here

where it can be observed, in edition mode, that points were adequately added to polyline (visualized as a new memory layer).

2

QgsGeometry has 2 methods to find the closest vertex and the closest segment to a given point. A point may be before or after an existing vertex, having the same distance to that vertex. But in that case the nearest segments are different.

In the following approach I use a combination of both methods to decide, if a new point has to be added as new first or new last vertex, or insert on the nearest segment. There is a parameter tolerance, with which to specify if the point is considered to be still on the first or last segment, or out of the polyline.

According documentation inserting a vertex before a vertex whose index is greater than the last actual vertex on the requested ring and item, it is assumed that the vertex is to be appended instead of inserted. Using Python this seems not to be the case. Therefore the last vertex is moved to the location of the point to be appended, and a new vertex is inserted at the position of the former last vertex.

To test the function vertex_add() I have included a demo which creates a memory layer and an initial polyline and that adds several vertices to this line. Simply run this script from the editor window.

def demo_line():
    lyr = QgsVectorLayer('LineString?crs=EPSG:4326', 'Line', 'memory')
    QgsMapLayerRegistry.instance().addMapLayer(lyr)
    prov = lyr.dataProvider()
    feats = []
    feat = QgsFeature()
    feat.setGeometry(QgsGeometry().fromPolyline([QgsPoint(-1, 1), QgsPoint(0, 0), QgsPoint(2, 0)]))
    feats.append(feat)
    prov.addFeatures(feats)
    lyr.triggerRepaint() 

def vertex_add(geom, x, y, tol=0.01):
    p1, at, b1, after, d1 = geom.closestVertex(QgsPoint(x, y))
    dist, p2, to = geom.closestSegmentWithContext(QgsPoint(x, y))
    if at == 0:
        if dist < tol:
            # insert into first segment
            geom.insertVertex(x, y, after)
        else:
            # insert before first vertex
            geom.insertVertex(x, y, 0)
    elif after == -1:
        if dist < tol:
            # insert after last vertex
            geom.insertVertex(x, y, at)
        else:
            # insert into last segment
            last = geom.vertexAt(at)
            geom.moveVertex(x, y, at)
            geom.insertVertex(last.x(), last.y(), at)
    else:
        # insert into any other segment
        geom.insertVertex(x, y, to)
    return geom

demo_line()

lyr = iface.activeLayer()
prov = lyr.dataProvider()
lyr.selectAll()
feat = lyr.selectedFeatures()[0]
geom = feat.geometry()

x = -1.1
y = 1.1
v = [(-0.9, 0.9), (-1.1, 1.1), (2, -1), (2.5, -1), (2.4, -1.1), (2.6, -1.1), (2.1, -1), (2, -0.9)]
for p in v:
    geom = vertex_add(geom, p[0], p[1], tol=0.015)

prov.changeGeometryValues({feat.id(): geom})
lyr.triggerRepaint() 

The result should look like the green line on the right. Tested with QGIS 2.18.4 on a Windows machine.

enter image description here

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.