Page 1 of 1 (7 posts)

  • talks about »
  • routing

Tags

Last update:
Mon Jul 28 23:35:09 2014

A Django site.

QGIS Planet

A routing script for the Processing toolbox

Did you know that there is a network analysis library in QGIS core? It’s well hidden so far, but at least it’s documented in the PyQGIS Cookbook. The code samples from the cookbook can be used in the QGIS Python console and you can play around to get a grip of what the different steps are doing.

As a first exercise, I’ve decided to write a Processing script which will use the network analysis library to create a network-based route layer from a point layer input. You can find the result on Github.

You can get a Spatialite file with testdata from Github as well. It contains a network and a routepoints1 layer:

points_to_route1

The interface of the points_to_route tool is very simple. All it needs as an input is information about which layer should be used as a network and which layer contains the route points:

points_to_route0

The input points are considered to be ordered. The tool always routes between consecutive points.

The result is a line layer with one line feature for each point pair:

points_to_route2

The network analysis library is a really great new feature and I hope we will see a lot of tools built on top of it.


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 ;)


A Look at PgRouting find_ node_by_nearest_link_within_distance

The function with the glorious name “find_node_by_nearest_link_within_distance” is part of pgRouting and can be found in matching.sql.

“This function finds nearest node as a source or target of the nearest link”
That means that we can use this function e.g. to find the best road network node for a given address.

The function returns an object of type link_point:

CREATE TYPE link_point AS (id integer, name varchar);

To access only the id value of the nearest node, you can use:

SELECT id(foo.x) 
FROM (
   SELECT find_node_by_nearest_link_within_distance(
	'POINT(14.111 47.911)',
	0.5,
	'nw_table')::link_point as x
) AS foo


A Closer Look at Alpha Shapes in pgRouting

Alpha shapes are generalizations of the convex hull [1]. Convex hulls are well known and widely implemented in GIS systems. Alpha shapes are different in that they capture the shape of a point set. You can watch a great demo of how alpha shapes work on François Bélair’s website “Everything You Always Wanted to Know About Alpha Shapes But Were Afraid to Ask” I borrowed the following pictures from that site:

Alpha shapes for different values of alpha. The left one equals the convex hull of the point set. The right picture represents the alpha shape for a smaller value of alpha

pgRouting comes with an implementation of alpha shapes. There is an alpha shape function: alphashape(sql text) and a convenience wrapper: points_as_polygon(query character varying). The weird thing is that you don’t get to set an alpha value. The only thing supplied to the function is a set of points. Let’s see what kind of results it produces!

Starting point for this experiment is a 10 km catchment zone around node #2699 in my osm road network. Travel costs to nodes are calculated using driving_distance() function. (You can find more information on using this function in Catchment Areas with pgRouting driving_distance().)

CREATE TABLE home_catchment10km AS
SELECT *
   FROM osm_nodes
   JOIN
   (SELECT * FROM driving_distance('
      SELECT gid AS id,
          source,
          target,
          meters AS cost
      FROM osm_roads',
      2699,
      10000,
      false,
      false)) AS route
   ON
   osm_nodes.id = route.vertex_id

After costs are calculated, we can create some alpha shapes. The following queries create the table and insert an alpha shape for all points with a cost of less than 1500:

CREATE TABLE home_isodist (id serial, max_cost double precision);
SELECT AddGeometryColumn('home_isodist','the_geom',4326,'POLYGON',2);

INSERT INTO home_isodist (max_cost, the_geom) (
SELECT 1500, ST_SetSRID(the_geom,4326)
FROM 
  points_as_polygon(
    'SELECT id, ST_X(the_geom) AS x, ST_Y(the_geom) AS y FROM home_catchment10km where cost < 1500'));

In previous posts, I’ve created catchment areas by first interpolating a cost raster and creating contours from there. Now, let’s see how the two different approaches compare!

The following picture shows resulting catchment areas for 500, 1000, 1500, and 2000 meters around a central node. Colored areas show the form of pgRouting alpha shape results. Black contours show the results of the interpolation method:

Comparison of pgRouting alpha shapes and interpolation method

At first glance, results look similar enough. Alpha shape results look like a generalized version of interpolation results. I guess that it would be possible to get even closer if the alpha value could be set to a smaller value. The function should then produce a finer, more detailed polygon.

For a general overview about which areas of a network are reachable within certain costs, pgRouting alpha shapes function seems a viable alternative to the interpolation method presented in previous posts. However, the alpha value used by pgRouting seems too big to produce detailed catchment areas.

