2

I have some LineStrings that overlap in some segments, but not necessarily with coinciding vertices, like shown in the example below. With PostGIS I would like to extract all "segments" and count how many LineStrings are overlapping here (1 would mean that nothing overlaps).

As an example, consider the following situation: Two LineStrings (LINESTRING(0 0, 100 40) and LINESTRING(50 20, 150 60)) and a point P (POINT(75 25)).

Easy to see, that the LineStrings overlap between B(50,20) and C(100,40) as shown here:

enter image description here

My approach is:

  • Extract all (distinct, in the case of duplicate vertices) vertices (A,B,C,D in a first step, see below) from the LineStrings and make a MULTIPOINT of it
  • ST_Split the LineStrings and count the geometries

I achieve this with the following query (sample data included):

with
    sample(geom, id) as (values
        (st_makeline(array[st_makepoint(0, 0), st_makepoint(100, 40)]), 1),
        (st_makeline(array[st_makepoint(50, 20), st_makepoint(150, 60)]), 2)),

    splitter(geom) as (select st_union(distinct x.geom) as splitter from (select (st_dumppoints(geom)).*  from sample) as x)

select count(y.geom), st_asText(y.geom) as segment from (select (st_dump(st_split(s.geom, splitter.geom))).* from splitter, sample as s) as y group by y.geom

returning the expected result (2 occurences of B-C):

count | segment
----------------------------
1     | LINESTRING(0 0,50 20)
2     | LINESTRING(50 20,100 40)
1     | LINESTRING(100 40,150 60)

To make things a little more interesting, I add a point P' to LineString B-D (cp. above, LINESTRING(50 20, 150 60)), which is the closest point on B-D to P.

Running the same query again (with the modifications to the sample data mentioned above)...

with 
    sample(geom, id) as (values
        (st_makeline(array[st_makepoint(0, 0), st_makepoint(100, 40)]), 1),
        (st_addpoint(st_makeline(array[st_makepoint(50, 20), st_makepoint(150, 60)]),
                     st_closestpoint(st_makeline(array[st_makepoint(50, 20), st_makepoint(150, 60)]), st_makepoint(75, 15)),1), 2)),

    splitter(geom) as (select st_union(distinct x.geom) as geom from (select (st_dumppoints(geom)).*  from sample) as x)

select count(y.geom), st_astext(y.geom) as segment from (select (st_dump(st_split(s.geom, splitter.geom))).* from splitter, sample as s) as y group by y.geom

... I expect (not knowing the exact coordinates of P) six segments, with two of them twice:

count | segment
----------------------------
1     | LINESTRING(0 0,50 20)
2     | LINESTRING(50 20, P'x P'y)
2     | LINESTRING(P'x P'y, 100 40)
1     | LINESTRING(100 40,150 60)

... but I get

count | segment
-------------------------------
1     | LINESTRING(0 0,50 20)
1     | LINESTRING(50 20,69.8275862068966 27.9310344827586)
1     | LINESTRING(50 20,100 40)
1     | LINESTRING(69.8275862068966 27.9310344827586,150 60)

So apparently, C does not split B-P'-D (but P' does) and P' does not split A-C (but B does).

Selection of the vertices works fine:

-- create sample data
select distinct st_astext(x.geom) as vertices from (select (st_dumppoints(geom)).*  from sample) as x

vertices
-------------------------
POINT(0 0)
POINT(100 40)
POINT(150 60)
POINT(50 20)
POINT(69.8275862068966 27.9310344827586)

So ST_Split works only when points lie exactly on the LineStrings, but I would expect that P' matches this condition for both LineStrings (or none).

Does adding P' modify the geometry in a way that ST_Split does not work as expected or what might I am doing wrong here? Or is the whole approach likely to be unsuitable? Or is ST_ClosestPoint unsuitable to create P'?

  • 1
    do you get different results if you guarantee both lines are equal by e.g. snapping them to a grid? after all, P' is created on/for one line only and 13 decimal digits leave a lot of room for inequality... – ThingumaBob Feb 6 '18 at 15:59
3

Your frustration with cutting lines with points is not unmarked. Check out this feature request from 2013 asking to add tolerance parameter to ST_Split for splitting lines by points.

The solution is to snap the lines to the points. Note that snapping the points to the lines does not work.

Your example with ST_Snap() added:

with
    sample(geom, id) as (values
        (st_makeline(array[st_makepoint(0, 0), st_makepoint(100, 40)]), 1),
        (st_addpoint(st_makeline(array[st_makepoint(50, 20), st_makepoint(150, 60)]),
                     st_closestpoint(st_makeline(array[st_makepoint(50, 20), st_makepoint(150, 60)]), st_makepoint(75, 15)),1), 2)),

    splitter(geom) as (select st_union(distinct x.geom) as geom from (select (st_dumppoints(geom)).*  from sample) as x)

select count(y.geom), st_astext(y.geom) as segment 
from (
    select (st_dump(st_split(ST_Snap(s.geom, splitter.geom, 0.00000001), splitter.geom))).* 
    from splitter, sample as s
) as y group by y.geom

Output:

count | segment
-------------------------------
1     | LINESTRING(0 0,50 20)
2     | LINESTRING(50 20,69.8275862068966 27.9310344827586)
2     | LINESTRING(69.8275862068966 27.9310344827586,100 40)
1     | LINESTRING(100 40,150 60)

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.