Scotland, distorted by drive time using FOSS tools

This map shows Scotland, distorted by driving time from Edinburgh.

The map is centred on Edinburgh, and the bearing (azimuth) to each point on the road network has been kept the same. The distance has been changed to reflect the drive-time from Edinburgh, following Motorways, A roads and B roads.

Note on interpreting : you can’t compare distances between two arbitrary points with this sort of chart. But you can compare distances to Edinburgh - e.g. by putting a finger on Edinburgh, and sweeping another finger around like a pair of compasses


The road data was OS Open Data (crown copyright an database rights 2015).

The data was imported into a PostgreSQL database with the postgis and pgrouting extensions added, using the shp2pgsql tool. pgRouting was used to create the network topology, and costs were weighted according to road types and lengths (cars go faster on Motorways than B-roads, for example)

Then, pg_drivingdistance was used to work out the time needed to reach each point in the network. The starting node number (13985) was discovered using the QGIS PgRouting Layer plugin.

CREATE TABLE drivecostsbytime AS (
    SELECT seq, id1 AS node, cost
        FROM pgr_drivingDistance(
                'SELECT gid as id, source, target, cost_time as cost FROM justscotland',
                13985, 6000000.0, false, false
        )
        ORDER BY cost ASC
);

Next, add in the geometry...

CREATE TABLE drivetimes2bytime AS (
    SELECT
        dt.node, dt.cost, verts.the_geom as geom
    FROM drivecostsbytime as dt, justscotland_vertices_pgr as verts
    WHERE dt.node = verts.id
);

Using a postgis query, I followed these steps for each point

  • draw a line from Edinburgh to the point using ST_MakeLine
  • Translated all points so that Edinburgh is at (0,0) with ST_Translate
  • scaled it down to unit length (length 1, but pointing in the correct direction) using ST_Scale
  • scaled it back up using the cost of the node from the drivetime output
  • took the endpoint of this line with ST_EndPoint
  • translated every point so that Edinburgh is in the middle of the OSGB grid
CREATE TABLE edi3 AS (
    SELECT  dt.node as id, 

    ST_Translate(
        -- translate to false origin 
        ST_EndPoint(
            -- scale to length of cost
            ST_Scale(
                ST_Scale(
                    -- move center point to origin
                    ST_Translate(
                        ST_MakeLine(verts.the_geom,dt.geom),
                        - ST_X(verts.the_geom),
                        - ST_Y(verts.the_geom)
                    ),
                    -- scale down to unit length
                    1.0/(ST_Length(ST_MakeLine(verts.the_geom,dt.geom))+.00001), -- avoids divide by zero
                    1.0/(ST_Length(ST_MakeLine(verts.the_geom,dt.geom))+.00001)         
                ),
                dt.cost/1000.0,
                dt.cost/1000.0
            )
        ),
        -- x and y coordinates of false origin
        500000

    ) as geom
FROM    
        drivetimes2bytime as dt,
        justscotland_vertices_pgr as verts
wHERE   
        verts.id = 13985
ORDER BY cost ASC
)

These steps meant I now had each point in its new place on the polar grid.

A final join between the original table and the warped point locations meant I could recreate the line segments with ST_Makeline

select
    ST_AsEWKT(ST_MakeLine(nodes1.geom,nodes2.geom))
from 
    EDI3 as nodes1,
    EDI3 as nodes2,
    justscotland as segs
where 
    nodes1.id = segs.source and nodes2.id = segs.target

For debugging, I found it useful to use ST_AsEWKT() on the geometries and export the results from pgAdminIII as CSVs, then import into QGIS as delimited layers. eWKT is useful as it can be read - you can see if you've made mistakes.