4

Referring to the below snapshot what is the best possible way to find the road distance from the red star node to all of the yellow nodes using PostGIS.

Get 1 to many Shortest Path

Note:

  1. I don't have the road distance data between these points.
  2. I want to achieve this as per road distance and not linear distance.
  3. The line data in the snapshot is road vector data which is saved in the database, and I have to achieve the shortest path using this road data.
  4. All Points are connected to the road vector data.

Sample Data I am working on - (Dont Use this dataset, I have updated another dataset below):

DROP TABLE IF EXISTS roads;
CREATE TABLE roads (id varchar, geom geometry);
INSERT INTO roads VALUES
    ('ROAD1', 'LINESTRING(-1.141319196255101 51.347918898902854,-1.141361485547053 51.34796102120045,-1.141541206242804 51.34814055776023,-1.141763074795662 51.34836221514714)'),
    ('ROAD2', 'LINESTRING(-1.143044819374352 51.34730362524181,-1.142182014405775 51.34761122033146,-1.141319196255101 51.347918898902854)'),
    ('ROAD3', 'LINESTRING(-1.141872927516498 51.347960873667525,-1.141797178454002 51.347748451769824)'),
    ('ROAD4', 'LINESTRING(-1.141319196255101 51.347918898902854,-1.140977707176789 51.34804381159354,-1.140714722149814 51.348139985313914,-1.14011024056886 51.34836115847704)'),
    ('ROAD5', 'LINESTRING(-1.141331098443016 51.34818591418458,-1.141458722853229 51.34805815888045)'),
    ('ROAD6', 'LINESTRING(-1.141328493341411 51.34774061970859,-1.141383846217907 51.347895844904535)');
    
DROP TABLE IF EXISTS origin;
CREATE TABLE origin (id varchar, geom geometry);
INSERT INTO origin VALUES
    ('ORIGIN1', 'POINT(-1.141566300914552 51.34816562854784)');

DROP TABLE IF EXISTS destination;
CREATE TABLE destination (id varchar, geom geometry);
INSERT INTO destination VALUES
    ('DEST1', 'POINT(-1.141872927516498 51.347960873667525)'),
    ('DEST2', 'POINT(-1.141331098443016 51.34818591418458)');

Updated Dataset which works with pgrouting:

DROP TABLE IF EXISTS roads;
CREATE TABLE roads (gid varchar, geom geometry);
INSERT INTO roads VALUES
    ('ROAD1', 'LINESTRING(-1.145755434320116 51.35288041913489,-1.145958533757461 51.353156152678686,-1.146161635617728 51.35343188586885,-1.146235049722674 51.353561544543304,-1.146248259363631 51.353584839902574,-1.146297066030328 51.35367088834458,-1.14643235592271 51.35390979970788)'),
    ('ROAD2', 'LINESTRING(-1.145755434320116 51.35288041913489,-1.145169320480553 51.35312342272015,-1.144583200458985 51.353366423334)'),
    ('ROAD3', 'LINESTRING(-1.144583200458985 51.353366423334,-1.144721848582808 51.35363269662706,-1.14486049829772 51.35389896975151)'),
    ('ROAD4', 'LINESTRING(-1.145691797718894 51.353985960013546,-1.145619669763756 51.354003417997525,-1.145547541754238 51.35402087593678,-1.145117361865354 51.35408679647938,-1.145043480663247 51.354098126822066,-1.14486049829772 51.35389896975151)'),
    ('ROAD5', 'LINESTRING(-1.146348470573612 51.352650950913706,-1.146051953184763 51.35276568540412,-1.145755434320116 51.35288041913489)'),
    ('ROAD6', 'LINESTRING(-1.146348470573612 51.352650950913706,-1.146518014664603 51.35280055232736,-1.146687559857114 51.352950153497005)'),
    ('ROAD7', 'LINESTRING(-1.147454238700725 51.353423312797986,-1.147244200155751 51.35313404352985,-1.147034164239104 51.35284477388327,-1.14686086224574 51.35289746381949,-1.146687559857114 51.352950153497005)');
    
