Drive Time Isochrones – An Example Using Finnish Airports

Site analyses can benefit greatly from using “drive-time” isochrones to define the study area. Drive time isochrones are often significantly different from simple buffer areas which disregard natural barriers such as rivers or slow roads.

Of course, creating drive time isochrones requires more input data and more compute-intensive algorithms than a simple buffer analysis. It is necessary to create a routable network graph with adequate weights to be used by the routing algorithm.

One of the most popular routing  applications in the open source world is pgRouting for PostGIS enabled databases. I’ve already shown how to create drive time isochrones for one network node based on pgRouting and QGIS.  Today, I’ll show how to create drive time isochrones for a set of points – in this case all airports in Finland.

The first step is to find the closest network node to every airport:

ALTER TABLE airport
   ADD COLUMN nearest_node integer;

CREATE TABLE temp AS
   SELECT a.gid, b.id, min(a.dist)
   FROM
     (SELECT airport.gid, 
             min(distance(airport.the_geom, node.the_geom)) AS dist
      FROM airport, node
      GROUP BY airport.gid) AS a,
     (SELECT airport.gid, node.id, 
             distance(airport.the_geom, node.the_geom) AS dist
      FROM airport, node) AS b
   WHERE a.dist = b. dist
         AND a.gid = b.gid
   GROUP BY a.gid, b.id;

UPDATE airport
SET nearest_node = 
   (SELECT id 
    FROM temp
    WHERE temp.gid = airport.gid);

Then, we can calculate drive times between network nodes and “airport nodes”. I am still looking for the most efficient way to perform this calculation. The trivial solution I used for this example was to calculate all drive time values separately for each airport node (as described in “Creating Catchment Areas with pgRouting and QGIS”).
I create the table with the first node:

CREATE TABLE catchment_airport AS
SELECT 
    id,
    the_geom,
    (SELECT sum(cost) FROM (
	   SELECT * FROM shortest_path('
	   SELECT gid AS id,
		  start_id::int4 AS source,
		  end_id::int4 AS target,
		  traveltime::float8 AS cost
	   FROM network',
	   5657,
	   id,
	   false,
	   false)) AS foo ) AS cost
FROM node;

Every following node is done with an INSERT:

INSERT INTO catchment_airport (
SELECT 
    id,
    the_geom,
    (SELECT sum(cost) FROM (
	   SELECT * FROM shortest_path('
	   SELECT gid AS id,
		  start_id::int4 AS source,
		  end_id::int4 AS target,
		  traveltime::float8 AS cost
	   FROM network',
	   123, -- put a new index here!
	   id,
	   false,
	   false)) AS foo ) AS cost
FROM node);

Afterwards, I combined the values to find the minimum drive time for each network node:

CREATE table catchment_airport_final AS
SELECT id, the_geom, min (cost) AS cost
FROM catchment_airport
GROUP By id, the_geom;

The resulting point layer was imported into QGIS. Using TIN interpolation (from Interpolation plugin), you can calculate a continuous cost surface. And Contour function (from GDALTools) finally yields drive time isochrones.

Drive time isochrones around airports in northern Finland - spatial data © National Land Survey of Finland 2011

Based on this analysis, it is possible to determine how many inhabitants live within one hour driving distance from an airport or how many people have to drive longer than e.g. ninety minutes to reach any airport.

