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

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:

• You should check OSRM or Graphhopper. 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 ? Jan 14 at 12:41
• @geozelot Could you help me out with this problem ? Jan 19 at 0:48
• 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. Jan 20 at 14:25

3 Answers

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';
``````
• @George is this approach working ok for you? Jan 21 at 14:04
• I just don't want to run the 1-many query 3000 times. I want the output in 1 query if it is faster than the loop. It is actually many to many but not like a CrossJoin more like an InnerJoin. I will add a snapshot of the output I want from the query in the question. Also I would like to add that you have been of great help from the beginning, thanks a lot, I have given you the well deserved bounty. Jan 24 at 11:57
• Ok, I understand what you are trying to achieve. The only way I see how you can get this is looping over all 1:Many pairs and using the same logic as up top. This will not make it slower. It is perfectly ok to do so. I will change the logic now and update you once done. Jan 24 at 13:24
• I updated my answer; this is tested with pgrouting 3.3.0. Jan 24 at 20:09
• Hi @GeorgeV did you have any luck with the refactored function? Jan 30 at 20:13

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.

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.

• 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 ? 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. Jan 15 at 7:48
• 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
• You just have to make sure vertices exist wherever edges touch or intersect each other. Jan 16 at 9:26