2

I try to remove the overlap of line geometries in one table based on a class attribute. My table consists of two columns containing line geometries and a number "class" from 1 to 3.

CREATE TABLE lines(class INT, geom geometry );

INSERT INTO  lines VALUES
    (1, 'SRID=4326;LINESTRING(1 1, 1 4)'),
    (2, 'SRID=4326;LINESTRING(1 3, 1 6)'),
    (3, 'SRID=4326;LINESTRING(1 5, 1 9)');

Now I'd like to extract the part of geometry 2 wich is not intersecting with geometry 1 and the part of geometry 3 wich is not intersecting with geometry 2. I did this already in Python:

class_1 = shp[(shp.Class == 1)]
class_2 = shp[(shp.Class == 2)]
class_3 = shp[(shp.Class == 3)]


part_1 = gpd.overlay(class_2,class_1, how='difference') 
part_2 = gpd.overlay(class_3,class_2, how='difference')
part_3 = class_1

but running following code, the query comes to no end:

SELECT
CASE WHEN lines.class = 2 THEN ST_Difference(lines.geom, (SELECT ST_Union(lines.geom) FROM lines WHERE class = 1 ))
     WHEN lines.class = 3 THEN ST_Difference(lines.geom, (SELECT ST_Union(lines.geom) FROM lines WHERE class = 2 ))
ELSE lines.geom
END
FROM lines;

How can I write the code most clever?

0

4 Answers 4

2

It works on a small subset. The issue is that for every row in the table, you are building the union of all lines, which takes forever on a large dataset.

The query needs to build the geometry union only of line segments that intersect the "current" one.

To simplify a bit, since you are calling the same table several times, use a different table name alias


WITH lines(class, geom ) AS (VALUES
    (1, 'SRID=4326;LINESTRING(1 1, 1 4)'::geometry),
    (2, 'SRID=4326;LINESTRING(1 3, 1 6)'),
    (3, 'SRID=4326;LINESTRING(1 5, 1 9)'))
SELECT
CASE 
    WHEN l1.class = 2 THEN 
        ST_Difference(l1.geom, 
        (SELECT ST_Union(l2.geom) FROM lines l2 
         WHERE l2.class = 1 AND st_intersects( l1.geom,l2.geom)))
    WHEN l1.class = 3 THEN 
    ST_Difference(l1.geom, 
    (SELECT ST_Union(l2.geom) FROM lines l2 
     WHERE l2.class = 2 AND st_intersects( l1.geom,l2.geom)))
ELSE l1.geom
END
FROM lines l1;

But you don't need the case conditions, you can just use the "current" class -1 or, if there is no intersection, the current geometry

WITH lines(class, geom ) AS (VALUES
    (1, 'SRID=4326;LINESTRING(1 1, 1 4)'::geometry),
    (2, 'SRID=4326;LINESTRING(1 3, 1 6)'),
    (3, 'SRID=4326;LINESTRING(1 5, 1 9)'))
SELECT coalesce(
        ST_Difference(l1.geom, 
        (SELECT ST_Union(l2.geom) FROM lines l2 
         WHERE l2.class = l1.class-1 AND st_intersects( l1.geom,l2.geom))),l1.geom)
FROM lines l1;

This can then be rewritten using a lateral join:

WITH lines(class, geom ) AS (VALUES
    (1, 'SRID=4326;LINESTRING(1 1, 1 4)'::geometry),
    (2, 'SRID=4326;LINESTRING(1 3, 1 6)'),
    (3, 'SRID=4326;LINESTRING(1 5, 1 9)'))
SELECT coalesce(
        ST_Difference(l1.geom, 
        ljoin.geom),l1.geom)
FROM lines l1
CROSS JOIN LATERAL 
  (SELECT ST_Union(l2.geom) as geom 
   FROM lines l2 
   WHERE l2.class = l1.class-1 
   AND st_intersects( l1.geom,l2.geom)) ljoin;

1
  • Fantastic solution, especially with the lateral join. I was able to reduce file sizes by a significant amount.
    – Stonetip
    Dec 8, 2021 at 1:59
0

Since you want to get the part of geom(ID = n) that do not intersect geom(ID = n-1) you can use the LEAD() function:

SELECT LEAD(class) OVER() as class, 
       ST_Difference(LEAD(geom) OVER(), geom) as geom 
FROM lines
1
  • You need to order the records in the over() clause, else the previous row is a random one.
    – JGH
    Aug 27, 2020 at 15:51
0

So, if I understood the translation correctly and pondered your text, and you need to remove the overlap as shown in the figure below, i.e to get dotted lines, then please visit me, if not then correct me or pass by 🙂

enter image description here

Run CTE on your table:

WITH
tbla AS (SELECT DISTINCT a.class, ST_Intersection(a.geom, LEAD(a.geom) OVER (ORDER BY a.class)) geom FROM lines a, lines b WHERE ST_Intersects(a.geom,b.geom)),
tblb AS (SELECT DISTINCT a.class, ST_Intersection(a.geom, b.geom) geom FROM tbla a, lines b WHERE a.class<b.class AND ST_Intersects(a.geom,b.geom))
SELECT a.class, (ST_Dump(ST_Difference(a.geom, b.geom))).geom AS geom FROM lines AS a CROSS JOIN LATERAL (SELECT ST_Collect(geom) AS geom FROM tblb WHERE ST_Intersects(a.geom, geom)) AS b

And check your result.

"If you come with gifts, give them if you forget to bring them another time"...🙂

Original solutions...

The script is called: ST_DashedLine

1
  • This is the answer to creating dotted lines in PostGIS, and if you have a question, I'll move it to another place, but for now let it be here ... 🙂 Sep 7, 2020 at 6:44
0

The literal translation of the Python code is something like this:

WITH lines(class, geom ) AS (VALUES
    (1, 'SRID=4326;LINESTRING(1 1, 1 4)'::geometry),
    (2, 'SRID=4326;LINESTRING(1 3, 1 6)'),
    (3, 'SRID=4326;LINESTRING(1 5, 1 9)')),
lines1 AS ( SELECT ST_Collect( geom ) AS geom FROM lines WHERE class = 1),
lines2 AS ( SELECT ST_Collect( geom ) AS geom FROM lines WHERE class = 2),
lines3 AS ( SELECT ST_Collect( geom ) AS geom FROM lines WHERE class = 3)
SELECT 
  ST_Difference( lines2.geom, lines1.geom ) AS part1,
  ST_Difference( lines3.geom, lines2.geom ) AS part2,
  lines1.geom AS part3
  FROM lines1, lines2, lines3;

This should provides a similar execution profile as the Python code as well.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.