18 comments
  1. How does this compare speed-wise against an ArcGIS implementation? I am looking for something faster than ArcGIS network analyst, and this could be promising.

    Thank You!

    • underdark said:

      As I have no access to ArcGIS, I can’t really answer your question. But it would certainly be interesting to perform such a comparison. A well-tuned PostGIS setup and a more sophisticated approach to merging catchment areas (if that’s of interest to your task) should be able to compete well.

      Have you checked for comparisons between Network Analyst and pgRouting yet?

      • I’m trying this out now on one of my virtual machines. Thanks for writing this- incredibly helpful

  2. allfab said:

    Interesting post !!

    Indeed, I don’t know how you realize your isochrones !!

    “The resulting point layer was imported into QGIS. Using TIN interpolation (from Interpolation plugin), you can calculate a continuous cost surface. And Contour function (from GDALTools) finally yields drive time isochrones.”

    What’s QGIS TIN interpolation ? How do you use Contour function with QGIS GDALTools ?!

    Thanks

    • underdark said:

      You can create a TIN (wikipedia article) using “Interpolation” plugin. (You might have to activate it in Plugin Manager.) “Contour” is located in Raster menu. Which step of creating contours from a raster causes problems?

  3. allfab said:

    Ok, Thanks for the response !!

    Finally, I find process after the post !!

  4. Taro said:

    Thank you for your excellent guidance of pgRouting on QGIS!!!

    Now, I have a question which is indirectly relation to your post.
    Using your posts as reference, we are trying to identify where is the most suitable point for firehouse in our town with OpenStreeMap, QGIS and pgrouting.
    As deeper as we understand the functional capability of this set of GIS system, we have found that they don’t have the function to determine where is the perfect coodinates :In other hand, they don’t have any optimization engine.

    After gathering information about optimization, the type of this problem is classified as p-median problem. We found an operations research locational analysis package named ‘orloca’ which work with R and we have already succeeded to find the point with it; however, this package provide only one-point answer(Fermat–Weber point or 1-median)… we need more than one-point answer.

    Now, let me ask you about the software or script to solve the problem of p-median problem(facility location problem). If you have some information about it, please let me know.

    Very sincerely.
    Taro.

    cf. Facility location
    http://en.wikipedia.org/wiki/Facility_location

    • Ilya said:

      Hi Taro,

      You could try out FLOWMAP wich has Service Location Modeling inculded. This is a free T-GIS package from Utrecht University. Personally I could not manage FLOWMAP to work with my data since it has quite specific requirements for dbf-files, but may be you would be luckier than I am. (http://flowmap.geog.uu.nl/download.php).

      Also, you could give a chance to an old LoLa package, which is unsupported already, but has some pre-compiled binaries that may help. As far as I remember it could solve p-median problems. http://www.mathematik.uni-kl.de/~lola/download.html

      Another option is to use generic MIP/LP-solver and feed it with adjacency matrix also as optimization function & constraints. Rather good one is GLPK which has a simplistic IDE called gusek http://gusek.sourceforge.net/gusek.html. It has AMPL-like syntax, so if you are familiar with Operations Research it would not be a problem.

      Good luck,
      Ilya.

  5. ludgerx said:

    Hi,
    may someone explain to me the step “Afterwards, I combined the values to find the minimum drive time for each network node”? Do I have to make a UNION with the resulting tables?
    Another question is if there is another solution than the “trivial solution”? It would be nice to make it in one step, without the calculation for each airport separatley.

    • underdark said:

      The query for determining minimum drive time for each network node is given in the post:

      CREATE table catchment_airport_final AS
      SELECT id, the_geom, min (cost) AS cost
      FROM catchment_airport
      GROUP By id, the_geom;

      I have been working on improving the “trivial” solution, but haven’t found the time to publish it yet.

      • ludgerx said:

        Thanks for your prompt response.

        What I don’t understand is the hollowing:
        If I calculate all drive time values separately for 3 airport nodes, then I receive 3 catchment-tables. But what is the next step? In your example you use the table ‘catchment_airport’. How did you get this table? I made a UNION with all catchment-tables and then I used the query you mentioned above. The result looks good to me, but I’m not shure.

        Thanks for your help

  6. Hans-Dieter said:

    Hi Anita,
    I would also be very interested in getting on missing information of your very interesting guidance.
    you describe following steps:
    – calculate drive times between network nodes and “airport nodes”(calculate all drive time values separately for each airport node)
    – Afterwards, I combined the values to find the minimum drive time for each network node
    Here you use a table “catchment_airport” as mentioned by user ludgerx before…
    How did you create this table. Could you perhaps post the sql-code you used to create this layer? this would be very helpful.
    Thanks!

    • underdark said:

      Hi Hans-Dieter,
      The SQL statement was always available in the linked post. I copied it here now. Maybe it is less confusing this way. The CREATE TABLE statement is run once. All other catchments are added manually (room for improvement, I know) using INSERT INTO. (I didn’t validate the syntax of this statement now, but you get the idea.)

      • Hans-Dieter said:

        Hi Anita,
        thanks a lot for this clarification. Now it makes sense to me :-)

        It would be great if it would be possible to calculate and insert the catchment-areas of all the nearest nodes ‘automatically’. This would combine the convieniece of a tool like the network-extension of gvsig with the power of postgis.
        I’ll stay tuned to your blog in case you’ll get some more spontaneous intuitions :-)

        best wishes
        and thanks for spending your time on expanding and maintaining this blog and for sharing your experience

  7. Trevor Johnson said:

    Hi. The map is now at Wikimedia Commons, and used in Isochrone map. I believe this is fine because it’s created using free software. But if you know something different, please say so. Thanks.

    • underdark said:

      Hi Trevor,

      Thanks for the notification. I’ve updated the image description page with further information. I’ve also changed the license to CC-BY-3.0 since the image is a map rather than a simple screenshot.

      Best wishes,
      Anita

  8. Trevor Johnson said:

    Great – that’s very helpful. Danke!