[1] http://biogeometry.duke.edu/software/alphashapes/


A Beginner’s Guide to pgRouting

The aim of this post is to describe the steps necessary to calculate routes with pgRouting. In the end, we’ll visualize the results in QGIS.

This guide assumes that you have the following installed and running:

  • Postgres with PostGIS and pgAdmin
  • QGIS with PostGIS Manager and RT Sql Layer plugins

Installing pgRouting

pgRouting can be downloaded from www.pgrouting.org.

Building from source is covered by pgRouting documentation. If you’re using Windows, download the binaries and copy the .dlls into PostGIS’ lib folder, e.g. C:\Program Files (x86)\PostgreSQL\8.4\lib.

Start pgAdmin and create a new database based on your PostGIS template. (I called mine ‘routing_template’.) Open a Query dialog, load and execute the three .sql files located in your pgRouting download (routing_core.sql, routing_core_wrappers.sql, routing_topology.sql). Congratulations, you now have a pgRouting-enabled database.

Creating a routable road network

The following description is based on the free road network published by National Land Survey of Finland (NLS). All you get is one Shapefile containing line geometries, a road type attribute and further attributes unrelated to routing.

pgRouting requires each road entry to have a start and an end node id. We’ll create those now:

First step is to load roads.shp into PostGIS. This is easy using PostGIS Manager – Data – Load Data from Shapefile.

Next, we create start and end point geometries. I used a view:

CREATE OR REPLACE VIEW road_ext AS
   SELECT *, startpoint(the_geom), endpoint(the_geom)
   FROM road;

Now, we create a table containing all the unique network nodes (start and end points) and we’ll also give them an id:

CREATE TABLE node AS
   SELECT row_number() OVER (ORDER BY foo.p)::integer AS id,
          foo.p AS the_geom
   FROM (
      SELECT DISTINCT road_ext.startpoint AS p FROM road_ext
      UNION
      SELECT DISTINCT road_ext.endpoint AS p FROM road_ext
   ) foo
   GROUP BY foo.p;

Finally, we can combine our road_ext view and node table to create the routable network table:

CREATE TABLE network AS
   SELECT a.*, b.id as start_id, c.id as end_id
   FROM road_ext AS a
      JOIN node AS b ON a.startpoint = b.the_geom
      JOIN node AS c ON a.endpoint = c.the_geom

(This can take a while.)

I recommend adding a spatial index to the resulting table.

Calculating shortest routes

Let’s try pgRouting’s Shortest Path Dijkstra method. The following query returns the route from node #1 to node #5110:

SELECT * FROM shortest_path('
   SELECT gid AS id,
          start_id::int4 AS source,
          end_id::int4 AS target,
          shape_leng::float8 AS cost
   FROM network',
1,
5110,
false,
false)

Final step: Visualization

With RT Sql Layer plugin, we can visualize the results of a query. The results will be loaded as a new layer. The query has to contain both geometry and a unique id. Therefore, we’ll join the results of the previous query with the network table containing the necessary geometries.

SELECT *
   FROM network
   JOIN
   (SELECT * FROM shortest_path('
      SELECT gid AS id,
          start_id::int4 AS source,
          end_id::int4 AS target,
          shape_leng::float8 AS cost
      FROM network',
      1,
      5110,
      false,
      false)) AS route
   ON
   network.gid = route.edge_id

In my case, this is how the result looks like:

Route from node #1 to node #5110


Routing Multiple Vehicles with pgRouting DARP Function

pgRouting has become even more powerful: A DARP (Dial-a-Ride-Problem) solver is now available in the “darp branch” of the pgRouting repository.

The Dial-a-Ride Problem (DARP) solver tries to minimize transportation cost while satisfying customer service level constraints (time windows violation, waiting and travelling times) and fleet constraints (number of cars and capacity, as well as depot location).

Documentation can be found at www.pgrouting.org/docs.


Calculating Shortest Path in QGIS Using RoadGraph Plugin

Today, Alexander Bruy announced a new QGIS plugin called RoadGraph. It is a C++ plugin that calculates the shortest path between two points on any polyline layer (e.g. Openstreetmap shapefiles).

More information can be found at GIS-Lab.

Binary files are available for Windows and Linux:

Read on: “Travelling through Brazil with Quantum GIS”


  • Page 1 of 1 ( 7 posts )
  • routing

Back to Top

Sponsors