Introduction to pgRouting
Regina Obe
Buy our books! https://postgis.us/page_buy_book
Latest books
|
 |
|
 |
What we'll cover
Will use pgRouting 3.8+ syntax, most should work in lower versions. Checkout the workshops pgRouting workshops
pgRouting deals with graphs
pgRouting is a bag of graph functions. Graphs are composed of edges connecting nodes (stop and starting positions).
Unweighted |
Directed and Weighted |
Install in PostgreSQL db
-- will also install postgis if you don't have it already
CREATE EXTENSION pgrouting CASCADE;
-- check version you are running
SELECT * FROM pgr_full_version();Prepping data for pgRouting
Tools that load OSM data into routable format
- osm2pgrouting
- https://workshop.pgrouting.org/dev/en/basic/data.html - osm2po https://osm2po.de/
Prepping data for routability
Loading any data e.g. using osm2pgsql, ogr2ogr, ogr_fdw postgresql extension, shp2pgsql, etc
Once loaded: Use following pgrouting functions to make it routable
- pgr_extractVertices - generate nodes table (from an edge geometry or source and target ids of edges table)
- pgr_separateCrossing break edges that cross
- pgr_separateTouching properly connect geometries at end points
Creating a topology for your road network
First add columns to hold node ids and costs
ALTER TABLE nyc_streets
ADD COLUMN source bigint
,ADD COLUMN target bigint
, ADD COLUMN cost float8
, ADD COLUMN reverse_cost float8;
pgRouting only works with LINESTRING
ALTER TABLE nyc_streets
ALTER COLUMN geom TYPE geometry(LINESTRING,26918) USING ST_geometryN(geom,1);Using pgr_ExtractVertices with geometry
CREATE TABLE nyc_streets_vertices_pgr AS
SELECT id, in_edges, out_edges, x, y,
geom::geometry(POINT, 26918) AS geom
FROM pgr_extractVertices(
'SELECT gid AS id, geom
FROM nyc_streets
ORDER BY id') ;
CREATE UNIQUE INDEX ux_nyc_streets_vertices_pgr ON nyc_streets_vertices_pgr(id);
CREATE INDEX ix_nyc_streets_vertices_pgr_geom
ON nyc_streets_vertices_pgr USING gist(geom);Update your edges table with source and target columns
/* -- set the source information */
WITH
out_going AS (
SELECT id AS vid, unnest(out_edges) AS eid
FROM nyc_streets_vertices_pgr
)
UPDATE nyc_streets
SET source = vid
FROM out_going WHERE gid = eid;
/* -- set the target information */
WITH
in_coming AS (
SELECT id AS vid, unnest(in_edges) AS eid
FROM nyc_streets_vertices_pgr
)
UPDATE nyc_streets
SET target = vid
FROM in_coming WHERE gid = eid;SELECT * INTO dc_ways_vertices_pgr_40
FROM pgr_extractVertices('SELECT gid AS id, source, target
FROM dc_ways ORDER BY gid');
UPDATE dc_ways_vertices_pgr_40 AS n
SET geom = o.the_geom, x = o.lon, y = o.lat
FROM dc_ways_vertices_pgr AS o
WHERE n.geom IS NULL AND o.id = n.id;
ALTER TABLE dc_ways_vertices_pgr_40
ADD CONSTRAINT pk_dc_ways_vertices_pgr_40 PRIMARY KEY(id);Output
SELECT * FROM dc_ways_vertices_pgr_40 LIMIT 2; id | in_edges | out_edges | x | y | geom
----+-----------+--------------------+-------------+------------+---------------------------------------
1 | {21305} | {7174,21306,22045} | -77.0034402 | 38.8842085 | 0101000020E6...
2 | {7,21306} | {7173,21307} | -77.003434 | 38.8833257 | 0101000020E...
(2 rows)
Checking graph connectedness
Use component family of functions. A component is a set of nodes that are reachable from each other.
- pgr_connectedComponents - nodes reachable from each other in an undirected graph.
- pgr_strongComponents nodes reachable in directed way
- Other component family of functions - https://docs.pgrouting.org/dev/en/components-family.html
Strongly connected
SELECT component, count(node) AS num_nodes
FROM pgr_strongComponents(
'SELECT gid AS id, source, target, cost, reverse_cost
FROM dc_ways'
)
GROUP BY component
ORDER BY num_nodes DESC
LIMIT 5; component | num_nodes
-----------+-----------
1 | 85454
19809 | 68
65644 | 54
7292 | 46
8121 | 45
(5 rows)
ALTER TABLE dc_ways_vertices_pgr_40 ADD COLUMN component bigint;
UPDATE dc_ways_vertices_pgr_40 AS n SET component = s.component
FROM pgr_strongComponents(
'SELECT gid AS id, source, target, cost, reverse_cost FROM dc_ways'
) AS s
WHERE s.node = n.id;Mark off your components for future use in edges
ALTER TABLE dc_ways_vertices_pgr_40 ADD COLUMN components bigint[];
-- source component
UPDATE dc_ways AS e SET components = ARRAY[n.component]
FROM dc_ways_vertices_pgr_40 AS n
WHERE n.id = e.source;
-- add target component if different from source
UPDATE dc_ways AS e
SET components = COALESCE(e.components, e.components || n.component)
FROM dc_ways_vertices_pgr_40 AS n
WHERE n.id = e.target AND COALESCE(e.components[1],0) <> n.component
AND n.component > 0;Basic routing using pgr_dijkstra
SELECT e.name, dj.seq, dj.cost, dj.agg_cost, e.the_geom AS geom
FROM pgr_dijkstra('SELECT gid AS id, source, target,
cost_s AS cost, reverse_cost_s AS reverse_cost
FROM dc_ways', 33998, 56395 ,directed=> true ) AS dj
INNER JOIN dc_ways AS e on dj.edge = e.gid
ORDER BY dj.seq;name | seq | cost | agg_cost ------------------------+-----+--------------------+-------------------- Potomac Heritage Trail | 1 | 598.6978157236285 | 0 Potomac Heritage Trail | 2 | 106.94874627527133 | 598.6978157236285 Potomac Heritage Trail | 3 | 58.0141171218561 | 705.646561998899 : : Mount Vernon Trail | 31 | 1.8059528318046105 | 1473.9401035418382 Mount Vernon Trail | 32 | 34.968050404314305 | 1475.74605637364
-- determine closest node ids connected
WITH a AS ( SELECT array_agg(ns.id ORDER BY f.sort) AS nodes
FROM ( VALUES ( 'Hotel', ST_Point(-77.23212, 38.9551, 4326), 1 )
, ('Pizza', ST_Point(-77.0842, 38.8392, 4326), 2 ) )
AS f(name, geom, sort)
CROSS JOIN LATERAL (SELECT id
FROM dc_ways_vertices_pgr_40 AS n
WHERE component = 1
ORDER BY n.geom <-> f.geom LIMIT 1) AS ns
)
SELECT e.name, dj.seq, dj.cost, dj.agg_cost, e.the_geom AS geom
FROM a, pgr_dijkstra('SELECT gid AS id, source, target,
cost_s AS cost, reverse_cost_s AS reverse_cost
FROM dc_ways', a.nodes[1], a.nodes[2],directed=> true ) AS dj
INNER JOIN dc_ways AS e on dj.edge = e.gid;With Points
With Points family of functions works with virtual nodes. Refer to https://docs.pgrouting.org/latest/en/withPoints-family.html
CREATE TABLE dc_pois AS
SELECT f.id, f.name, NULL::bigint AS edge_id, NULL::float AS fraction
, NULL::character(1) AS side
, NULL::float8 AS distance, f.geom,
NULL::bigint AS pid
FROM ( VALUES
('Hotel', ST_Point(-77.23212, 38.9551, 4326), 1 ),
('Pizza', ST_Point(-77.0842, 38.8392, 4326), 2 ) )
AS f(name, geom, id) ;Compute Virtual Nodes
Use pgr_findCloseEdges, which is new in pgRoutng 3.8. cap is max number of virtual nodes you want returned for each point.
UPDATE dc_pois AS p SET (edge_id, fraction, side, distance,
pid) =
(ce.edge_id, ce.fraction, ce.side, ce.distance,
CASE
WHEN ce.fraction = 0 THEN e.source
WHEN ce.fraction = 1 THEN e.target
ELSE -p.id END )
FROM
pgr_findCloseEdges(
'SELECT gid AS id, the_geom AS geom FROM dc_ways WHERE 1 = ANY(components)',
ARRAY(SELECT geom FROM dc_pois),
0.06, cap => 1) AS ce LEFT JOIN dc_ways AS e ON ce.edge_id = e.gid
WHERE ce.geom = p.geom;Routing withPoints
SELECT *
FROM
pgr_withPoints(
'SELECT gid AS id, source, target, cost_s AS cost, reverse_cost_s AS reverse_cost
FROM dc_ways
WHERE 1 = ANY(components) ORDER BY gid',
'SELECT pid, edge_id, fraction, side FROM dc_pois WHERE pid < 0',
-- note for nodes not on network the pid will be negative
(SELECT pid FROM dc_pois WHERE name = 'Hotel'),
(SELECT pid FROM dc_pois WHERE name = 'Pizza'), details => true);Find 2 Answers and aggregate using pgr_Ksp
SELECT dj.path_id, SUM(dj.cost) AS cost
, ST_NPoints(ST_MakeLine(e.the_geom ORDER BY dj.path_seq)) AS npoints
FROM a, pgr_ksp('SELECT gid AS id, source, target,
cost_s AS cost, reverse_cost_s AS reverse_cost
FROM dc_ways', 33998, 56395, 2, -- max number of results
directed=> true ) AS dj
INNER JOIN dc_ways AS e on dj.edge = e.gid
GROUP BY dj.path_id;path_id | cost | npoints
---------+-------------------+---------
1 | 1510.714106777957 | 745
2 | 1511.106210155594 | 752
(2 rows)
Computing service areas
- pgr_DrivingDistance - for nodes on network
- pgr_withPointsDD - if you need to include virtual nodes
Follow this with using PostGIS ST_ConcaveHull or postgis_sfcgal ST_AlphaShape
equicost => true will return at most one node for each item in an array.It will return the item with the lowest cost.
SELECT *
FROM pgr_drivingDistance(
'SELECT gid AS id, source, target, cost_s AS cost, reverse_cost_s AS reverse_cost
FROM dc_ways WHERE 1 =ANY(components)',
array[33998], 60*20, equicost => true);seq | depth | start_vid | pred | node | edge | cost | agg_cost -----+-------+-----------+-------+-------+--------+---------------------+-------------------- 1 | 0 | 33998 | 33998 | 33998 | -1 | 0 | 0 2 | 1 | 33998 | 33998 | 23290 | 25632 | 598.6978157236285 | 598.6978157236285 3 | 2 | 33998 | 23290 | 27888 | 30503 | 106.94874627527133 | 705.6465619988999 4 | 3 | 33998 | 27888 | 32754 | 35477 | 58.0141171218561 | 763.660679120756 : :
SELECT ST_ConcaveHull(ST_Collect(n.the_geom), 0.7, true)
FROM pgr_drivingDistance(
'SELECT gid AS id, source, target, cost_s AS cost,
reverse_cost_s AS reverse_cost
FROM dc_ways WHERE 1 =ANY(components)',
array[33998], 60*20, equicost => true) AS dd
INNER JOIN dc_ways_vertices_pgr AS n ON dd.node = n.id;