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

postgissimplifysnapping

I have a PostGIS table of singlepart line segments that I've run the ST_ChaikinSmoothing function on. Smoothing all of the line segments worked great, but these simplified geometries created small gaps between the other line segments in the table. (It's a network of flowlines) I need to snap the endpoints of each linestring to the closest segment of another line in the table.

I was directed to this answer about snapping lines, but I think it is snapping features in one layer (line) to features in another layer (input). I always struggle to comprehend PostGIS examples that use abstract text geometries in the queries as opposed to referencing tables of existing features.

My query below is returning an empty table with point geometry. I can successfully run each of these three functions individually to locate the points on the lines, but I can't figure out what i'm missing to snap the lines and create line geometries.

"geometry" is the name of my geometry column in the smoothed_lines table

create table snapped_lines as (
SELECT ogc_fid, name, (ST_LineInterpolatePoint(geometry,ST_LineLocatePoint(geometry,ST_EndPoint(geometry)))) as geometry from smoothed_lines);

snapped lines example

UPDATE: Sample dataset: https://www.dropbox.com/s/wz2265symhzsenl/smoothed_lines.pgsql?dl=0

The answer from @gabriel-de-luca below does snap the smoothed lines to the original lines based on matching id's, but ST_Snap moves all vertices in a linestring, not just the end nodes and wrecks the smoothness of the input features.

example after running this ST_Snap query:
enter image description here

Best Answer

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)
);
Related Question