Using pgRouting on an OSM database imported by osm2pgsql

If you've tried to use pgRouting on an OSM database imported by osm2pgsql you've probably noticed it doesn't work.

It doesn't work because one of the requirements for pgRouting is that streets are split at intersections, meaning it expects each osm_id from planet_osm_line to only contain street segments.
However, in planet_osm_line, a street is represented by a LINESTRING, most of the times the intersection with other streets being just another POINT in the street's LINESTRING.

A visual representation of the long explanation above:







The good news is that it's fairly easy to split the OSM streets into constituent segments using PLPGSQL.
One way to achieve that is the one listed below (note the filtering on "Sector 1" - I am restricting table entries within the administrative boundary of the first sector of Bucharest because things get really slow on the entire OSM road network):



drop table if exists network;
create table network(gid serial, osm_id integer, name varchar, the_geom geometry, source integer, target integer, length float);


CREATE OR REPLACE FUNCTION compute_network() RETURNS text as $$
DECLARE
streetRecord record;
wayRecord record;
pointCount integer;
pointIndex integer;
geomFragment record;
BEGIN
-- for each street
FOR streetRecord in select way, osm_id, name from planet_osm_line where highway is not null and st_contains((select way from planet_osm_polygon where boundary = 'administrative' and name like 'Sector 1' limit 1), way) LOOP
-- for each street in the region of interest
SELECT * from planet_osm_ways where id = streetRecord.osm_id into wayRecord;
FOR pointIndex in array_lower(wayRecord.nodes, 1)..array_upper(wayRecord.nodes,1)-1 LOOP
RAISE NOTICE 'Inserting name % source %, target %', streetRecord.name, wayRecord.nodes[pointIndex], wayRecord.nodes[pointIndex+1];
select st_makeline(st_pointn(streetRecord.way, pointIndex), st_pointn(streetRecord.way, pointIndex+1)) as way into geomFragment;
insert into network(osm_id, name, the_geom, source, target, length) values(streetRecord.osm_id, streetRecord.name, geomFragment.way, wayRecord.nodes[pointIndex], wayRecord.nodes[pointIndex+1], st_length(geomFragment.way));
END LOOP;
END LOOP;
return 'Done';
END;
$$ LANGUAGE 'plpgsql';

select * from compute_network();
select assign_vertex_id('network', 1, 'the_geom', 'gid');


The PLPGSQL procedure works by creating a street segment for each pair of nodes from planet_osm_ways; the geometry is created by simply creating a new line from the points corresponding to the nodes (see the st_makeline call).

Once assign_vertex_id has been run, you can start using shortest_path:



osm=# select n.name, p.edge_id, o.highway from network n, (select edge_id from shortest_path('
osm'# SELECT gid as id,
osm'# source,
osm'# target,
osm'# length as cost
osm'# FROM network',
osm(# (select source from network where name='Intrarea Jugului' limit 1),(select source from network where name='Strada Sofia' limit 1),false,false)) p,
osm-# planet_osm_line o
osm-# where p.edge_id = n.gid and n.osm_id = o.osm_id;
name | edge_id | highway
------------------------------+---------+-------------
Intrarea Jugului | 2330 | residential
Strada Elena Caragiani | 808 | residential
Strada Borșa | 2016 | residential
Strada Borșa | 2015 | residential
Strada Borșa | 2014 | residential
Strada Borșa | 2013 | residential
| 1992 | residential
Șoseaua Pipera | 8098 | secondary
Strada Nicolae Caramfil | 1828 | secondary
Strada Nicolae Caramfil | 1827 | secondary
Strada Nicolae Caramfil | 1826 | secondary
Strada Nicolae Caramfil | 1825 | secondary
Strada Nicolae Caramfil | 1824 | secondary
Strada Nicolae Caramfil | 1823 | secondary
Strada Nicolae Caramfil | 1822 | secondary
Strada Nicolae Caramfil | 1821 | secondary
Bulevardul Beijing | 1035 | secondary
Bulevardul Beijing | 1036 | secondary
Bulevardul Beijing | 7883 | secondary
Bulevardul Beijing | 7766 | secondary
Bulevardul Beijing | 1474 | residential
Bulevardul Beijing | 1473 | residential
Bulevardul Beijing | 1472 | residential
Bulevardul Beijing | 1471 | residential
Bulevardul Beijing | 1470 | residential
Bulevardul Beijing | 1469 | residential
Bulevardul Beijing | 1468 | residential
Bulevardul Beijing | 1467 | residential
Piața Charles de Gaulle | 859 | residential
Piața Charles de Gaulle | 858 | residential
Piața Charles de Gaulle | 8918 | residential
Piața Charles de Gaulle | 8917 | residential
Piața Charles de Gaulle | 8916 | residential
Piața Charles de Gaulle | 8915 | residential
Piața Charles de Gaulle | 8914 | residential
Piața Charles de Gaulle | 8913 | residential
Piața Charles de Gaulle | 8912 | residential
Calea Dorobanților | 264 | primary
| 9999 | footway
Calea Dorobanților | 8116 | primary
Calea Dorobanților | 8117 | primary
Calea Dorobanților | 8118 | primary
Calea Dorobanților | 8119 | primary
Strada Emil Pangrati | 3489 | residential
Strada Emil Pangrati | 3412 | residential
Strada Emil Pangrati | 3413 | residential
Strada Emil Pangrati | 3414 | residential
Cpt. Av. Demetriade Gheorghe | 3446 | residential
Strada Andrei Muresanu | 8086 | residential
(49 rows)

Time: 77.807 ms


If you render the above (I used QGIS), the (start of) the route looks like this:



Note how pgRouting has found multiple edges (2016, 2015, 2014, 2013) alongside "Strada Borșa"; that's because we split the street using the PLPGSQL procedure described at the beginning of the post:


osm=# select osm_id, name from network where gid in (2016, 2015, 2014, 2013);
osm_id | name
----------+--------------
23730998 | Strada Borșa
23730998 | Strada Borșa
23730998 | Strada Borșa
23730998 | Strada Borșa
(4 rows)