A while back I looked at generating alphashapes using pgRouting and displaying in QGIS. It was a pretty convoluted hack though and my SQL skills need some honing.
After seeing a post by smathermather that generated isochrones using PostgreSQL functions and the pgRouting algorithms I thought there must be a better way to generate my alphashapes. A use case of identifying access to defibrillator sets presented itself.
I looked at the code used by the pgRouting Layer QGIS plugin and stripped out the basic functionality. I used this to generate a temporary table of nodes representing the source and target nodes on the network.
CREATE TEMPORARY TABLE node AS SELECT id, ST_X(geometry) AS x, ST_Y(geometry) AS y, geometry FROM ( SELECT source AS id, ST_Startpoint(geometry) AS geometry FROM itn_network --my network table UNION SELECT target AS id, ST_Startpoint(geometry) AS geometry FROM itn_network --my network table ) AS node;
This table is then used in a function which can be passed any number of input points which it then uses to generate individual alphashapes for. The input points need some preparation as each point needs to be associated to the nearest node on the network. Taken from Anita Graser's work and adapted to my structure:
ALTER TABLE yourpoints ADD COLUMN nearest_node integer; CREATE TEMPORARY TABLE temp AS SELECT a.gid, b.id, min(a.dist) FROM (SELECT yourpoints.gid, min(distance(yourpoints.the_geom, network_vertice_pgr.the_geom)) AS dist FROM yourpoints, network_vertice_pgr GROUP BY yourpoints.gid) AS a, (SELECT yourpoints.gid, network_vertice_pgr.id, distance(yourpoints.the_geom, network_vertice_pgr.the_geom) AS dist FROM yourpoints, network_vertice_pgr) AS b WHERE a.dist = b. dist AND a.gid = b.gid GROUP BY a.gid, b.id; UPDATE yourpoints SET nearest_node = (SELECT id FROM temp WHERE temp.gid = yourpoints.gid);
This assigns each input point a node id from the pgrouting vertices table. The function then uses the nearest_node value for each point as the starting point for the alphashape.
CREATE OR REPLACE FUNCTION public.make_isochrones() RETURNS integer AS $BODY$ DECLARE cur_aed CURSOR FOR SELECT nearest_node FROM yourpoints WHERE nearest_node <> 35819; v_nn cor_aed_sites.nearest_node%TYPE; v_geom geometry; v_sql varchar(120); BEGIN v_sql:='TRUNCATE TABLE public.yourisochrones'; EXECUTE v_sql; -- I truncate the table each time I run the function OPEN cur_aed; LOOP FETCH cur_aed INTO v_nn; EXIT WHEN NOT FOUND; SELECT ST_SetSRID(ST_MakePolygon(ST_AddPoint(foo.openline, ST_StartPoint(foo.openline))),27700) AS geometry --build polygon from line FROM ( SELECT ST_Makeline(points ORDER BY id) AS openline --build line from points FROM ( SELECT row_number() over() AS id, ST_MakePoint(x, y) AS points --get points from alphashape FROM pgr_alphashape(' SELECT * FROM node --get nodes making up alphashape JOIN (SELECT * FROM pgr_drivingDistance('' SELECT gid AS id, source::int4 AS source, target::int4 AS target, cost_time::float8 AS cost, rcost_time::float8 AS reverse_cost FROM itn_network'', '||v_nn||', --loop through each input point 300, --adjust cost as necessary. Mine is in seconds. true, true)) AS dd ON node.id = dd.id1'::text) ) AS a ) AS foo INTO v_geom; --save the result geometry into variable INSERT INTO public.isochrones(geometry) VALUES (v_geom); --insert into table END LOOP; RETURN 1; CLOSE cur_aed; END; $BODY$ LANGUAGE plpgsql VOLATILE COST 100;
I created the table to hold the isochrones beforehand - a simple structure of id and geometry with a primary key and spatial index.
The results of a few hours work distilled down into 15 seconds runtime for 42 input points and a network of 72000 edges.
Approximately 5 minutes travel time to the nearest Automated External Defibrillator [AED]
Approximately 10 minutes travel time to the nearest Automated External Defibrillator [AED]
Provided you have someone at the scene of the cardiac arrest who is proficient in CPR and who can keep going for 20 minutes then if you're in the pink then you may be able to get the AED to save someone's life.
Note1: the list of AED devices shown here is not complete and not all are publicly accessible 24/7.
Note2: Road network is Ordnance Survey's ITN dataset © Crown copyright and database right 2016 10023404.
Note3: There must be a better way to do this...as ever.
Note4: I have worked out a better way of doing it by consolidating more of the process into the function and using some dynamic substitution (just discovered how to do it - amazing!). That's for the next post though as I am on holiday.