Wednesday 30 November 2016

pgrouting - Using pgr_drivingDistance for isochrones?


I have tried to calculate Isochrones in 500m walking distance along paths around planned subway stations with pgrouting.


For the setup I used an OSM-road.shp clipped to the borders of Vienna and a point shp. with the 14 stations which I want to analyse.



I worked my way through several tutorial since I am new to pgrouting for setting up a street-network and finding the closest nodes from the stations.


So far so good…but when I to calculate the isochrones with pgr_DrivingDistance the output shows very few results.


Either I set up the network wrong or I have been using the DrivingDistance function wrong. It seems like the function is not able to use every path.


How do I improve the query?


PostgreSQL 9.5.11 Pgrouting 2.4.1


CREATE OR REPLACE VIEW road_ext AS
SELECT *, ST_StartPoint(geom), ST_EndPoint(geom)
FROM road;

CREATE TABLE node AS

SELECT row_number() OVER (ORDER by foo.p)::integer AS id,
foo.p AS geom
FROM (
SELECT DISTINCT road_ext.ST_StartPoint AS p FROM road_ext
UNION
SELECT DISTINCT road_ext.ST_EndPoint AS p FROM road_ext
) foo
GROUP BY foo.p;

I had to change the geometry type since OSM provided multistring



ALTER TABLE road
ALTER COLUMN geom TYPE geometry(LineString, 31256)
USING ST_GeometryN(geom, 1);

Setup network


CREATE TABLE network AS 
SELECT a.*, b.id as start_id, c.id as end_id
FROM road_ext AS a
JOIN node AS b ON a.ST_Startpoint = b.geom
JOIN node AS c ON a.ST_EndPoint = c.geom;



ALTER TABLE network add column len_km double precision

ALTER TABLE network add column time_walking double precision

UPDATE network set len_km = (ST_Length(geom))/1000

UPDATE network set time_walking = ((60/5::double precision)*len_km)


Setup the point shp


ALTER TABLE ubahnpointsmgi
ADD Column nearest_node integer


--Node in clostest distance is enough accuracy
CREATE TABLE temp AS
SELECT
DISTINCT ON (a.gid) a.gid, b.id AS start_id,
ST_Distance(a.geom, b.geom) AS dist

FROM ubahnpointsmgi a, node b
ORDER BY a.gid, ST_Distance(a.geom, b.geom), b.id

UPDATE ubahnpointsmgi
SET nearest_node =
(SELECT start_id
FROM temp
WHERE temp.gid = ubahnpointsmgi.gid)

Calculating Line geometry



CREATE TABLE ISO_500 AS
SELECT seq, cost, agg_cost, node, edge, l.geom FROM pgr_drivingDistance(
'SELECT gid as id, start_id as source, end_id as target, len_km as
cost FROM network',
41217,0.5, false) AS foo JOIN network l ON l.gid = foo.edge

Calculating Point geometry


CREATE TABLE ISO_500 AS
SELECT seq, cost, node, edge, l.geom FROM pgr_drivingDistance(
'SELECT gid as id, start_id as source, end_id as target, len_km as

cost FROM network',
41217,0.5, false) AS foo LEFT JOIN node l ON l.id = foo.node

Output Table:enter image description here
Output QGIS:enter image description here



Answer



The problem with the calculation was the network based on the osm.shp. After a closer look on the lines I realised that a lot of the lines provided from osm a very long segments. Therefore fewer nodes were created and fewer connection between the paths were available.


The solution for my problem was to use a road network which uses shorter line segments with way more connections between them. Luckily such a shp is provided by the open government data for Vienna. If such data would have been not available, I would have to break up the long lines of the OSM data in order to generate shorter segments.


In the image below the results for the same point can be seen. Sadly there are not as many paths in general in the new network in comparison to the OSM shp., but this can be handled with a simple merge of the missing paths.


enter image description here



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...