QGIS – How to Snap Multiple Nodes of a Polygon to Nodes of a Line

digitizingnodespolygonqgissnapping

I am working with vector polygons that have derived from hand drawn areas from when mapping was carried out with tracing paper over paper maps. Therefore many of the polygons don't exactly line up to the boundaries on the ground.

Another dataset of vector lines I have is accurate and so I have been manually snapping nodes of the polygons to the nodes of the lines. But this is a lengthy process.

Here is a screenshot…

The green line is the accurate dataset. You can see how the polygon boundaries (orange) are roughly correct. Bottom left I have started to edit the nodes, with the snapping options set up, and this works to manually move the nodes of the polygons to snap to the nodes of the line. (Yes, I could have enabled topographical editing to move both the polygon boundaries together, and I do do that).

Is there a method for snapping multiple nodes of the polygon to the nodes of the lines, or is there perhaps a better method of doing this? Using Python perhaps? I am looking for something akin to the 'Enable tracing' function (which is great) but for moving polygon boundaries.

I am aware that this may in places be a duplicate as there are many similar questions, but I haven't quite found anything that answers my particular situation? In particular I have looked at…

QGIS how to align existing polygon borders

How to snap polygons to lines?

How to snap lines to nearest point?

https://gis.stackexchange.com/search?q=how+to+snap+many+polygons+to+lines

Best Answer

Preamble

I propose a solution using PyQGIS. This is a complex issue and I'm sure that there are (many) specific cases that I didn't consider in my approach (so, take this answer as a starting point rather than a complete solution).

I didn't deeply test the code since I was not sure about the specific issue (see my unanswered comment where I requested some explanations), so please do some tests for validating it: if you are going to use the first attached code, create a backup copy of your layers because the edits are immediately committed to the geometry; instead, if you want to do some test without affecting the original layer, try to use the last attached code (both codes are the same, the unique difference is that the first one commits the edits on the layer itself, while the second one commits them to a new memory layer).

Finally, let's assume to start from this situation, where the red line is the guideline and the two polygons are two separate layers (each one having one feature):

Workflow explanation

The solution requires a line layer as a guideline and a set of polygon layers that need to be snapped. My approach iterates over the nodes of every polygon feature and tries to make a collage with the line geometry by deleting the unnecessary vertices and adding new ones.

This approach is quite simple to understand, but there are some things to consider:

while the polygon vertices are automatically ordered clockwise when digitizing (i.e. there is only one possible direction), a line layer can be digitized in two different directions according to which side you start from (the final direction is obviously unique, but you don't know if it follows the same order of the polygon vertices);
as a consequence of the previous point, it is very easy creating features that are geometrically wrong and it basically depends on how you are able to consider all the different situations that may occur (this is the part where my solution needs to be improved).

Having said that, my workflow could be summarized in this way:

read the current line feature geometry as a sequence of coordinates (they are returned in a list);
read the current polygon feature geometry as a sequence of coordinates (they are returned in a list);
from the current polygon feature, find the closest vertex to both line start and line end;
edit the list from point 2 by assigning the coordinates within the list from point 1 in the proper position.

As you can understand, the main issues are in the last step because it didn't seem too easy enumerating all the possible cases (it would be easier if you provide a sample dataset).

Solution with edits applied to the original layer

You may run the following code as a new script from the Processing Toolbox:

##input_guideline=vector line
##input_polygons=multiple vector

from qgis.core import *
from qgis.utils import iface

polygons = input_polygons.split(';')
line = QgsVectorLayer(input_guideline, input_guideline, 'ogr')

for linefeat in line.getFeatures():
    line_geom = linefeat.geometry()
    geom_to_add = line_geom.asPolyline()
    (first_line_point, end_line_point) = (line_geom.interpolate(0), line_geom.interpolate(-1))

    for polygon in polygons:
        poly = QgsVectorLayer(polygon, polygon, 'ogr')
        for feat in poly.getFeatures():
            attrs = feat.attributes()
            geom = feat.geometry()
            coords = geom.asPolygon()[0]

            first_vertex=geom.closestVertex(first_line_point.asPoint())[0]
            last_vertex=geom.closestVertex(end_line_point.asPoint())[0]
            first_indices = [i for i, x in enumerate(coords) if x == first_vertex]
            first_index = first_indices[0]
            second_index = coords.index(last_vertex)
            if first_index < second_index:
                coords[first_index:second_index] = geom_to_add
                new_geom = QgsGeometry.fromPolygon([coords])
                if not new_geom.isGeosValid():
                    first_index = first_indices[-1]
                    if first_index < second_index:
                        coords = geom.asPolygon()[0]
                        coords[second_index:first_index] = (geom_to_add)
                        new_geom = QgsGeometry.fromPolygon([coords])
                    else:
                        coords = geom.asPolygon()[0]
                        coords[second_index:first_index] = reversed(geom_to_add)
                        new_geom = QgsGeometry.fromPolygon([coords])  
            else:
                coords[first_index:second_index] = reversed(geom_to_add)
                new_geom = QgsGeometry.fromPolygon([coords])
                if not new_geom.isGeosValid():
                    first_index = first_indices[-1]
                    coords = geom.asPolygon()[0]
                    coords[second_index:first_index] = reversed(geom_to_add)
                    new_geom = QgsGeometry.fromPolygon([coords])

            poly.dataProvider().changeGeometryValues({feat.id(): new_geom})
iface.mapCanvas().refreshAllLayers()

and you will get this result:

If everything works, the edited polygons will perfectly match the guideline and have both the edges and the vertices perfectly aligned to it.

Solution with edits applied to a new layer

This solution differs from the previous one because saves the edits to a new memory layer instead of committing the changes to the original layer (it could be useful for testing):

##input_guideline=vector line
##input_polygons=multiple vector

from qgis.core import *

polygons = input_polygons.split(';')
line = QgsVectorLayer(input_guideline, input_guideline, 'ogr')

for linefeat in line.getFeatures():
    line_geom = linefeat.geometry()
    geom_to_add = line_geom.asPolyline()
    (first_line_point, end_line_point) = (line_geom.interpolate(0), line_geom.interpolate(-1))

    for polygon in polygons:
        poly = QgsVectorLayer(polygon, polygon, 'ogr')

        # Create the output layer
        crs = poly.crs().toWkt()
        outLayer = QgsVectorLayer('Polygon?crs='+ crs, 'snapped' , 'memory')
        prov = outLayer.dataProvider()
        fields = poly.pendingFields() # Fields from the input layer
        prov.addAttributes(fields) # Add input layer fields to the outLayer
        outLayer.updateFields()


        for feat in poly.getFeatures():
            attrs = feat.attributes()
            geom = feat.geometry()
            coords = geom.asPolygon()[0]

            first_vertex=geom.closestVertex(first_line_point.asPoint())[0]
            last_vertex=geom.closestVertex(end_line_point.asPoint())[0]
            first_indices = [i for i, x in enumerate(coords) if x == first_vertex]
            first_index = first_indices[0]
            second_index = coords.index(last_vertex)
            if first_index < second_index:
                coords[first_index:second_index] = geom_to_add
                new_geom = QgsGeometry.fromPolygon([coords])
                if not new_geom.isGeosValid():
                    first_index = first_indices[-1]
                    if first_index < second_index:
                        coords = geom.asPolygon()[0]
                        coords[second_index:first_index] = (geom_to_add)
                        new_geom = QgsGeometry.fromPolygon([coords])
                    else:
                        coords = geom.asPolygon()[0]
                        coords[second_index:first_index] = reversed(geom_to_add)
                        new_geom = QgsGeometry.fromPolygon([coords])  
            else:
                coords[first_index:second_index] = reversed(geom_to_add)
                new_geom = QgsGeometry.fromPolygon([coords])
                if not new_geom.isGeosValid():
                    first_index = first_indices[-1]
                    coords = geom.asPolygon()[0]
                    coords[second_index:first_index] = reversed(geom_to_add)
                    new_geom = QgsGeometry.fromPolygon([coords])

            outGeom = QgsFeature()
            outGeom.setAttributes(attrs)
            outGeom.setGeometry(new_geom)
            prov.addFeatures([outGeom])

        # Add the layer to the Layers Panel
        QgsMapLayerRegistry.instance().addMapLayer(outLayer)

Related Solutions

[GIS] Snapping start- and end nodes of lines to other lines in PostGIS

I've managed to solve this, without using the mentioned GRASS tools or topological functions.

Basically I take all start- and endnodes, put them in a new, temporary table, put a buffer around them, union the buffer objects, and move all found nodes in each buffer to the centroid of the buffer.

When that's done I move the original begin and end points to the new location.

Easier than expected, and still fast, but I expected PostGIS to have some built-in function for this - that would be even quicker.

Edit: in the interest of giving back to the community, this is my (quite crappy) code for now.

drop table if exists nodes;
drop table if exists nodes2;
drop table if exists buffers;

-- Get Start and End nodes
select ST_StartPoint(wkb_geometry) startnode,  ST_EndPoint(wkb_geometry) endnode,    ogc_fid into nodes  from sourceTable;
-- Combine all nodes into one table for easier queries
select startnode node, ogc_fid into nodes2 from nodes;
insert into nodes2 select endnode node, ogc_fid from nodes;

-- Some indexes to speed everything up
CREATE INDEX nodesstart_idx ON nodes USING gist  (startnode);
CREATE INDEX nodesend_idx ON nodes USING gist  (endnode);
CREATE INDEX nodes2_idx ON nodes2 USING gist  (node);
CREATE INDEX nodes_ogcfid_idx ON nodes USING btree (ogc_fid ASC NULLS LAST);

-- Create buffers, combine them, split combined objects again
select (ST_Dump(ST_Union(ST_Buffer(node, 1)))).geom geom into buffers from nodes2;
CREATE INDEX buffers_idx ON buffers USING gist  (geom);

-- Update start/end nodes table
UPDATE nodes SET startnode = ST_Centroid((select geom from buffers WHERE geom && startnode));
UPDATE nodes SET endnode = ST_Centroid((select geom from buffers WHERE geom && endnode));
-- Update original points
update sourceTable set wkb_geometry = ST_SetPoint(
ST_SetPoint(wkb_geometry, 0, (select startnode from nodes where ogc_fid=sourceTable.ogc_fid)), 
ST_NumPoints(wkb_geometry) - 1, (select endnode from nodes where ogc_fid=sourceTable.ogc_fid));

DROP TABLE nodes;
DROP TABLE nodes2;
DROP TABLE buffers;

[GIS] QGIS snapping not working for me

To close this off. I'll provide my own answer as per commentator's suggestion.

To recap, I could not get nodes on adjacent polygons to snap together when I hovered one node over the other. It turns out that the tolerance was set to zero by default.

In Settings - Snapping options, set the tolerance to 6 pixels (see screenshot). A magenta 'plus' sign appears over the node when it gets in range to snap to another node. After experimenting a while now, I think this is a good value to start with on any project.

Ensure the snapping tolerance is set to a sensible value