pgRouting computation time & code efficiency

Question

I am trying to find the closest landfill by distance from various locations, however, the computation time is 8 hr 36 min to find the shortest path from 3602 locations to 5 landfill locations.

Before using pgRouting I used QNEAT3 on another larger data set and it didn't take nearly as long.

As a test, I scaled down the data to 9 locations and 5 landfill locations and it completed the shortest path analysis in ~10 min.

Is there a problem with the code and its efficiency?

-- build road topology
ALTER TABLE landfill.road ADD COLUMN "source" integer; -- add source column for pg routing
ALTER TABLE landfill.road ADD COLUMN "target" integer; -- add target column for pg routing
SELECT pgr_createTopology('landfill.road', 0.001, 'geom', 'id'); -- prepare roads layer by building network topology based on geometry information.

-- add pgr networking information
ALTER TABLE landfill.road ADD COLUMN length float8; -- add length field for pgr cost
ALTER TABLE landfill.road ADD COLUMN cost float8; -- add cost column
ALTER TABLE landfill.road ADD COLUMN reverse_cost float8; -- add reverse_cost column
UPDATE landfill.road SET length = ST_Length(geom::geography);

-- replace null values with 0
UPDATE landfill.road SET oneway=0 WHERE oneway IS null; 

-- create cost and reverse costs for road network based on ONEWAY field
UPDATE landfill.road
SET cost = ST_Length(geom::geography) 
WHERE oneway IN ('0','FT'); -- both direction, road direction and the digitizing direction are the same

UPDATE landfill.road
SET cost = 999999999
WHERE oneway = 'TF'; -- The road direction and the digitizing direction are opposite.

UPDATE landfill.road
SET reverse_cost = ST_Length(geom::geography)
WHERE oneway IN ('0', 'TF');

UPDATE landfill.road
SET reverse_cost = 999999999
WHERE oneway = 'FT';

-- Begin Dijkstra algorithm shortest path query
CREATE TABLE landfill.pgr_routes AS 
WITH all_pairs AS (
  -- all pairs of start and end geometries with IDs
  -- that get carried through so the routing results
  -- match with the pt IDs you know.
  SELECT f.id AS fid, f.geom as fgeom,
         t.id as tid, t.geom as tgeom
    FROM public.from_pts AS f,
         landfill.to_pts AS t
), vertices AS (
  SELECT fid, tid,
       (SELECT id -- proximity search for closest from vertex
          FROM landfill.road_vertices_pgr
         ORDER BY the_geom <-> fgeom
         LIMIT 1) as fv,
       (SELECT id -- proximity search for closest to vertex
          FROM landfill.road_vertices_pgr
         ORDER BY the_geom <-> tgeom
         LIMIT 1) as tv
  FROM all_pairs
), pgr_result AS (
  SELECT fid, tid, pgr_Dijkstra(
    'SELECT id, source, target, length AS cost, reverse_cost FROM landfill.road',
    fv, tv, 
    directed := true
  ) from vertices
)
SELECT fid, tid, (pgr_dijkstra).* FROM pgr_result -- to unpack all fields in the composite column with each as a seperate column. 
WHERE (pgr_dijkstra).edge = -1;

-- min distance value of agg_cost results
CREATE TABLE landfill.pgr_routes_min AS
SELECT DISTINCT ON (fid) fid, tid, agg_cost
FROM landfill.pgr_routes
WHERE agg_cost > 0  -- landfill must be greater than 0 m from location
ORDER BY fid, agg_cost;

-- Left join min distance results to locations and original landfill locations table
CREATE TABLE landfill.nearest_landfill AS
SELECT public.from_pts.*,landfill.pgr_routes_min.*,landfill.to_pts.name
FROM landfill.pgr_routes_min 
LEFT JOIN public.from_pts
ON landfill.pgr_routes_min.fid = public.from_pts.id 
LEFT JOIN landfill.to_pts
ON landfill.pgr_routes_min.tid = landfill.to_pts.id;

Answer 1

Instead of calling pgr_dijkstra multiple times you should be invoking it with 2 arrays for sources and targets and not single values (thank you @DavidP for highlighting this). This way the graph will be loaded into memory one time only. I tested the following many-to-many query on a 4 core and 16GB laptop on a topology containing 3.2M edges (OSM North California processed with osm2po) and it took 22 seconds to compute. The only difference is the array_agg() call. You can read more about it here. There are a few additional handy speedups you can implement listed here.

First of all I am creating some random points in the target area.

CREATE TABLE from_pts (
  id serial primary key,
  geom geometry(point, 4326)
);
INSERT INTO from_pts(geom)
VALUES 
(ST_MakePoint(-122.62390136718749, 38.24842651622814)),
(ST_MakePoint(-122.684326171875, 38.285624966683756)),
(ST_MakePoint(-122.78114318847655, 38.25543637637947)),
(ST_MakePoint(-122.89512634277345, 38.313107227858914));

CREATE TABLE to_pts (
  id serial primary key,
  geom geometry(point, 4326)
);
INSERT INTO to_pts(geom)
VALUES 
(ST_MakePoint(-121.58020019531249, 39.12792964388499)),
(ST_MakePoint(-121.38519287109375, 38.989302551359515)),
(ST_MakePoint(-121.56372070312499, 39.51887357127223)),
(ST_MakePoint(-122.00042724609374, 39.191819549771694));

Afterwards I'm basically using your logic with the small change mentioned up top.

WITH all_pairs AS (
  SELECT f.id AS fid, f.geom as fgeom,
         t.id as tid, t.geom as tgeom
    FROM from_pts AS f,
         to_pts AS t
),
vertices AS (
  SELECT fid, tid,
     (SELECT id
        FROM norcal_2po_4pgr
       ORDER BY geom_way <-> fgeom
       LIMIT 1) as fv,
     (SELECT id
        FROM norcal_2po_4pgr
       ORDER BY geom_way <-> tgeom
       LIMIT 1) as tv
  FROM all_pairs
), pgr_result AS (
  SELECT pgr_dijkstra(
    'SELECT id, source, target, cost, reverse_cost FROM norcal_2po_4pgr',
    array_agg(fv), array_agg(tv), 
    directed := true
  ) from vertices
)
SELECT (pgr_dijkstra).* FROM pgr_result;

To be able to obtain the user defined source and destination IDs, I would suggest to use a temporary lookup table and then join it with the pgr results.

CREATE TEMP TABLE vertices_lookup 
AS             
WITH all_pairs AS (
  SELECT f.id AS fid, f.geom as fgeom,
         t.id as tid, t.geom as tgeom
    FROM from_pts AS f,
         to_pts AS t
),
vertices AS (
  SELECT fid, tid,
     (SELECT id
        FROM norcal_2po_4pgr
       ORDER BY geom_way <-> fgeom
       LIMIT 1) as fv,
     (SELECT id
        FROM norcal_2po_4pgr
       ORDER BY geom_way <-> tgeom
       LIMIT 1) as tv
  FROM all_pairs
)
SELECT * FROM vertices;


WITH pgr_result AS (
  SELECT pgr_dijkstra(
    'SELECT id, source, target, cost, reverse_cost FROM norcal_2po_4pgr',
    array_agg(fv), array_agg(tv), 
    directed := true
  ) FROM vertices_lookup
)
SELECT (pgr_dijkstra).*, a.fid, a.tid FROM pgr_result
JOIN vertices_lookup a
ON (pgr_dijkstra).start_vid = a.fv
AND (pgr_dijkstra).end_vid = a.tv;