7

I have a line network and I am placing points at equal intervals along each line. These points will be used to cut the lines into segments. I have found a way to create the points along the lines:

CREATE TABLE split_pt as

WITH line AS 
    (SELECT
        id,
        geom 
    FROM line_table),
linemeasure AS
    (SELECT
        ST_AddMeasure(line.geom, 0, ST_Length(line.geom)) AS linem,
        generate_series(0, ST_Length(line.geom)::int, 160.9) AS i
    FROM line),
geometries AS (
    SELECT
        i,
        (ST_Dump(ST_GeometryN(ST_LocateAlong(linem, i), 1))).geom AS geom 
    FROM linemeasure)

SELECT
    row_number() over() as id,
    i,
    ST_SetSRID(ST_MakePoint(ST_X(geom), ST_Y(geom)), 3857) AS geom
FROM geometries;

However, I have two problems:

  1. I'm getting around 400 duplicate points
  2. About 60% of the points are not actually on the line so running the subsequent ST_Split() with these points against the lines doesn't work in most instances.

I am looking for help to eliminate duplicate points and to ensure that the points I create in the above query will return True on ST_Intersects(). I'd like to be able to incorporate these two aspects in the the above query and not in subsequent steps, if possible.

15
  • 2
    could you add a simple input geometry and the result of the intermediate results – Ian Turton Oct 8 '19 at 16:43
  • 3
    Due to the limitations of IEEE floating-point values, it might not be possible for a point to land exactly on a line. – Vince Oct 8 '19 at 17:36
  • 1
    doubling point removal: DELETE from geometries WHERE NOT id IN( SELECT min(id) as idpnt from geometries group by geom); – Cyril Mikhalchenko Oct 8 '19 at 19:06
  • 1
    It's a pity that this is actually the case, because users assume that the point should be exactly on the line, and that's a lie, it would still be a lie to understand what the concept of a line in postgis is - it's hair, it's nano hair, how thick the line is, etc., etc. ? I've seen the JGH answer, but it doesn't solve your problem... – Cyril Mikhalchenko Oct 9 '19 at 16:36
  • 2
    PostGIS interpolates between vertices, it's functional minimal unit for most geometries (e.g. a line) is a point array, and trigonometry; floating point error is the issue here. ST_QuantizeCoordinates can help, but really you should be looking at ST_LineSubstring. – geozelot Oct 9 '19 at 21:43
12

Better look at using ST_LineSubstring; no issues with coordinate precision, but somewhat tricky to set up...


Edit: I then packed it into a more sophisticated set of functions a while ago; there's a C function add-on as well, for those who like to build from source.

I packed that functionality into a (naive) function a few years ago:

CREATE OR REPLACE FUNCTION ST_LineSubstringByLength(

  geom       GEOMETRY(LINESTRING),
  length     FLOAT8,
  use_meter  BOOLEAN DEFAULT TRUE

) RETURNS SETOF geometry_dump AS

  $$
  DECLARE

    len_frac  FLOAT8;
    s_frac    FLOAT8;
    e_frac    FLOAT8;

  BEGIN

    IF ($3)
      THEN  len_frac := $2 / ST_Length(geom::GEOGRAPHY);
      ELSE  len_frac := $2 / ST_Length(geom);
    END IF;

    FOR n IN 0..CEIL(1.0 / len_frac)
      LOOP
        s_frac := len_frac * n;
        IF (s_frac >= 1.0)
          THEN
            EXIT;
        END IF;
        e_frac := len_frac * (n + 1);
        IF (e_frac > 1.0)
          THEN
            e_frac := 1.0;
        END IF;
        RETURN NEXT (ARRAY[n + 1], ST_LineSubstring($1, s_frac, e_frac));
      END LOOP;

    RETURN;

  END
  $$
  LANGUAGE plpgsql;

The function accepts:

  • a LineString GEOMETRY
  • a Length value by which the given LineString should be divided
    (until the last element which will be of the rest length)
  • a Switch to set the given length to be interpreted as meter; defaults to TRUE
    (input geometry needs to be in EPSG:4326; it really just casts to GEOGAPHY for the fraction calculation)

The function returns a SETOF geometry_dump (a RECORD just like e.g. ST_Dump) having:

  • a path value (INTEGER[]) indicating the position of the current row in the extraction starting from the beginning of the input LineString
  • a geom member (GEOMETRY)

A record looks like:

({1}, 01020000000200000000000000000000000000000000000000A65F8237663FC13F0000000000000000)
({2}, 010200000002000000A65F8237663FC13F0000000000000000A65F8237663FD13F0000000000000000)

Usage:

SELECT * FROM ST_LineSubstringByLength('LINESTRING(0 0, 5 0)'::GEOMETRY, 15000);
SELECT * FROM ST_LineSubstringByLength('LINESTRING(0 0, 5 0)'::GEOMETRY, 1, FALSE);

SELECT id, (dmp).path[1], (dmp).geom
FROM   <table_with_geom_and_id_columns>,
       LATERAL ST_LineSubstringByLength(geom, <length_in_meter?>[, TRUE|FALSE]) AS dmp
;
1
  • 1
    It's trivial btw. to alter this function to accept an integer value of segments it should divide the input LineString into. I use a customized PostGIS version on my machine where I add all those things that I can directly into the C function base; it's trivial to have these functions then to return segments for any geometry (using e.g. the boundaries)...and it's about 10 to 100 times faster. if you don't mind building from source...? – geozelot Oct 24 '19 at 13:48

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.