We’ve got some quick and dirty pgRouting-based code up on github. I say quick and dirty because it directly references the table names in both of the functions. I hope to fix this in the future.
The objective with this code is to input a point, use a nearest neighbor search to find the nearest intersection, and from that calculate a drive time alpha shape.
First, the nearest neighbor search:
CREATE OR REPLACE FUNCTION pgr.nn_search (geom geometry) RETURNS int8 AS $$ SELECT id FROM pgr.roads_500_vertices_pgr AS r ORDER BY geom <#> r.the_geom LIMIT 1; $$ LANGUAGE SQL VOLATILE;
This function takes one argument– pass it a point geometry, and it will do a knn search for the nearest point. Since we are leveraging pgRouting, it simply returns the id of the nearest intersection point. We wrote this as a function in order to run it, e.g. on all the points in a table. As I stated earlier, it directly references a table name, which is a little hackerish, but we’ll patch that faux pax later.
Now that we have the ID of the point in question, we can do a driving distance calculation, wrap that up in an alpha shape, and return our polygon. Again, we write this as a function:
CREATE OR REPLACE FUNCTION pgr.alpha_shape (id integer, minutes integer) RETURNS geometry AS $$ WITH alphashape AS(SELECT pgr_alphaShape('WITH DD AS ( SELECT seq, id1 AS node, cost FROM pgr_drivingDistance(''SELECT id, source, target, cost FROM pgr.roads_500'',' || id || ', ' || minutes || ', false, false)), dd_points AS ( SELECT id_ AS id, x, y FROM pgr.roads_500_vertices_pgr v, DD d WHERE v.id = d.node) SELECT * FROM dd_points')), alphapoints AS ( SELECT ST_Makepoint((pgr_alphashape).x, (pgr_alphashape).y) FROM alphashape), alphaline AS ( SELECT ST_MakeLine(ST_MakePoint) FROM alphapoints) SELECT ST_MakePolygon(ST_AddPoint(ST_MakeLine, ST_StartPoint(ST_MakeLine))) AS the_geom FROM alphaline $$ LANGUAGE SQL VOLATILE;
Finally, we’ll use these functions in conjunction with a set of park feature access points to map our our 5-minute drive time. Overlaps in 5-minute zones we’ve rendered as a brighter green in the image above.
CREATE TABLE pgr.alpha_test AS WITH dest_ids AS ( SELECT pgr.nn_search(the_geom) AS id FROM pgr.dest_pts ) SELECT pgr.alpha_shape(id::int, 5)::geometry, id FROM dest_ids;
Oh, and I should point out that for the driving distance calculation, we have pre-calculated the costs of the roads based on the speed limit, so cost is really travel time in minutes.