Public transport isochrones with pgRouting

This post covers a simple approach to calculating isochrones in a public transport network using pgRouting and QGIS.

For this example, I’m using the public transport network of Vienna which is loaded into a pgRouting-enable database as network.publictransport. To create the routable network run:

select pgr_createTopology('network.publictransport', 0.0005, 'geom', 'id');

Note that the tolerance parameter 0.0005 (units are degrees) controls how far link start and end points can be apart and still be considered as the same topological network node.

To create a view with the network nodes run:

create or replace view network.publictransport_nodes as
select id, st_centroid(st_collect(pt)) as geom
from (
	(select source as id, st_startpoint(geom) as pt
	from network.publictransport
	) 
union
	(select target as id, st_endpoint(geom) as pt
	from network.publictransport
	) 
) as foo
group by id;

To calculate isochrones, we need a cost attribute for our network links. To calculate travel times for each link, I used speed averages: 15 km/h for buses and trams and 32km/h for metro lines (similar to data published by the city of Vienna).

alter table network.publictransport add column length_m integer;
update network.publictransport set length_m = st_length(st_transform(geom,31287));

alter table network.publictransport add column traveltime_min double precision;
update network.publictransport set traveltime_min = length_m  / 15000.0 * 60; -- average is 15 km/h
update network.publictransport set traveltime_min = length_m  / 32000.0 * 60 where "LTYP" = '4'; -- average metro is 32 km/h

That’s all the preparations we need. Next, we can already calculate our isochrone data using pgr_drivingdistance, e.g. for network node #1:

create or replace view network.temp as
 SELECT seq, id1 AS node, id2 AS edge, cost, geom
  FROM pgr_drivingdistance(
    'SELECT id, source, target, traveltime_min as cost FROM network.publictransport',
    1, 100000, false, false
  ) as di
  JOIN network.publictransport_nodes pt
  ON di.id1 = pt.id;

The resulting view contains all network nodes which are reachable within 100,000 cost units (which are minutes in our case).

Let’s load the view into QGIS to visualize the isochrones:

isochrone_publictransport_1

The trick is to use data-defined size to calculate the different walking circles around the public transport stops. For example, we can set up 10 minute isochrones which take into account how much time was used to travel by pubic transport and show how far we can get by walking in the time that is left:

1. We want to scale the circle radius to reflect the remaining time left to walk. Therefore, enable Scale diameter in Advanced | Size scale field:

scale_diameter

2. In the Simple marker properties change size units to Map units.
3. Go to data defined properties to set up the dynamic circle size.

datadefined_size

The expression makes sure that only nodes reachable within 10 minutes are displayed. Then it calculates the remaining time (10-"cost") and assumes that we can walk 100 meters per minute which is left. It additionally multiplies by 2 since we are scaling the diameter instead of the radius.

To calculate isochrones for different start nodes, we simply update the definition of the view network.temp.

While this approach certainly has it’s limitations, it’s a good place to start learning how to create isochrones. A better solution should take into account that it takes time to change between different lines. While preparing the network, more care should to be taken to ensure that possible exchange nodes are modeled correctly. Some network links might only be usable in one direction. Not to mention that there are time tables which could be accounted for ;)