DROP TABLE IF EXISTS origin;
CREATE TABLE origin (gid varchar, geom geometry);
INSERT INTO origin VALUES
    ('ORIGIN1', 'POINT(-1.145841376721831 51.35324471667814)');

DROP TABLE IF EXISTS destination;
CREATE TABLE destination (gid varchar, geom geometry);
INSERT INTO destination VALUES
    ('DEST1', 'POINT(-1.141872927516498 51.347960873667525)'),
    ('DEST2', 'POINT(-1.145441417998335 51.352888872841156)'),
    ('DEST3', 'POINT(-1.146277431045813 51.35276500193711)'),
    ('DEST4', 'POINT(-1.146846127121788 51.35296278541605)'),
    ('DEST5', 'POINT(-1.145622016050505 51.35410975452462)');

Queries Executed: (as per Shortest paths between points over multiple linestrings )

1)

CREATE TABLE road_segments AS
WITH dumps AS ( 
  SELECT gid, ST_DumpPoints(geom) AS pt FROM roads), 
pts AS (
  SELECT gid, (pt).geom, (pt).path[1] AS vert FROM dumps) 
SELECT a.gid as source_id, row_number() over() as gid, ST_MakeLine(ARRAY[a.geom, b.geom]) AS geom
FROM pts a, pts b 
WHERE a.gid = b.gid AND a.vert = b.vert-1 AND b.vert > 1;
ALTER TABLE road_segments ADD COLUMN IF NOT EXISTS source integer;
ALTER TABLE road_segments ADD COLUMN IF NOT EXISTS target integer;
SELECT pgr_createTopology('road_segments', 0.0001, 'geom', 'gid');
WITH start_point AS (
    SELECT geom from origin
),
destination_point AS (
    SELECT geom from destination
    WHERE gid = 'DEST1'
),
dijkstra AS (
    SELECT * FROM pgr_dijkstra(
        'SELECT gid as id, source, target, ST_Length(geom) AS cost FROM road_segments',
        -- finding the source node spatially
        (SELECT a.id FROM road_segments_vertices_pgr a, start_point b
            ORDER BY a.the_geom <-> b.geom LIMIT 1),
        -- finding the target node spatially
        (SELECT a.id FROM road_segments_vertices_pgr a, destination_point b
            ORDER BY a.the_geom <-> b.geom LIMIT 1),
        false
    )
),
route_geom AS (
    SELECT dijkstra.seq, dijkstra.cost,
    CASE
        WHEN dijkstra.node = road_segments.source THEN geom
        ELSE ST_Reverse(geom)
    END AS route_geom
    FROM dijkstra JOIN road_segments
    ON (edge = gid) ORDER BY seq
)
SELECT ST_AsText(ST_Collect(route_geom)) FROM route_geom;

Above query no 4) gives me the output for 1 to 1 records. But as there are thousands of records, 1 to 1 execution is not feasible. I want to connect 1 source to all its destinations at once. That is I want each source to destination path in individual row. (1 to many)

Below is the sample output I expect from this dataset: Sample Output

What changes should I make to this query to achieve this ?

Also the broad view of my task is to connect many to many points, where Source1 will connect Dest1-Dest5 and Source2 will connect Dest6-Dest12 and so on. I was thinking of implementing 1-many first and then looping it on each SourceNodes. But if even this (many-many)[conditional] is possible in a single query reducing the overall execution time then it will be really helpful.

4
  • 1
    You should check OSRM or Graphhopper.
    – pLumo
    Jan 14 at 12:32
  • @pLumo I have the vector road data as you can see in the snapshot. But its in bits and pieces and I have to use that data only and not through any other API. Is there a way to achieve this using PostGIS only ?
    – George V
    Jan 14 at 12:41
  • @geozelot Could you help me out with this problem ?
    – George V
    2 days ago
  • 1
    I could, but in your case I would use pgr_withPoints, which adds the need for a different table setup and different route aggregation methodology, and its implementation as an answer is beyond my capabilities in terms of time.
    – geozelot
    yesterday
