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);
?>
The length error stems from this portion of the sql statement:
length(rt.the_geom) AS lengthThis should be changed to justlength. I'm not really sure whatlength(rt.the_geom) AS lengthdoes, 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.