7

I have a few linestrings that together form a long line. The start- and endpoints of the linestrings are equal with their direct neighbour, but when tested the ST_Distance between any 2 lines is like 1.0x10^-11.

This seems to be the known issue with encoding floats. They are the same but not the same.

However when I update 2 neighbouring lines with their existing Astext coordinates they are encoding by the same algorithm and their distance is 0. But then they have a small distance with their outside neighbours (that are not updated this way).

What I want to do is snap a line left and right with its neighbours so the ST_Distance with both neighbours is 0. I think ST_Snap only moves a geometry to another geometry. That is not what I want. I only want the start- and endpoint of a line to snap to the neighbour. The rest of the line should stay where it is.

How do I do that?

With the given answer I fixed both ends of the conflicting line like this:

update network set the_geom = ST_RemovePoint(the_geom, (select st_npoints(the_geom)-1 from network where gid=6909)) where gid=6909;
update network set the_geom = ST_AddPoint(the_geom, (select st_endpoint(the_geom) from network where gid=6123), (select st_npoints(the_geom) from network where gid=6909)) where gid=6909;
update network set the_geom = ST_RemovePoint(the_geom, 0) where gid=6055;
update network set the_geom = ST_AddPoint(the_geom, (select st_startpoint(the_geom) from network where gid=6909), 0) where gid=6055;
  • you could try extract start and end points, snap them and replace old start and end points with data you get from st_snap – simplexio Mar 8 '13 at 8:38
  • 1
    There is ST_AddPoint(line, point, index), ST_RemovePoint(line, index) but you can use only ST_ReplacePoint (line,index,point) then you need ST_PointN(line, index) and ST_NumPoints(line). Add St_Snap twice and replace you old end points with new ones – simplexio Mar 8 '13 at 8:47
  • simplexia - I used St_RemovePoint and ST_AddPoint as a way to replace a vertex (as I couldn't find the ST_ReplacePoint you mentioned). It worked. If you upgrade your comment to an answer I will accept it as the answer I was looking for. – mrg Mar 8 '13 at 14:24
  • Thanks you. Could you add example SQL to asnwer ? – simplexio Mar 11 '13 at 7:33
  • A weird thing is that ST_RemovePoint and ST_AddPoint are 0 indexed but ST_PointN is 1 indexed. – mrg Mar 11 '13 at 7:55
4

You can use ST_RemovePoint(line, index) and ST_AddPoint(line, point, index). First, get first and last point with ST_pointN(line, index) and ST_NumPoints(line), then use ST_Snap to those points. Then use returned points to replace the old ones.

These functions should be in Postgis < 2.0 ST_ReplacePoint hould be in Postgis 2.0

0

I realize this is an ancient question, but I just landed here and want to share another option to clean up lines with tiny gaps between them. I had the same issue as the original question - I had multilinestrings made up of linestrings that had minuscule gaps between them (less than a cm). This was preventing me from merging these multilinestrings into linestrings for linear referencing. The solution below removes the gaps automatically.

My answer is taken from this post Cleaning Road Networks with PostGIS Topology

First, add the topology extension to your postgis database. You can either create a new database or add the extension to an existing database.

CREATE DATABASE "roadsdb";
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;
SET search_path = topology,public;

Next, create topology for your lines layer - roads_topo is the name of the new topology schema that will be created, 2154 is the epsg code for the spatial reference:

SELECT topology.CreateTopology('roads_topo', 2154);
SELECT topology.AddTopoGeometryColumn('roads_topo', 'public', 'roads', 'topo_geom', 'LINESTRING');

Convert the lines table to nodes and edges, setting tolerance to 1m. Adjust the 1.0 value below to the tolerance that is acceptable to you (the largest gaps you need to close). If you're not using a projected coordinate system you will need to convert this value to degrees. In WGS84, 1 degree is around 100km:

UPDATE roads SET topo_geom = topology.toTopoGeom(wkb_geometry, 'roads_topo', 1, 1.0);

Only if you encounter errors and toTopoGeom fails, you may need to insert records one at a time and skip errors.

DO $$DECLARE r record;
BEGIN
  FOR r IN SELECT * FROM roads LOOP
    BEGIN
      UPDATE roads SET topo_geom = topology.toTopoGeom(wkb_geometry, 'roads_topo', 1, 1.0)
      WHERE ogc_fid = r.ogc_fid;
    EXCEPTION
      WHEN OTHERS THEN
        RAISE WARNING 'Loading of record % failed: %', r.ogc_fid, SQLERRM;
    END;
  END LOOP;
END$$;

You can then export a new 'clean' layer using the new topo_geom column:

CREATE TABLE roads_clean AS (
    SELECT ogc_fid, topo_geom::geometry
    FROM roads
);
  • Thx for sharing! – mrg Oct 21 at 19:24

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.