trouble getting route distance using pgrouting
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 layer 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);
?>
No comments:
Post a Comment