TechQA.

trouble getting route distance using pgrouting

1.2k views Asked by At

I started with the pgrouting-workshop files, but I'm having trouble getting an accurate route length using multiple methods.

The first method is just adding up the length variables returned from the ajax call to the pgrouting.php script. The length returned from the routing db seems to be obtained from the distance formula: sqrt((x2-x1)^2 + (y2-y1)^s) where the x coordinates are longitude and y coordinates are latitude in epsg:4326. Since degrees of latitude and longitude are different lengths when projected onto the surface of the earth, I'm not really sure what good this value is, but I'm new to this game so...

Since this first method didn't provide an accurate total distance, I decided to just sum up the lengths of each of the segments returned from the ajax call using the haversine formula for the length calculation. However, I compared this length with the "as the crow flies" (ATCF) distance between the starting and ending points and found that this length was less than the ATCF lenth. So I added each of these route segments to a separate vector layer in openlayers and found that these segments did not cover the entire route. A number of segments are missing particularly on curved portions of the route.

So then I thought, well I'll just sum up the gaps by getting the distance between the beginning of one segment and the end of the previous one. However, I found that these segments are not returned in order.

I'm stymied. How can I get an accurate route length using pgrouting? Here are the html and php scripts I used:

routing-final05.html:

<html>
<head>

<title>A Basic GeoExt Page</title>
<script src="/scripts/openlayers/OpenLayers.js" type="text/javascript"></script>
<style type="text/css">
  html, body { height: 100%; }
  body {margin: 0px;}
  #map {
    width: 100%;
    height: 90%;
  }
</style>

<script src="DrawPoints.js" type="text/javascript"></script>

<script src="proj4js/lib/proj4js.js" type="text/javascript"></script>

<script type="text/javascript">

        // global projection objects (uses the proj4js lib)
        var epsg_4326 = new OpenLayers.Projection("EPSG:4326"),
             epsg_900913 = new OpenLayers.Projection("EPSG:900913");

        function pgrouting(layer, method) {

                if (layer.features.length == 2) 
                {

                    route_layer.removeAllFeatures();
                    crow_layer.removeAllFeatures();
                    verify1_layer.removeAllFeatures();
                    verify2_layer.removeAllFeatures();

                    var startpoint = layer.features[0].geometry.clone();
                    startpoint.transform(epsg_900913, epsg_4326);
                    var finalpoint = layer.features[1].geometry.clone();
                    finalpoint.transform(epsg_900913, epsg_4326);

                    var result = {
                            startpoint: startpoint.x + " " + startpoint.y,
                            finalpoint: finalpoint.x + " " + finalpoint.y,
                            method: method,
                    };

//                  layer.removeFeatures(layer.features);

                    var params = OpenLayers.Util.getParameterString(result);

                    OpenLayers.loadURL("./php/pgrouting05.php?" + params,
                                                         '',
                                                         null,
                                                         drawRoute);


             }
        }

        function init(){

                map = new OpenLayers.Map('map', {
                        controls: [
                                new OpenLayers.Control.Navigation(),
                                new OpenLayers.Control.PanZoomBar(),
                                new OpenLayers.Control.LayerSwitcher({'ascending':true}),
                                //new OpenLayers.Control.Permalink(),
                                //new OpenLayers.Control.ScaleLine(),
                                //new OpenLayers.Control.Permalink('permalink'),
                                new OpenLayers.Control.MousePosition(),
                                new OpenLayers.Control.OverviewMap(),
                                new OpenLayers.Control.KeyboardDefaults(),
                                new OpenLayers.Control.Scale({'geodesic':true})
                        ],
                        //numZoomLevels: 6,
                        maxResolution: 156543.0339,
                        units: 'm',
                        projection: new OpenLayers.Projection("EPSG:900913"),
        //                 displayProjection: new OpenLayers.Projection("EPSG:4326"),
                        displayProjection: new OpenLayers.Projection("EPSG:900913"),
                        maxExtent: new OpenLayers.Bounds(-20037508.34, -20037508.34, 20037508.34, 20037508.34)

                });

                // TileLite by default mounts on localhost port 8000 
                var tiles_url = "http://localhost:8000/";

                var tilelite_layer = new OpenLayers.Layer.OSM("Mapnik locally via TileLite", tiles_url + '${z}/${x}/${y}.png');

                map.addLayer(tilelite_layer);

                var lonLat = new OpenLayers.LonLat(-104.99323, 39.74259).transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913"));

                map.setCenter (lonLat, 14);

                // create the layer where the route will be drawn
                route_layer = new OpenLayers.Layer.Vector("route", {
                        styleMap: new OpenLayers.StyleMap(new OpenLayers.Style({
                                strokeColor: "#ff9933",
                                strokeWidth: 3
                        }))
                });

                // create the layer where the "as the crow flies" route will be drawn
                crow_layer = new OpenLayers.Layer.Vector("crow", {
                        styleMap: new OpenLayers.StyleMap(new OpenLayers.Style({
                                strokeColor: "#000000",
                                strokeWidth: 3
                        }))
                });

                // create a verification layer 
                verify1_layer = new OpenLayers.Layer.Vector("verify1", {
                        styleMap: new OpenLayers.StyleMap(new OpenLayers.Style({
                                strokeColor: "#FF0000",
                                strokeWidth: 3
                        }))
                });

                // create a verification layer 
                verify2_layer = new OpenLayers.Layer.Vector("verify2", {
                        styleMap: new OpenLayers.StyleMap(new OpenLayers.Style({
                                strokeColor: "#00FF00",
                                strokeWidth: 3
                        }))
                });

                // create the layer where the start and final points will be drawn
                points_layer = new OpenLayers.Layer.Vector("points");

                // when a new point is added to the layer, call the pgrouting function
                points_layer.events.on({
                        featureadded: function() {
                //                 pgrouting(store, points_layer, method.getValue());
                                pgrouting(points_layer, method);
                        }
                });

                // add the layers to the map
                map.addLayers([points_layer, route_layer, crow_layer, verify1_layer, verify2_layer]);

                // create the control to draw the points (see the DrawPoints.js file)
                var draw_points = new DrawPoints(points_layer);

                // create the control to move the points
                var drag_points = new OpenLayers.Control.DragFeature(points_layer, {
                        autoActivate: true
                });

                // when a point is moved, call the pgrouting function
                drag_points.onComplete = function() {
                //               pgrouting(store, points_layer, method.getValue());
                //               pgrouting(store, points_layer, method);
                                            pgrouting(points_layer, method);
                };

                // add the controls to the map
                map.addControls([draw_points, drag_points]);




                method = "SPS";


        };  // end init()

        function change_method()
        {
//          console.log(document);
            method = document.forms["f1"]["method_select"].value
            pgrouting(points_layer, method);
//          calc_dist();
        }

        function calc_dist()
        {
                console.log("entered calc_dist()")
                console.log(map.layers[2])
                console.log(map.layers[2].features[0])

        }

        function haversine_dist(lon1,lat1,lon2,lat2,unit)
        {

            var radius_km = 6371
            var radius_mi = 3958.75 

            var dlat = (lat2-lat1)*Math.PI/180.0;
            var dlon = (lon2-lon1)*Math.PI/180.0;
            var a = Math.sin(dlat/2) * Math.sin(dlat/2) + Math.cos(lat1*Math.PI/180.0) * Math.cos(lat2*Math.PI/180.0) * Math.sin(dlon/2) * Math.sin(dlon/2)
            var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a))
            if (unit == 'mi')
            {
                var d = radius_mi * c
            }
            else if (unit == 'km')
            {
                var d = radius_km * c
            }

            return d;

        }

        function drawRoute(response)
        {
            // interpret geojson data coming in from ajax response
            var geojson_format900913 = new OpenLayers.Format.GeoJSON({
                                                        internalProjection: epsg_900913,
                                                        externalProjection: epsg_4326
                                                });

            var geojson_format4326 = new OpenLayers.Format.GeoJSON();

            var geom4326 = geojson_format4326.read(response.responseText)
            var geom900913 = geojson_format900913.read(response.responseText)

//          console.log(geom900913);
//          console.log(geom4326);

            // add features
            route_layer.addFeatures(geom900913);

            // calculate distances between current segment start and previous segment end
            var lon1,lat1,lon2,lat2,seg_len;
            var pointList = [];
            var p1,p2;
            var dist2 = 0;
            for (i=0;i<geom900913.length;i++)
            {
                    if (i == 0)
                    {

                            var startpoint = points_layer.features[0].geometry.clone();
//                          var startpoint = points_layer.features[1].geometry.clone();

                            lon1 = startpoint.x;
                            lat1 = startpoint.y;
                            lon2 = geom900913[i].geometry.components[0].x;
                            lat2 = geom900913[i].geometry.components[0].y;

                    }
                    else
                    {
                            lon1 = geom900913[i-1].geometry.components[1].x
                            lat1 = geom900913[i-1].geometry.components[1].y
                            lon2 = geom900913[i].geometry.components[0].x
                            lat2 = geom900913[i].geometry.components[0].y


                    }

                    seg_len = haversine_dist(lon1,lat1,lon2,lat2,'mi');

                    if ( seg_len > 0 )
                    {
                            dist2 += seg_len;   

                            var pointList = [];
                            var p1 = new OpenLayers.Geometry.Point(lon1,lat1);
                            var p2 = new OpenLayers.Geometry.Point(lon2,lat2);

                            pointList.push(p1);
                            pointList.push(p2);

                            var route_seg = new OpenLayers.Feature.Vector(
                                        new OpenLayers.Geometry.LineString(pointList),null,null);

                            verify2_layer.addFeatures([route_seg]);

                    }

            }
//          verify2_layer.addFeatures(geom4326.transform(epsg_4326, epsg_900913));

            var sp900913 = points_layer.features[0].geometry.clone();
            var sp4326 = points_layer.features[0].geometry.clone();
            sp4326.transform(epsg_900913, epsg_4326);
            var fp900913 = points_layer.features[1].geometry.clone();
            var fp4326 = points_layer.features[1].geometry.clone();
            fp4326.transform(epsg_900913, epsg_4326);

//          console.log(sp4326);
//          console.log(fp4326);

            var lengths = [];
            var tot_length = 0.0;
//          console.log(tot_length)
            var lon1,lat1,lon2,lat2;
            for (i=0;i<geom4326.length;i++)
            {
//              lengths.push(geom4326[i].data.length);
                lon1 = geom4326[i].geometry.components[0].x
                lat1 = geom4326[i].geometry.components[0].y
                lon2 = geom4326[i].geometry.components[1].x
                lat2 = geom4326[i].geometry.components[1].y

                tot_length += parseFloat(haversine_dist(lon1,lat1,lon2,lat2,'mi'))
//              console.log(tot_length)

                var pointList = [];
                var p1 = new OpenLayers.Geometry.Point(lon1,lat1).transform(epsg_4326, epsg_900913);
                var p2 = new OpenLayers.Geometry.Point(lon2,lat2).transform(epsg_4326, epsg_900913);
                pointList.push(p1);
                pointList.push(p2);
    //          console.log(pointList)

                var route_seg = new OpenLayers.Feature.Vector(
                      new OpenLayers.Geometry.LineString(pointList),null,null);

    //          console.log(crow_line)

                verify1_layer.addFeatures([route_seg]);

            }

//          console.log(lengths);

//          console.log("sum of the parts: " + tot_length);

            // as the crow flies
            lon1 = sp4326.x
            lat1 = sp4326.y
            lon2 = fp4326.x
            lat2 = fp4326.y

            var line_length_mi = haversine_dist(lon1,lat1,lon2,lat2,'mi')
            var line_length_km = haversine_dist(lon1,lat1,lon2,lat2,'km')

//          console.log("straight line dist: " + line_length);

            document.getElementById("dist").innerHTML = "Distance: " + Math.round(tot_length*1000)/1000 + " miles"
            document.getElementById("straight_dist").innerHTML = "As the crow flies: " + Math.round(line_length_mi*1000)/1000 + " miles, " + Math.round(line_length_km*1000)/1000 + " km"

            var pointList = [];
            pointList.push(new OpenLayers.Geometry.Point(sp900913.x,sp900913.y));
            pointList.push(new OpenLayers.Geometry.Point(fp900913.x,fp900913.y));
//          console.log(pointList)

            var crow_line = new OpenLayers.Feature.Vector(
                new OpenLayers.Geometry.LineString(pointList),null,null);

//          console.log(crow_line)

            crow_layer.addFeatures([crow_line]);



        }

</script>
</head>
<body onload="init()">
<form name="f1">
    <table width="100%">
        <tr>
            <td align="center">
                <select name="method_select" onchange="change_method();">
                    <option value="SPD">Shortest Path Dijkstra</option>
                    <option value="SPA">Shortest Path A*</option>
                    <option value="SPS">Shortest Path Shooting*</option>
                </select>
            </td>
            <td align="center">
                <p id="dist"></p>
            </td>
            <td align="center">
                <p id="straight_dist"></p>
            </td>
        </tr>
    </table>
</form>
<div id="map" ></div>
<div id="gxmap"></div>
<div id="method"></div>


</body>
</html>

pgrouting.php (I didn't change anything here, this is straight from the workshop files)

<?php

        // Database connection settings
        define("PG_DB"  , "routing");
        define("PG_HOST", "localhost"); 
        define("PG_USER", "postgres");
        define("PG_PORT", "5432"); 
        define("TABLE",   "ways");

        //echo($_REQUEST)

        // Retrieve start point
        $start = preg_split('/ /',$_REQUEST['startpoint']);
        $startPoint = array($start[0], $start[1]);

        // Retrieve end point
        //$end = split(' ',$_REQUEST['finalpoint']);
        $end = preg_split('/ /',$_REQUEST['finalpoint']);
        $endPoint = array($end[0], $end[1]);




        //echo($end[0]);

?>

<?php

   // Find the nearest edge
   $startEdge = findNearestEdge($startPoint);
   $endEdge   = findNearestEdge($endPoint);

   // FUNCTION findNearestEdge
   function findNearestEdge($lonlat) {

      // Connect to database
      $con = pg_connect("dbname=".PG_DB." host=".PG_HOST." user=".PG_USER);

      $sql = "SELECT gid, source, target, the_geom, 
                  distance(the_geom, GeometryFromText(
                       'POINT(".$lonlat[0]." ".$lonlat[1].")', 4326)) AS dist 
                 FROM ".TABLE."  
                 WHERE the_geom && setsrid(
                       'BOX3D(".($lonlat[0]-0.1)." 
                              ".($lonlat[1]-0.1).", 
                              ".($lonlat[0]+0.1)." 
                              ".($lonlat[1]+0.1).")'::box3d, 4326) 
                 ORDER BY dist LIMIT 1";

      $query = pg_query($con,$sql);  

      $edge['gid']      = pg_fetch_result($query, 0, 0);  
      $edge['source']   = pg_fetch_result($query, 0, 1);  
      $edge['target']   = pg_fetch_result($query, 0, 2);  
      $edge['the_geom'] = pg_fetch_result($query, 0, 3);  

      // Close database connection
      pg_close($con);

      return $edge;
   }

?>

<?php

   // Select the routing algorithm
   switch($_REQUEST['method']) {

      case 'SPD' : // Shortest Path Dijkstra 

        $sql = "SELECT rt.gid, ST_AsGeoJSON(rt.the_geom) AS geojson, 
                     length(rt.the_geom) AS length, ".TABLE.".gid 
                  FROM ".TABLE.", 
                      (SELECT gid, the_geom 
                          FROM dijkstra_sp_delta(
                              '".TABLE."',
                              ".$startEdge['source'].",
                              ".$endEdge['target'].",
                              0.1)
                       ) as rt 
                  WHERE ".TABLE.".gid=rt.gid;";
        break;

      case 'SPA' : // Shortest Path A* 

        $sql = "SELECT rt.gid, ST_AsGeoJSON(rt.the_geom) AS geojson, 
                       length(rt.the_geom) AS length, ".TABLE.".gid 
                    FROM ".TABLE.", 
                        (SELECT gid, the_geom 
                            FROM astar_sp_delta(
                                '".TABLE."',
                                ".$startEdge['source'].",
                                ".$endEdge['target'].",
                                0.1)
                         ) as rt 
                    WHERE ".TABLE.".gid=rt.gid;";  
        break;

      case 'SPS' : // Shortest Path Shooting*

        $sql = "SELECT rt.gid, ST_AsGeoJSON(rt.the_geom) AS geojson, 
                       length(rt.the_geom) AS length, ".TABLE.".gid 
                    FROM ".TABLE.", 
                        (SELECT gid, the_geom 
                            FROM shootingstar_sp(
                                '".TABLE."',
                                ".$startEdge['gid'].",
                                ".$endEdge['gid'].",
                                0.1, 'length', true, true)
                         ) as rt 
                    WHERE ".TABLE.".gid=rt.gid;";
        break;   

   } // close switch

   // Connect to database
   $dbcon = pg_connect("dbname=".PG_DB." host=".PG_HOST." user=".PG_USER);

   // Perform database query
   $query = pg_query($dbcon,$sql); 

?>

<?php

   // Return route as GeoJSON
   $geojson = array(
      'type'      => 'FeatureCollection',
      'features'  => array()
   ); 

   // Add edges to GeoJSON array
   while($edge=pg_fetch_assoc($query)) {  

      $feature = array(
         'type' => 'Feature',
         'geometry' => json_decode($edge['geojson'], true),
         'crs' => array(
            'type' => 'EPSG',
            'properties' => array('code' => '4326')
         ),
         'properties' => array(
            //'osm_id' => $edge['osm_id'],
            'length' => $edge['length']
            //'x1' => $edge['x1']
         )
      );

      // Add feature array to feature collection array
      array_push($geojson['features'], $feature);
   }


   // Close database connection
   pg_close($dbcon);

   // Return routing result
   header('Content-type: application/json',true);
   echo json_encode($geojson);

?>
php ajax openlayers mapnik pgrouting
Original Q&A
2

There are 2 answers

0
nomadicME On BEST ANSWER

The length error stems from this portion of the sql statement: length(rt.the_geom) AS length This should be changed to just length. I'm not really sure what length(rt.the_geom) AS length does, but as I stated in my original question this seems to use the distance formula.

There are still gaps in the route if you use the route segments returned by geojson, but the total route distance (using the modified length) is now accurate to within 8 meters when compared to the "as the crow flies distance" described in the original question (along a straight route). This must mean that the x and y values in geometry for some segments is wrong, thus leading to gaps.

1
mattesCZ On

I didn't checked up pgRouting workshop, but the "length returned from the routing db" is probably cost, because cost is returned from pgRouting SQL call, see the docs. Sum of cost values equals to the real length only if you use real length as cost in the SQL parameter of pgRouting functions.

Assuming you have correct lengths precomputed in your database of ways and use length as cost, you can get accurate length by using your first method.

If length is not equal to cost, you can use edge_id column of your result and somehow get sum of lengths of all ways that have id in edge_id. Something like this (not tested):

SELECT sum(length) FROM my_ways WHERE id IN (
   SELECT edge_id FROM shortest_path_astar("SQL", source, target, false, false)
) as edge_ids

P. S. I don't know which app you use to create database of ways, but I'm sure osm2po computes correct length of ways.

Related Questions in PHP

Related Questions in AJAX

Related Questions in OPENLAYERS

Related Questions in MAPNIK

Related Questions in PGROUTING

Popular Questions

Popular Tags

javascript python java c# php android html jquery c++ css ios sql mysql r reactjs node.js arrays c asp.net json python-3.x ruby-on-rails .net sql-server swift django angular objective-c pandas excel

Trending Questions