5

I'd like to know a way to detect overshoot lines (such as flyovers, over/underpasses) in road network through queries in PostGIS to validate topological connections.

This is a sample network:

CREATE TABLE net (id integer, geom geometry(MultiLineString, 4326));
INSERT INTO net (id, geom)
VALUES (1,ST_GeomFromText('MULTILINESTRING((136.63964778201967 36.58031552554121,136.64078637948796 36.57968013023401,136.6414207216891 36.579313967665655,136.64203092428795 36.578959650067986,136.6428838673968 36.57843301697034,136.6438347098042 36.577844992552514))', 4326)),
       (2,ST_GeomFromText('MULTILINESTRING((136.64039880046448 36.57933335255234,136.64078637948796 36.57968013023401,136.64119407544638 36.580076445271914,136.64165541506554 36.58051798901414))', 4326)),
       (3,ST_GeomFromText('MULTILINESTRING((136.64039880046448 36.57933335255234,136.64044439789086 36.578869185464725))', 4326)),
       (4,ST_GeomFromText('MULTILINESTRING((136.6409003730539 36.57873995108804,136.6414207216891 36.579313967665655,136.64222404380473 36.580108753416425))', 4326)),
       (5,ST_GeomFromText('MULTILINESTRING((136.6416232292288 36.57843409435816,136.64203092428795 36.578959650067986))', 4326)),
       (6,ST_GeomFromText('MULTILINESTRING((136.64203092428795 36.578959650067986,136.64244666738 36.57935489221467,136.64283290461503 36.57974259264671))', 4326)),
       (7,ST_GeomFromText('MULTILINESTRING((136.64250969906357 36.57802376968141,136.6428838673968 36.57843301697034,136.64326608196427 36.578854108330574))', 4326)),
       (8,ST_GeomFromText('MULTILINESTRING((136.64364829653186 36.577513284810436,136.6438347098042 36.577844992552514))', 4326)),
       (9,ST_GeomFromText('MULTILINESTRING((136.6438347098042 36.577844992552514,136.64412438862985 36.57833070559775))', 4326));

As illustrated in the following image, this network has 9 lines where line 2, 4 and 7 do not topologically connected to line 1. I need to generate points and export it as a new table at the corossing points of those links as shown in the right-hand side figure.

I tried to search for the solutions but could not find so far.

enter image description here

Actually I posted a similar question Detecting overshoot lines in road network (flyovers, over/underpasses) in QGIS but I would like to be able to manipulate it in PostGIS as well.

HSJ
  • 355
  • 1
  • 7

3 Answers3

4

pgRouting has a comprehensive set of functions to address most of the issues you are facing, with regards to your series of questions about network topologies:

I'll try to answer them all with a set of queries to create a routable network for most input lineworks.


  1. Dump your MultiLineStrings:

    A proper topology should consist of simple LineStrings - for pgRouting it is mandatory:

    CREATE TABLE <network_simple> AS (
      SELECT
        ROW_NUMBER() OVER() AS id,
        id AS old_id,
        dmp.path[1] AS old_id_seq,
        <attribute_1>,
        ...
        dmp.geom::GEOMETRY(LINESTRING, 4326) AS geom
      FROM
        <network>,
        LATERAL ST_Dump(geom) AS dmp
    );
    

    ALTER TABLE <network_simple> ADD PRIMARY KEY (id) ;

  2. Create a noded network:

    Split the resulting linework at intersections:

    SELECT
      pgr_nodeNetwork(
        '<network_simple>',
        0.000001,
        'id',
        'geom'
      )
    ; 
    
  3. Create a pgRouting topology:

    Add source/target information to your <network_simple>_noded table:

    SELECT
      pgr_createTopology(
        '<network_simple>_noded',
        0.000005,
        'geom',
        'id',
        'source',
        'target'
      )
    ;
    

Notes:

  • the resulting network topology in <network_simple>_noded maintains a reference to <network_simple> via <network_simple>_noded.old_id = <network_simple>.id
  • the gap issue is addressed by providing the tolerance parameter in pgr_createTopology; while it does not alter the geometries, it treats those within distance of the given tolerance as connected

  • while it is possible to brute force a (interconnected) line merge (i.e. simplification) and noding operation via
    SELECT
      ROW_NUMBER() OVER() AS id,
      dmp.geom
    FROM
      (
        SELECT
          ST_Node(ST_LineMerge(ST_Union(geom))) AS geom
        FROM
          <network>
      ) AS nu,
      LATERAL ST_Dump(nu.geom) AS dmp
    ;
    
    you'd loose all (unique) references to your <network> table: if you are okay with that, just execute 3. from above afterwards
  • if you were to run a graph contraction operation (pgr_contraction) on the <network_simple>_noded, you would likely solve the gaps and line merge issues, but you have to be aware that this algorithm is meant for edge/vertex compression, and the resulting topology would likely not look like your input network due to altered edges
