TransWikia.com

Intersecting line at junction in PostgreSQL

Geographic Information Systems Asked by Nihcas on September 30, 2021

I have roads of a country having 899371 features, which were in shapefile and I had imported to PostGIS. I am trying to intersect line at every junction. I have used following commands:

1st command:

CREATE TABLE line_intersection AS
SELECT DISTINCT (ST_DUMP(ST_INTERSECTION(a.geom, b.geom))).geom AS ix
from roads_with_type a
INNER JOIN roads_with_type b
ON ST_INTERSECTS(a.geom, b.geom)
WHERE geometrytype(st_intersection(a.geom, b.geom)) = 'POINT';

2nd command:

CREATE INDEX ON line_intersection USING gist(ix);

3rd command:

CREATE TABLE line_union AS
SELECT ST_UNION(geom) as geom
FROM roads_with_type;

4th command:

CREATE INDEX on line_union USING gist(geom);

All above commands run successfully. How ever the last command I run, it takes more then 12 hours and still running. My last command is :

CREATE TABLE line_segments AS 
SELECT(ST_DUMP(ST_SPLIT(a.geom,b.ix))).geom
FROM line_union a
INNER JOIN  line_intersection b
ON ST_INTERSECTS(a.geom, b.ix);

Is there any faster way to do this? I am using PostgreSQL 12 and PostGIS 3.0.

2 Answers

Here's a solution that makes use of the fact that in PostGIS running ST_Difference on two intersecting LineStrings (or MultiLineStrings) has the effect of noding the first line. So it's possible to loop over every line in the table, and difference it with the collection of lines that intersect it, to effectively node it. The noded line is computed as a MultiLineString. This could be ST_Dumped to extract the individual lines with their original attributes.

ST_Difference of course also removes common line sections, but these can be "added back" by unioning with the intersection of the lines.

This uses a LATERAL JOIN to loop over the lines. A spatial index should be used to speed up the query for intersecting lines.

WITH data(id, geom) AS (VALUES
   ( 1, 'LINESTRING (10 10, 90 90)'::geometry )
  ,( 2, 'LINESTRING  (10 30, 30 10)'::geometry )
  ,( 3, 'LINESTRING  (60 30, 10 80)'::geometry )
  ,( 4, 'LINESTRING  (20 60, 80 60)'::geometry )
  ,( 5, 'LINESTRING  (80 60, 80 90)'::geometry )
  ,( 6, 'LINESTRING  (80 40, 80 70)'::geometry )
),
noded AS (SELECT id, ST_Union(
          ST_Difference( d1.geom, noder.geom ),
          ST_Intersection( d1.geom, noder.geom )  ) AS geom
  FROM data AS d1
  JOIN LATERAL
    (SELECT ST_Collect(d2.geom) AS geom
      FROM data AS d2
      WHERE d1.id != d2.id AND ST_Intersects(d1.geom, d2.geom)
    ) AS noder ON true
)
SELECT * FROM noded;

For the example arrangement of lines, the output is:

 id |                                         noded                                          
----+----------------------------------------------------------------------------------------
  1 | MULTILINESTRING((10 10,20 20),(20 20,45 45),(45 45,60 60),(60 60,80 80),(80 80,90 90))
  2 | MULTILINESTRING((10 30,20 20),(20 20,30 10))
  3 | MULTILINESTRING((60 30,45 45),(45 45,30 60),(30 60,10 80))
  4 | MULTILINESTRING((20 60,30 60),(30 60,60 60),(60 60,80 60))
  5 | MULTILINESTRING((80 70,80 80),(80 80,80 90),(80 60,80 70))
  6 | MULTILINESTRING((80 40,80 60),(80 60,80 70))

It would be nice to hear if this works and how long it takes.

Answered by dr_jts on September 30, 2021

Another possible solution using ST_Split is given in this PostGIS User list thread.

However, it relies on the lines having no coincident linework (since in that case ST_Split errors out).

Answered by dr_jts on September 30, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP