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. For reference, the shortest path predicted by Google is as under. 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:  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.

bumped to the homepage by Community♦yesterday

This question has answers that may be good or bad; the system has marked it active so that they can be reviewed.

• try either without ST_MakeLine, or ORDER BY pt.seq. – ThingumaBob Sep 9 at 5:15
• @ThingumaBob It works without ST_MakeLine. I could not figure out where to add the ORDER BY pt.seq. Would be better if I can just return one geometry into several rows that are returned without the st_makeline. Can you explain why this happens? – Chintan Pathak Sep 9 at 5:34
• ST_MakeLine stitches those joined geometries together in, I guess, order of rd.id, which doesn't represent the order of the path; it creates those straight lines between edges that are thus not connected. you could ST_MakeLine(pt.geom ORDER BY pt.seq) to force path order. – ThingumaBob Sep 9 at 6:57
• I had to use rd.geom as pt does not have geom. Using ORDER BY pt.seq gave the same result as before, using pt.edge gives a slightly better result, still with some bad connections. – Chintan Pathak Sep 9 at 7:31
• ah, sure, rd.geom is what I meant. I wrote to hastly, I know where the mistake is, I'll write up an answer with all this... – ThingumaBob Sep 9 at 7:36

There is an official solution to the directionality problem (pointed to by Daniel Kastl). Exercise 8 in Geometry Directionality from the PostGIS tutorial shows how to do it.

Code below in case the link does not work:

WITH
dijkstra AS (
SELECT * FROM pgr_dijkstra(
'SELECT gid AS id, * FROM vehicle_net',
(SELECT id FROM ways_vertices_pgr WHERE osm_id = 252643343),
(SELECT id FROM ways_vertices_pgr WHERE osm_id = 302057309))
),
get_geom AS (
SELECT dijkstra.*, ways.name,
CASE
WHEN dijkstra.node = ways.source THEN the_geom
ELSE ST_Reverse(the_geom)
END AS route_geom
FROM dijkstra JOIN ways ON (edge = gid)
ORDER BY seq)
SELECT seq, name, cost,
degrees(ST_azimuth(ST_StartPoint(route_geom), ST_EndPoint(route_geom))) AS azimuth,
ST_AsText(route_geom),
route_geom
FROM get_geom
ORDER BY seq;

The important construct in the above query is:

CASE
WHEN dijkstra.node = ways.source THEN the_geom
ELSE ST_Reverse(the_geom)
END AS route_geom

i.e. reverse the geometry if the Dijkstra node is not the source and leave as is otherwise. I then used it for ST_MakeLine() and things appear peachy for now.