1

You can make use of the dijkstraNear function recently added to pgRouting for one-to-many, many-to-one or many-to-many queries. I have written up a PL/pgSQL function based on some logic of the other thread you already posted. Note: you will have to run all steps listed up top, i.e. a topology will be required. I have commented inline for a better understanding.

CREATE OR REPLACE FUNCTION compute_many_to_many()
RETURNS SETOF RECORD AS $$
DECLARE
  
  route record;
  origins BIGINT[];
  origins_text TEXT[];
  destinations BIGINT[];
  destinations_text TEXT[];
 
BEGIN
  
    -- this finds the closest nodes to your origin points
    -- and saves both the nodes and them to an array
    -- the ids array will be input for the many-to-many query
    SELECT INTO origins_text, origins
    array_agg(origins.gid)::TEXT[] AS origins,
    array_agg(vertices.id)::BIGINT[] as ids
    FROM
      origin AS origins
    CROSS JOIN LATERAL
      (SELECT id, the_geom
        FROM road_segments_vertices_pgr
        ORDER BY origins.geom <-> the_geom
        LIMIT 1) AS vertices;
  
    -- does the same for our destinations and is used in our
    -- many-to-many query
    SELECT
    INTO destinations_text, destinations
    array_agg(destinations.gid)::TEXT[] AS destinations,
    array_agg(vertices.id)::BIGINT[] as ids
    FROM
      destination AS destinations
    CROSS JOIN LATERAL
      (SELECT id, the_geom
        FROM road_segments_vertices_pgr
        ORDER BY destinations.geom <-> the_geom
        LIMIT 1) AS vertices;
  
    -- RAISE NOTICE '% % % %', origins, origins_text, destinations, destinations_text;

    FOR route in
        SELECT 
        SUM(dn.cost)::DOUBLE PRECISION AS cost,
        start_vid::INT AS origin_node,
        -- we want the origin as text again hence the index lookup
        origins_text[array_position(origins, start_vid)] AS origin, 
        end_vid::INT AS destination_node,
        -- we want the destination as text again hence the index lookup
        destinations_text[array_position(destinations, end_vid)] AS destination,  
        -- the edge sequence
        array_agg(edge) AS edge_ids,
        -- the node sequence
        array_agg(node) AS node_ids,
        -- we collect the geometries into a multilinestring
        ST_AsText(
          ST_Collect(et.geom)
        ) AS geometries
        -- https://docs.pgrouting.org/3.2/en/pgr_dijkstraNear.html
        -- we set cap to -1 to get all routes
        FROM pgr_dijkstraNear( 
          'SELECT gid AS id, source, target, ST_Length(geom) AS cost, ST_Length(geom) AS reverse_cost FROM road_segments'::TEXT,
            origins,
            destinations,
            directed => false,
            global => true,
            cap => -1
        ) AS dn
        -- join with the road segments again to derive the geometry
        JOIN road_segments et 
        ON dn.edge = et.gid
        -- we want to group the individual steps of the route together to aggregate
        -- the information required
        GROUP BY end_vid, start_vid
    LOOP
        RETURN NEXT route;
    END LOOP;


END $$ LANGUAGE 'plpgsql';

This function returns a set of records, so you can call it this way:

SELECT * FROM compute_many_to_many() AS (
  cost DOUBLE PRECISION, 
  origin_node INT, 
  origin_text TEXT, 
  destination_node INT, 
  destination_text TEXT, 
  edge_ids BIGINT[], 
  node_ids BIGINT[], 
  geometries TEXT
);

This will output records looking like:

cost             | 0.0018078542184695255
origin_node      | 19
origin_text      | ORIGIN2
destination_node | 1
destination_text | DEST2
edge_ids         | {21,22,23,19,18,16,17}
node_ids         | {19,20,21,17,16,14,15}
geometries       | MULTILINESTRING((-1.147244200155751 51.35313404352985,-1.147034164239104 51.35284477388327),(-1.147034164239104 51.35284477388327,-1.14686086224574 51.352897463819
49),(-1.14686086224574 51.35289746381949,-1.146687559857114 51.352950153497005),(-1.146518014664603 51.35280055232736,-1.146687559857114 51.352950153497005),(-1.146348470573612 51.35
2650950913706,-1.146518014664603 51.35280055232736),(-1.146348470573612 51.352650950913706,-1.146051953184763 51.35276568540412),(-1.146051953184763 51.35276568540412,-1.145755434320
116 51.35288041913489))
1
  • @George is this approach working ok for you? 5 hours ago
0

Looking at your network example I would use pgRouting's Driving Distance function, which returns all distances from one start node to all other nodes.

If your network is large, then you should select a part of it, for example using a bounding box. you just need to make sure all the relevant nodes are included.

16
  • 1
    I basically answered a very similar question here with sql code @George V - check it here gis.stackexchange.com/questions/420360/… Jan 14 at 14:30
  • @TimothyDalton I tried to use your solution, followed all steps, there are records in both segment and segment_vertices_pgr tables, but still I am getting null as output in the final query. What could be the reason ?
    – George V
    Jan 15 at 7:24
  • @TimothyDalton I have also updated my original question with some sample data. If you can give me some working sample queries for shortest path to those 2 destination points from the origin point using the road linestrings it will be very helpful.
    – George V
    Jan 15 at 7:48
  • 1
    The reason this is not working is because your roads input data is not processed in a way that pgrouting can build a sound topology. What I would suggest is that you take your roads input file and use this handy command line tool which will break the roads into segments where they intersect each other pgrouting.org/docs/tools/osm2pgrouting.html here an example why this is not working i.stack.imgur.com/ul7iW.png the roads touch each other but where this happens the in red highlighted road has to be broken into 2 segments. Jan 15 at 10:34
  • 1
    You just have to make sure vertices exist wherever edges touch or intersect each other. Jan 16 at 9:26
0

After being confused for a lot of time. I finally managed to get the desired output for 1-many routing using the below query.

WITH start_point AS (
    SELECT geom from origin
),
 dest_vert As (
    (with vertices as (
    select ST_Union(a.the_geom) as the_geom from road_segments_vertices_pgr a)
    select b.gid,ST_AsText(geom),
    (select id from road_segments_vertices_pgr 
    where ST_AsText(the_geom)=ST_AsText(ST_ClosestPoint(vertices.the_geom,geom) ))
    from destination b,vertices
    order by b.gid)
),
dijkstra AS (
    SELECT * FROM pgr_dijkstra(
        'SELECT gid as id, source, target, ST_Length(geom) AS cost FROM road_segments',
        -- finding the source node spatially
        (SELECT a.id FROM road_segments_vertices_pgr a, start_point b
            ORDER BY a.the_geom <-> b.geom LIMIT 1),
        -- finding the target node spatially
        (SELECT array_agg(a.id) FROM dest_vert a),
        false
    )
),
route_geom AS (
    SELECT dijkstra.*,dijkstra.seq, dijkstra.cost,
    CASE
        WHEN dijkstra.node = road_segments.source THEN geom
        ELSE ST_Reverse(geom)
    END AS route_geom
    FROM dijkstra JOIN road_segments
    ON (edge = gid) ORDER BY seq
)
SELECT gid as destination_id,ST_AsText(ST_Collect(route_geom)) as route_geom,
ST_Length(ST_Collect(route_geom)::geography) as route_length 
FROM route_geom inner join dest_vert on dest_vert.id=route_geom.end_vid
group by gid order by gid

It works, although I am not exactly sure if this is the right way to do it. I really hope this helps somebody.

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.