29 comments
  1. You can also use the clamp function to avoid the WHEN.

    (10 – clamp(0,”cost”, 10)) * 100 * 2

    That will return 10 if “cost” is larger then 10 or else it returns “cost”

    • Thanks! I didn’t know the clamp() function. Not sure it’s easier to read – but definitely is shorter :D

  2. Regina said:

    Anita,

    Great posts by the way. Can you tag this one PostGIS too so it appears on Planet PostGIS.

    Thanks,
    Regina

  3. Regina said:

    Anita,
    Disregard my comment. I just realized you have both a postgis tag and postgis category. I changed to the feed to use your postgis category since that looks more complete.

    • Hi Regina,
      Thanks for linking to my blog. Sorry, I realized too late that I’ve navigated myself into quite a mess with the tags and categories on this blog :)
      The PostGIS category should be quite complete.

  4. stephinova said:

    Dear Anita,

    is it correct, that the network view creates a table, where geom = NULL?
    I always received the messege, that the entry in geometry_columns is missing. So I changed the command to

    create or replace view public.publictransport5_nodes as
    select id, st_centroid(st_collect(pt))::geometry(point, 4326) as geom…

    However, when I continue with the drivingdistance command, and change the node number from 1 to 354 (which is an existing node in my publictransport table), the result table is empty as well. i.e. all fields are NULL.
    I used the following command:

    create or replace view public.temp as
    SELECT seq, id1 AS node, id2 AS edge, cost, geom
    FROM pgr_drivingdistance(
    ‘SELECT id, source, target, traveltime_min as cost FROM public.publictransport5’,
    354, 100000, false, false
    ) as di
    JOIN public.publictransport5_nodes pt
    ON di.id1 = pt.id;

    Is there anything I’m missing?

    • stephinova said:

      Hi! I’ve overlooked to change the pgr_drivingdistance to ‘SELECT gid as id, source …
      So the command pgr_drivingdistance is working.
      However I can not load the publictransport_nodes to view the isochrones..

      • What happens if you try to load the nodes? It’s hard to guess from here what might be wrong.

  5. stephinova said:

    Dear Anita!

    It took some time for me as a newby to figure it out ;)
    ST_StartPoint(multilinestring) returns NULL for PostGis 2.0. Thus I had to stick to the workaround mentioned by David William Bitner and Sandro Santilli (http://postgis.17.x6.nabble.com/ST-StartPoint-for-MultiLineStrings-td5000023.html) and changed the network view command into:

    –etc
    (Select source as id, ST_GeometryN(ST_Multi(geom),1) as pt
    from public.publictransport5
    )
    union
    (Select target as id, ST_boundary(ST_Multi(geom)) as pt
    from public.publictransport5
    — etc

    Now it works! Thanks for your great blog!

  6. madusanka said:

    I want to use pgrouting Travelling Salesman Problem (TSP) function with PHP.

    This function expects geo locations in lat/lng format and will give the optimum path for a given set of locations according to the TSP algorithm.

    I have implemented shortest path function following a tutorial on this page.

    I’m looking for a similar solution. Can somebody help me?

  7. Jonas said:

    Hi Anita,

    Thanks for a great blog! I’m studying GIS at Stockholm University where they use ArcGIS. Without your blog and other great sources on the internet I would never be able to learn gis at home with QGIS and a Mac. I’m also totally impressed with the fact that you find the time to reply to all these questions!

    I’m trying to do exactly what your doing in this post, except for the fact that I already have the speed as column SPEED and the cost as column COST calculated in my attributes for the road network. (Yes, I’m not looking a public transportation but roads). My biggest problem right now is the fact that psql is a new language and I feel totally lost here.

    1. ‘pgr_createTopology’ went great.
    2. Creating a view with the network nodes in as well. (Didn’t modify anything except my table name etc)
    3. Since I already have COST and SPEED I skipped the step where you calculated traveltimes.
    4. For the last ‘driving distance’ step my modified code was:

    create or replace view net.temp as
    SELECT seq, id1 AS node, id2 AS edge, cost, geom
    FROM pgr_drivingdistance(
    ‘SELECT id, source, target, cost as cost FROM net.vagar’,
    1, 100000, false, false
    ) as di
    JOIN net.vagar_nodes pt
    ON di.id1 = pt.id;

    From what I can tell from my map it’s working but how do I “simply update the definition of the view network.temp” to decide where to measure this?

    I wish more stuff would be possible in QGIS by visually clicking instead of psql, that way QGIS could really compete with Arc. Love you experimental plug pgRoutingLayer, thats really a step towards ordinary people using QGIS, to bad the export button isn’t working on alpha shapes yet. My plan B since I need to calculate Alpha Shapes for my project is to use pgRoutingLayer to draw the polygon and then just manually draw a new polygon on top of it from a new layer. Do you have any better suggestions?

    Many thanks again!
    Kind regards,
    Jonas

    • To update a view, simply rerun the CREATE OR REPLACE VIEW sql statement with new settings for the source node id and the max cost value.

      Enabling the export for alpha shapes is definitely on my todo list, see https://github.com/sanak/pgRoutingLayer/issues/2.

      • Jonas said:

        Had the feeling it was something really easy, but for some reason I cant figure it out. I wanna measure from source id: 51145.

        I feel so stupid. What should I change?:

        create or replace view net.temp as
        SELECT seq, id1 AS node, id2 AS edge, cost, geom
        FROM pgr_drivingdistance(
        ‘SELECT id, source, target, cost as cost FROM net.vagar’,
        1, 100000, false, false
        ) as di
        JOIN net.vagar_nodes pt
        ON di.id1 = pt.id;

      • Jonas said:

        Never mind, I sorted it out. Here’s how you do it if you’re stuck with the same problem:

        create or replace view net.test as SELECT seq, id1 AS node, id2 AS edge, cost, geom FROM pgr_drivingdistance( ‘SELECT id, source, target, cost as cost FROM net.vagar’, 51145, 100000, false, false ) as di JOIN net.vagar_nodes pt ON di.id1 = pt.id;

        How much do I have to change to create a new view with an alpha shape?

      • If you want a new view, you just have to specify a different name, so something else than net.test, e.g. net.test2 would do.

  8. Jonas said:

    Just one final question, sorry for posting this much:

    Do you have any suggestions how to this: https://anitagraser.com/2011/09/25/a-closer-look-at-alpha-shapes-in-pgrouting/ today? It seems as if the pgRouting commands have changed since 2011? And just replacing ‘points_as_polygon’ with ‘pgr_pointsaspolygon’ doesn’t do the trick.

    Have you studied GIS or are you self-taught?

    Regards,
    Jonas

    • Don’t worry, Jonas. I did a Masters in Geomatics. All open source GIS is self-taught.

      PgRouting has changed because they released 2.0 which is quite different from 1.x versions. For more recent posts check https://anitagraser.com/tag/pgrouting/.

  9. romi said:

    Hey,
    thank you for all answers from before, I solved problem with instalation. But now I’m stucked at step 5. When I cilck on retrive columns i get this error:

    syntax error at or near “;”
    LINE 7: ON di.id2 = pt.id ;

    when I click on load as new layer nothing happend, can you help me with this ? I hope you understand me

    • The post states ON di.id1 = pt.id;. Maybe that’s the problem.

      • romi said:

        I tryed that, but still the same error :
        syntax error at or near “;”
        LINE 5: ON di.id1 = pt.id;
        I would be grateful if you could help me with this.

      • Well, there is no syntax error on this line, maybe somewhere else in the query …

      • romi said:

        Ok thanks anyway, probably I did something wrong.

      • It’s only possible to help if you post the whole query.

      • romi said:

        Thanks for helping,
        This post should be reply on “Public transport isochrones with pgRouting”. I posted wrong, sorry for that…But I just copy/pasted your SQL query

        SELECT seq, id1 AS node, id2 AS edge, cost, geom
        FROM pgr_dijkstra(
        ‘SELECT id, source, target, st_length(geom) as cost FROM network.publictransport’,
        1, 3000, false, false
        ) as di
        JOIN network.publictransport pt
        ON di.id2 = pt.id;

        and still have the same error…sory for my bad english.

      • Check if the single quotes (‘) copied as straight single quotes. Try typing them again.

      • romi said:

        Typed all query but no change…thanks again for trying help.

  10. Hola! I’ve been reading your site for some time now
    and finally got the bravery to go ahead and give you a shout out from Porter Texas!
    Just wanted to say keep up the excellent job!

    • Thanks for your kind words. Best wishes