Friday, 29 November 2019

postgis - pgr_dijkstra gives wacky routes sometimes with undirected graph


I am using the function below to find the shortest path between several origin-destination pairs, and it seems to work correctly most of the time.


CREATE OR REPLACE FUNCTION public.sp_od(
orig integer,
dest integer)
RETURNS TABLE(shortest_path geometry)
LANGUAGE 'sql'


AS $BODY$

SELECT st_makeline(geom) as shortest_path
FROM pgr_dijkstra(
'SELECT id, source, target, st_length(geom, true) as cost FROM public."WA_roads"',
(SELECT source FROM public."WA_roads"
ORDER BY ST_StartPoint(geom) <->
(select ST_SetSRID(ST_MakePoint(CAST(ocentx as double precision), CAST(ocenty as double precision)), 4326) from all_trips_non_zero where origin = orig LIMIT 1) ASC
LIMIT 1),

(SELECT source FROM public."WA_roads"
ORDER BY ST_StartPoint(geom) <->
(select ST_SetSRID(ST_MakePoint(CAST(dcentx as double precision), CAST(dcenty as double precision)), 4326) from all_trips_non_zero where destination = dest LIMIT 1) ASC
LIMIT 1), directed := false
) as pt
JOIN public."WA_roads" rd ON pt.edge = rd.id;

$BODY$;

The above function takes an origin and destination zip code and using the corresponding lat, long for the zip code from 'all_trips_non_zero' table, finds the nearest "source" node in our road network closest to the origin and destination points to use for shortest path calculation.



One of the problematic path is shown in the image below. Origin in red and destination is yellow. Brown is the network and blue is the shortest path from the above function.


enter image description here


For reference, the shortest path predicted by Google is as under.


enter image description here


I cannot explain some of the straight lines, not on the network, in the pgr_dijkstra path above. These are somehow part of the geometry and I get a total path length of around 360 miles, which is more than twice the distance predicted by Google. When I use directed:= true, I get no path for this OD pair. The "WA_roads" shapefile and "all_trips_non_zero" file are here. I further transformed the shapefile to EPSG:4326 as described here. The pgr_analyzegraph output is below:


NOTICE:  pgr_analyzeGraph('WA_roads',1e-06,'geom','id','source','target','true')
NOTICE: Performing checks, please wait ...
NOTICE: Analyzing for dead ends. Please wait...
NOTICE: Analyzing for gaps. Please wait...
NOTICE: Analyzing for isolated edges. Please wait...

NOTICE: Analyzing for ring geometries. Please wait...
NOTICE: Analyzing for intersections. Please wait...
NOTICE: ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE: Isolated segments: 0
NOTICE: Dead ends: 214
NOTICE: Potential gaps found near dead ends: 0
NOTICE: Intersections detected: 1
NOTICE: Ring geometries: 2

Successfully run. Total query runtime: 1 secs 204 msec.

1 rows affected.


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