Calculating shortest route by road from one origin point to multiple destination points at once using PostGIS (pgrouting)

Question

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.

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:

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-Dest3 and Source2 will connect Dest4-Dest5 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.

Many-to-Many Dataset Changes:

INSERT INTO origin VALUES
    ('ORIGIN2', 'POINT(-1.146456933398265 51.35255457699849)');
    
DROP TABLE IF EXISTS destination;
CREATE TABLE destination (gid varchar, geom geometry,origin_id varchar); 

INSERT INTO destination VALUES
    ('DEST1', 'POINT(-1.141872927516498 51.347960873667525)','ORIGIN1'),
    ('DEST2', 'POINT(-1.145441417998335 51.352888872841156)','ORIGIN1'),
    ('DEST3', 'POINT(-1.146277431045813 51.35276500193711)','ORIGIN1'),
    ('DEST4', 'POINT(-1.146846127121788 51.35296278541605)','ORIGIN2'),
    ('DEST5', 'POINT(-1.145622016050505 51.35410975452462)','ORIGIN2');

Expected Many to Many Output:

Answer 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))

### Update to accommodate OP expectation of mappings

Different to the above approach, this will conduct a set of one-to-many calls as not all combinations are required. It basically uses the same logic as before but adds one extra loop iterating over all origin points.

SELECT * FROM compute_many_one_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
);

CREATE OR REPLACE FUNCTION compute_many_one_to_many()
RETURNS SETOF RECORD AS $$
DECLARE
  route record;
  
  ogn_id TEXT;
  ogn_nn_node_id BIGINT;

  destinations BIGINT[];
  destinations_text TEXT[];
  
BEGIN
    
  -- in a loop of distinct origins
  -- this finds the closest node to each
  -- this is the "one"
  FOR ogn_id, ogn_nn_node_id IN 
    SELECT 
      DISTINCT(gid),
      vertices.id::BIGINT   
    FROM origin
    CROSS JOIN LATERAL
      (SELECT id, the_geom
        FROM road_segments_vertices_pgr
        ORDER BY origin.geom <-> the_geom
        LIMIT 1) AS vertices

  LOOP

    -- does the same for our destinations 
    -- and is used as the "many"
    SELECT
    INTO destinations_text, destinations
    array_agg(d.gid)::TEXT[] AS destinations,
    array_agg(vertices.id)::BIGINT[] as ids
    FROM
      destination AS d
    CROSS JOIN LATERAL
      (SELECT id, the_geom
        FROM road_segments_vertices_pgr
        ORDER BY d.geom <-> the_geom
        LIMIT 1) AS vertices
    -- this makes sure of the mappings
    WHERE d.origin_id = ogn_id;


    -- RAISE NOTICE '% % % %', ogn_id, ogn_nn_node_id, destinations, destinations_text;

    FOR route IN
        SELECT 
        SUM(dn.cost)::DOUBLE PRECISION AS cost,
        start_vid::INT AS origin_node,
        ogn_id 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,
            ogn_nn_node_id,
            destinations,
            directed => false,
            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 LOOP;

END $$ LANGUAGE 'plpgsql';

Answer 2

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.

Answer 3

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.