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.
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.