pgRouting in OpenLayers3

This collection of notes covers my foray into using pgRouting in OpenLayers3 in EPSG:27700, or British National Grid, rather than the usual Google or Spherical Web Mercator projections. This part looks at simply tweaking the tutorial to work in EPSG:27700.

The pgRouting workshop (http://workshop.pgrouting.org/) includes a section on using your newly built network in an OpenLayers3 web application. I've combined this with an OL3 tutorial from OSGIS by Matt Walker, of Astun Technology. Matt's tutorial (http://astuntechnology.github.io/osgis-ol3-leaflet/) shows us how to build an OpenLayers3 web map using EPSG:27700 (British National Grid) web services.

For best results start with the pgRouting workshop. In chapter 8 of the pgRouting workshop you create two functions in the database to create the geometries in the format we need. We'll be working in EPSG:27700 becuse that's how I loaded my data and how I consume it.

The following function simplifies (and sets default values) when it calls the shortest path Dijkstra function. First:

-- Function: routing.pgr_dijkstra(character varying, integer, integer)
-- DROP FUNCTION routing.pgr_dijkstra(character varying, integer, integer);
CREATE OR REPLACE FUNCTION routing.pgr_dijkstra(IN tbl character varying, IN source integer, IN target integer, OUT seq integer, OUT gid integer, OUT geom geometry)
  RETURNS SETOF record AS
$BODY$
DECLARE
    sql text;
    rec record;
BEGIN
    seq := 0;
    sql := 'SELECT gid,geometry FROM ' ||
            'pgr_dijkstra(''SELECT gid as id, source::int, target::int, '
                            || 'cost_time::float AS cost FROM '
                            || quote_ident(tbl) || ''', '
                            || quote_literal(source) || ', '
                            || quote_literal(target) || ' , false, false), '
                    || quote_ident(tbl) || ' WHERE id2 = gid ORDER BY seq';

    FOR rec IN EXECUTE sql
    LOOP
        seq     := seq + 1;
        gid     := rec.gid;
        geom    := rec.geometry;
        RETURN NEXT;
    END LOOP;
    RETURN;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100
  ROWS 1000;
ALTER FUNCTION routing.pgr_dijkstra(character varying, integer, integer)
  OWNER TO postgisadmin;
COMMENT ON FUNCTION routing.pgr_dijkstra
  IS 'Shortest path function for OL3 web application';

You can test this by running a query against your network:

SELECT * FROM routing.pgr_dijkstra('itn_network',30,60);

The second function takes lat/lon points as input parameters and returns a route that can be displayed in QGIS or WMS services such as Mapserver and Geoserver:

-- Function: routing.pgr_fromatob(character varying, double precision, double precision, double precision, double precision)
-- DROP FUNCTION routing.pgr_fromatob(character varying, double precision, double precision, double precision, double precision);
CREATE OR REPLACE FUNCTION routing.pgr_fromatob(IN tbl character varying, IN x1 double precision, IN y1 double precision, IN x2 double precision, IN y2 double precision, OUT seq integer, OUT gid integer, OUT name text, OUT heading double precision, OUT cost double precision, OUT geom geometry)
  RETURNS SETOF record AS
$BODY$
DECLARE
    sql     text;
    rec     record;
    source  integer;
    target  integer;
    point   integer;

BEGIN
    -- Find nearest node
    EXECUTE 'SELECT id::integer FROM itn_network_vertices_pgr 
            ORDER BY the_geom <-> ST_GeometryFromText(''POINT(' 
            || x1 || ' ' || y1 || ')'',27700) LIMIT 1' INTO rec;
    source := rec.id;

    EXECUTE 'SELECT id::integer FROM itn_network_vertices_pgr 
            ORDER BY the_geom <-> ST_GeometryFromText(''POINT(' 
            || x2 || ' ' || y2 || ')'',27700) LIMIT 1' INTO rec;
    target := rec.id;

    -- Shortest path query (TODO: limit extent by BBOX) 
        seq := 0;
        sql := 'SELECT gid, geometry, cost, source, target, 
                ST_Reverse(geometry) AS flip_geom FROM ' ||
                        'pgr_dijkstra(''SELECT gid as id, source::int, target::int, '
                                        || 'cost_time::float AS cost FROM '
                                        || quote_ident(tbl) || ''', '
                                        || source || ', ' || target 
                                        || ' , false, false), '
                                || quote_ident(tbl) || ' WHERE id2 = gid ORDER BY seq';

    -- Remember start point
        point := source;

        FOR rec IN EXECUTE sql
        LOOP
        -- Flip geometry (if required)
        IF ( point != rec.source ) THEN
            rec.geometry := rec.flip_geom;
            point := rec.source;
        ELSE
            point := rec.target;
        END IF;

        -- Calculate heading (simplified)
        EXECUTE 'SELECT degrees( ST_Azimuth( 
                ST_StartPoint(''' || rec.geometry::text || '''),
                ST_EndPoint(''' || rec.geometry::text || ''') ) )' 
            INTO heading;

        -- Return record
            seq     := seq + 1;
            gid     := rec.gid;
            --name    := rec.name;
            cost    := rec.cost;
            geom    := rec.geometry;
            RETURN NEXT;
        END LOOP;
        RETURN;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100
  ROWS 1000;
ALTER FUNCTION routing.pgr_fromatob(character varying, double precision, double precision, double precision, double precision)
  OWNER TO postgisadmin;
COMMENT ON FUNCTION routing.pgr_fromatob
  IS 'Function returns shortest path route geometry for OL3 web application';

You can test function by running this query against your network:

SELECT * FROM pgr_fromatob('itn_network',325000,725000,350000,750000);

If this works and you get some results back without error then we are good to go on building the web map.

The WMS layer

Set up a new Geoserver workspace, store and layer to hold the information returned by pgRouting. Follow along in http://workshop.pgrouting.org/chapters/geoserver.html#wms-server-with-geoserver and take a note of the name of the new layer you create. You'll use it in the next section.

The web map

From Matt's tutorial we get the basics of the HTML page to hold the OpenLayers3 map. I downloaded the repository and extracted it into a suitable folder on the web server.

<!DOCTYPE html>
<html>
  <head>
    <meta charset="utf-8" />
    <title>pgRouting in OL3</title>
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no, width=device-width">
    <link rel="stylesheet" href="lib/ol3/ol.css" />
    <link rel="stylesheet" href="lib/ol3-popup/src/ol3-popup.css" />
    <link rel="stylesheet" href="ol3.css" />
  </head>
  <body>
    <div id="map"></div>
    <div class="overlay title"><h2>pgRouting in OL3</h2><h3>click start | end points</h3></div>
    <div class="overlay button"><button id="clear">reset</button></div>
    <script src="../common/lib/reqwest.min.js"></script>
    <script src="lib/proj4.js"></script>
    <script src="lib/ol3/ol.js"></script>
    <script src="lib/ol3-popup/src/ol3-popup.js"></script>
    <script src="ol3-complete.js"></script>
  </body>
</html>  

It's the last Javascript file in the code above that we are interested in: ol3-complete.js

// START Matt Walker's code
// Extent of the map in units of the projection (these match our base map)
var extent = [-3276800, -3276800, 3276800, 3276800];

// Fixed resolutions to display the map at (pixels per ground unit (meters when
// the projection is British National Grid))
var resolutions = [1600,800,400,200,100,50,25,10,5,2.5,1,0.5,0.25,0.125,0.0625];

// Define British National Grid Proj4js projection (copied from http://epsg.io/27700.js)
proj4.defs("EPSG:27700","+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs");

// Define an OL3 projection based on the included Proj4js projection
// definition and set it's extent.
var bng = ol.proj.get('EPSG:27700');
bng.setExtent(extent);

// Define a TileGrid to ensure that WMS requests are made for
// tiles at the correct resolutions and tile boundaries
var tileGrid = new ol.tilegrid.TileGrid({
    origin: extent.slice(0, 2),
    resolutions: resolutions
});

var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
            source: new ol.source.TileWMS({
                url: 'http://yourserver/geoserver/base/wms?',
                attributions: [
                    new ol.Attribution({html: 'Angus Council 100023404 &copy; Ordnance Survey'})
                ],
                params: {
                    'LAYERS': 'oslayers',
                    'FORMAT': 'image/png',
                    'TILED': true
                },
                tileGrid: tileGrid
            })
        })
    ],
    view: new ol.View({
        projection: bng,
        resolutions: resolutions,
        center: [350000, 750000],
        zoom: 5
    }),
    controls: ol.control.defaults({
      attributionOptions: {
        collapsible: false
      }
    })
});

// END Matt Walker's code
// START pgRouting workshop code
// pgRouting layer values
// This is the layer you created in Ch. 9 of the pgRouting workshop
var params = {
    'LAYERS': 'routing:itn_routing',
    'FORMAT': 'image/png'
};

// The "start" and "destination" features.
var startPoint = new ol.Feature();
var destPoint = new ol.Feature();

// The vector layer used to display the "start" and "destination" features.
var vectorLayer = new ol.layer.Vector({
  source: new ol.source.Vector({
    features: [startPoint, destPoint]
  })
});
map.addLayer(vectorLayer);

// Note changes to EPSG:27700 from EPSG:4326 and EPSG:3857
// This bit could probably be removed
// A transform function to convert coordinates from EPSG:27700
// to EPSG:27700.
var transform = ol.proj.getTransform('EPSG:27700', 'EPSG:27700');

// Register a map click listener.
map.on('click', function(event) {
  if (startPoint.getGeometry() == null) {
    // First click.
    startPoint.setGeometry(new ol.geom.Point(event.coordinate));
  } else if (destPoint.getGeometry() == null) {
    // Second click.
    destPoint.setGeometry(new ol.geom.Point(event.coordinate));
    // Transform the coordinates from the map projection (EPSG:27700)
    // to the server projection (EPSG:27700).
    var startCoord = transform(startPoint.getGeometry().getCoordinates());
    var destCoord = transform(destPoint.getGeometry().getCoordinates());
    var viewparams = [
      'x1:' + startCoord[0], 'y1:' + startCoord[1],
      'x2:' + destCoord[0], 'y2:' + destCoord[1]
    ];
    params.viewparams = viewparams.join(';');
    result = new ol.layer.Image({
      source: new ol.source.ImageWMS({
        url: 'http://yourserver/geoserver/routing/wms?',    // URL of your 27700 WMS
        params: params
      })
    });
    map.addLayer(result);
  }
});

var clearButton = document.getElementById('clear');
clearButton.addEventListener('click', function(event) {
  // Reset the "start" and "destination" features.
  startPoint.setGeometry(null);
  destPoint.setGeometry(null);
  // Remove the result layer.
  map.removeLayer(result);
}); 
// END pgRouting workshop code

The web map

Once you have finished editing the HTML and JS files, open the HTML file in your browser. This should have got your Geoserver delivering a shortest path route to your OpenLayers3 web map in EPSG:27700. Happy routing!

References

  1. http://workshop.pgrouting.org/ (Ch. 8 - 10)
  2. http://astuntechnology.github.io/osgis-ol3-leaflet/ol3/README.html