geozelot
  • 30,050
  • 4
  • 32
  • 56
  • It helps me a lot. One thing is that I do not need to split flyover links when those links do not intersect with other links in actual situation (I haven't tried yet but 2. Create a noded network splits all links at all intersections, isn't it? Sorry if I misunderstood).

    After simplifing and merging all links, I want to check location of those remaining overshoot links (i.e., extracting those links or generating points at intersections) so that I can check whether those links are grade-separated or not in real road network.

    Sorry if my explanation is not sufficient.

    – HSJ Jun 10 '22 at 09:24
  • 1
    @HSJ gotcha. It sounded like you wanted a fully noded network. If at this point only detection is what you need, see dr_jts's answer. – geozelot Jun 10 '22 at 09:30
  • I am trying to understand how I can utilize your suggestions but have two quwstions:

    #1. pgr_contraction returns a table with columns type, id, contracted_vertices, source, target, cost but it doesn't contain geom. How to extract the contracted network in shp file from this?

    #2. silimar to #1. pgr_createTopology returns <network>_noded_vertices_pgr but it is just a set of nodes with geom. I do not understand how to use this pgRouting topology for my usage (export the result into shp file).

    – HSJ Jun 11 '22 at 09:18
  • I realized the above two questions are out of the original scope. Please ignore it. I will make a new post.if necessary. – HSJ Jun 11 '22 at 10:56
  • 1
    @HSJ the actual utilization of pgr_contraction is indeed rather complex - there is a guide in the docs, though. The <network_simple>_noded table holds your edges, you can export them. – geozelot Jun 11 '22 at 13:46
3

If your question is as in your title, you just want to detect the lines intersected by line 1 (but not the ones just touching it) :

SELECT a.*
FROM net a, net b
WHERE a.id != 1 AND b.id = 1 AND ST_Intersects(a.geom, b.geom)
   AND ST_Touches(a.geom, b.geom) IS NOT TRUE

If you want to generate points on theses intersections as shown in your image, you have to extract a point (ST_CollectionExtract(geom,1)) from the intersection of the lines (ST_Intersection(line1,line2))

SELECT ST_CollectionExtract(ST_Intersection(a.geom, b.geom),1) as geom
FROM net a, net b
WHERE a.id != 1 AND b.id = 1 AND ST_Intersects(a.geom, b.geom)
   AND ST_Touches(a.geom, b.geom) IS NOT TRUE
Cupain
  • 695
  • 4
  • 17
  • It works well in the network that I posted as an example, but the actual network that I am using in my research contains dozens of flyovers... – HSJ Jun 10 '22 at 09:26
  • 1
    And so ? It should work for all the flyovers. But clearly, @dr_jts' answer is more synthetic and efficient – Cupain Jun 10 '22 at 09:30
3

ST_Crosses tests if two lines cross rather than touch. This can be run to check all unordered pairs using a "triangle join". The intersection point is computed with ST_Intersection.

WITH data(id, geom) AS (VALUES
  (1,ST_GeomFromText('MULTILINESTRING((136.63964778201967 36.58031552554121,136.64078637948796 36.57968013023401,136.6414207216891 36.579313967665655,136.64203092428795 36.578959650067986,136.6428838673968 36.57843301697034,136.6438347098042 36.577844992552514))', 4326)),
  (2,ST_GeomFromText('MULTILINESTRING((136.64039880046448 36.57933335255234,136.64078637948796 36.57968013023401,136.64119407544638 36.580076445271914,136.64165541506554 36.58051798901414))', 4326)),
  (3,ST_GeomFromText('MULTILINESTRING((136.64039880046448 36.57933335255234,136.64044439789086 36.578869185464725))', 4326)),
  (4,ST_GeomFromText('MULTILINESTRING((136.6409003730539 36.57873995108804,136.6414207216891 36.579313967665655,136.64222404380473 36.580108753416425))', 4326)),
  (5,ST_GeomFromText('MULTILINESTRING((136.6416232292288 36.57843409435816,136.64203092428795 36.578959650067986))', 4326)),
  (6,ST_GeomFromText('MULTILINESTRING((136.64203092428795 36.578959650067986,136.64244666738 36.57935489221467,136.64283290461503 36.57974259264671))', 4326)),
  (7,ST_GeomFromText('MULTILINESTRING((136.64250969906357 36.57802376968141,136.6428838673968 36.57843301697034,136.64326608196427 36.578854108330574))', 4326)),
  (8,ST_GeomFromText('MULTILINESTRING((136.64364829653186 36.577513284810436,136.6438347098042 36.577844992552514))', 4326)),
  (9,ST_GeomFromText('MULTILINESTRING((136.6438347098042 36.577844992552514,136.64412438862985 36.57833070559775))', 4326))
)
SELECT a.id, b.id, ST_Intersection(a.geom, b.geom) AS geom 
FROM data a
JOIN data b ON ST_Intersects(a.geom, b.geom)
WHERE a.id < b.id AND ST_Crosses(a.geom, b.geom);

See also this answer: https://gis.stackexchange.com/a/431247/14766

dr_jts
  • 5,038
  • 18
  • 15