[GIS] Snapping line segments together to remove gaps between them

arcgis-desktoplinesnapping

I have a series of line segments that were created based on the location of an existing line. There are small gaps of around 0.1 feet between the different line segments. I need to somehow snap the 2 ends of the line segments together to remove these small gaps. Manually snapping them together is not a option as there are literally 100s of these gaps in my dataset. If arcgis has a tool that will do this that is great but I'm willing to try python or VB whatever will get the job done. I have included a screencap that illustrates my problem. The green line is the original and the black segments are what I'm trying to stretch. I tried geometric networking but it is moving the whole line segment around. I just need to stretch the line. Line Snapping

Best Answer

1) Create a topology and this will fix them with the tolerance setting if you have the correct license (http://webhelp.esri.com/arcgisserver/9.3/java/index.htm#geodatabases/topology_in_arcgis.htm).

2) ET Geowizards provides a nifty little clean polyline tool that works.

3) Maybe you can just do feature class to feature class if they are all tiny and just set the environment tolerance to 1m or similar.

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)

[GIS] PostGIS snap line segment endpoint to closest other line segment

ST_LineLocatePoint( geometry, ST_EndPoint( geometry)) don't have much sense, because is searching for the closest point of a line to the endpoint of the same line, and will always return 1.
And ST_LineInterpolatePoint( geometry, 1) will always return the endpoint.

Since you want to snap geometry vertices, use ST_Snap instead. Just find the distance threshold to search a vertex to snap.

In this example, I am snapping the geometries of the lines table to the vertices of the geometries in the smoothed_lines table, that are at most to 0.5 meters away.

Both tables have an id column, and a geometry column with the geometry( LINESTRING, 54004) type. I am casting the ST_Snap() output to the same type.

The inner join is returning only the candidate pairs, and from that pairs ST_Snap() is returning all lines geometries, modified or not.

DROP TABLE IF EXISTS snapped_lines;

CREATE TABLE IF NOT EXISTS snapped_lines AS ( 
  SELECT 
    l.id AS id_l,
    s.id AS id_s,
    ST_Snap(
      l.geometry,
      s.geometry,
      0.5
    )::geometry( LineString, 54004) AS geometry
  FROM
    lines AS l
  INNER JOIN
    smoothed_lines AS s
  ON
    ST_DWithin( l.geometry, s.geometry, 0.5)
);

UPDATE:
If you need to snap lines within the same layer, you can perform the inner join to quickly exclude the pairs of lines that are far from each other, snap between each pair of selected lines, but filter the results for those lines that are more than zero apart from their pair (to avoid the lines that are generated by snapping to themselves or their adjacents).

One criterion to choose which lines must be snapped to which other, is the distance (<= 0.5) between them, adjust the value of 0.5 (reference system of the geometry column units) to the one that fits your needs.

When the lines from and to snapping, belong to the same layer, the same criterion that is applied to snapping line 1 to line 2, is applied to snapping line 2 to line 1. So let's add one more criterion than be deterministic to decide which line should snap to the other, for example, that the shortest line must snaps to the longest one.

DROP TABLE IF EXISTS snapped_lines;

CREATE TABLE IF NOT EXISTS snapped_lines AS ( 
  SELECT 
    l.id AS id_l,
    s.id AS id_s,
    ST_Snap(
      l.geometry,
      s.geometry,
      0.5
    )::geometry( LineString, 54004) AS geometry
  FROM
    smoothed_lines AS l
  INNER JOIN
    smoothed_lines AS s
  ON
    ST_DWithin( l.geometry, s.geometry, 0.5)
  WHERE
    ST_Distance( l.geometry, s.geometry) > 0
    AND
    ST_Length( l.geometry) < ST_Length( s.geometry)
);

Best Answer

Related Solutions

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

Preamble

Workflow explanation

Solution with edits applied to the original layer

Solution with edits applied to a new layer

[GIS] PostGIS snap line segment endpoint to closest other line segment

Related Question