Monday, August 27, 2012


For every segment, turn, turn, turn, turn, ...

So, it's been a few weeks of  mildly complicated stuff but not a lot that really seemed blog worthy. But, these last three days have been spent on TURNING DIRECTIONS for the trails routing.

It wasn't good enough that  our directions read: Start on Blackberry Loop Trail, 0.25 miles, Red Wing Hiking Path, 1.84 miles, Soda Springs Maintenance Road, 7.25 miles. Turning directions are really nice: TURN RIGHT onto Red Wing, SLIGHT LEFT onto Soda Springs.

Now, pgRouting's dijkstra_sp() et al return a bunch of line segments and the segments' GIDs. And you can do a simple join to fetch the segments' other info (its name, whether it allows bicycles) in one go. But, turning directions are something else entirely, and must be computed based on the relative azimuths of the two segments. (the WHAT? tell you later)

The basic query to fetch segments
ST_AsText(ST_Transform(route.the_geom,4326)) AS wkt,
ST_Length(route.the_geom) AS length,, hiking_trails.duration
(SELECT gid, the_geom FROM astar_sp_delta_directed('hiking_trails',?,?,2640,true,true)) AS route
This fetches the route, and does an implicit join to get the segments' name. Note that we also fetched the Well Known Text (WKT) for each segment -- that comes in handy later when we have to do geometry stuff in PHP.

Iterate over the segments and build directions

The strategy here is that we iterate over the segments in PHP, and keep on tallying up the distance and the time. When we detect that the new segment is no longer on the same road as the previous segment, we call this a transition and we add a new step to the directions. As part of this step, we calculate text directions including whether you turn right or left.

Obvious question: Why didn't we aggregate the trail segments by name within the database? Answer: the trail names are really a mess, and it's really only feasible using more complex expressions. In the code below, the Trailpiece::trailContainsSameName() function encapsulates this -- it returns true to indicate that the current trail segment should be considered a continuation of the previous step (e.g. we are still on Blackberry Loop, even though the text name may imply otherwise) and a false indicates that this segment is not a continuation but should be considered a new trail and a new step.

The turning directions are the largest part of the code here. The strategy is as follows: if we are at a transition, where this segment is a new step from the previous segment, calculate the azimuth of both this segment and of the prior segment. (azimuth: the compass heading, e.g. 45 is northeast, 210 is south-southwest) The difference between these azimuths is the direction you would turn: if the previous segment was on a heading 90 (east) and the current segment is on a heading 180 (south) then it's a +90 difference or a right turn. Easy, huh?
    $steps      = array();
    $stepnumber = 0;
    $current_step_distance = 0;
    $current_step_duration = 0;
    $current_step_name     = $segments[0]->name;
    for ($i=0; $i<$howmanysegments; $i++) {
        $segment       = $segments[$i];
        $is_transition = ! Trailpiece::trailContainsSameName($current_step_name,$segment->name);
        $is_end        = $i == sizeof($segments)-1;

        // if this is a transition from some other step, zero out the counters for this step and append the step
        //but half of the code is just to determine Right or Left
        if ($is_transition) {
            // phase 1: fetch the first and last vertices for the previous and current segment, using regular expressions
            $thisvert  = preg_split('/,\s*/', preg_replace( array('/^[\w\s]+\(+/', '/[\)]+$/') , array('',''), $segment->wkt ) );
            $thisvert1 = explode(' ', $thisvert[0]);
            $thisvert2 = explode(' ', $thisvert[ sizeof($thisvert)-1]);
            $prevvert = preg_split('/,\s*/', preg_replace( array('/^[\w\s]+\(+/', '/[\)]+$/') , array('',''), $segments[$i-1]->wkt ) );
            $prevvert1 = explode(' ', $prevvert[0] );
            $prevvert2 = explode(' ', $prevvert[ sizeof($prevvert)-1 ] );
            $thislon1 = $thisvert1[0]; $thislat1 = $thisvert1[1];
            $thislon2 = $thisvert2[0]; $thislat2 = $thisvert2[1];
            $prevlon1 = $prevvert1[0]; $prevlat1 = $prevvert1[1];
            $prevlon2 = $prevvert2[0]; $prevlat2 = $prevvert2[1];

            // phase 2: either/both of the line segments may need to be flipped, depending on the distance, since the endpoints may be the two touching ends, the two far ends, or any combination
            // the vertices as listed above, may give the azimuth from the segment's end to its start, backwards!
            // strategy: find which combination of endpoints is closest together, and that would be the two touching endpoints
            // remember, "1" indicates the start of a segment and "2" indicates the end of a segment, so $dx12 means the distance from previous seg start to current seg end
            // two segments should meet with previous2 touching current1 ($dx21 is smallest), for the previous to END where the current one STARTS
            // if this is not the case, then one or both of the segments needs to have its vertices swapped
            $dx11 = ($thislon1 - $prevlon1) * ($thislon1 - $prevlon1) + ($thislat1 - $prevlat1) * ($thislat1 - $prevlat1); // distance (squared) between $thisvert1 and $prevvert1
            $dx22 = ($thislon2 - $prevlon2) * ($thislon2 - $prevlon2) + ($thislat2 - $prevlat2) * ($thislat2 - $prevlat2); // distance (squared) between $thisvert2 and $prevvert2
            $dx12 = ($thislon1 - $prevlon2) * ($thislon1 - $prevlon2) + ($thislat1 - $prevlat2) * ($thislat1 - $prevlat2); // distance (squared) between $thisvert1 and $prevvert2
            $dx21 = ($thislon2 - $prevlon1) * ($thislon2 - $prevlon1) + ($thislat2 - $prevlat1) * ($thislat2 - $prevlat1); // distance (squared) between $thisvert2 and $prevvert1
            $whichdx = min(array($dx11,$dx22,$dx12,$dx21));
            switch ($whichdx) {
                case $dx11:
                    // previous segment's start meets current segment start; flip the previous segment
                    list($prevvert1,$prevvert2) = array($prevvert2,$prevvert1);
                    $prevlon1 = $prevvert1[0]; $prevlat1 = $prevvert1[1];
                    $prevlon2 = $prevvert2[0]; $prevlat2 = $prevvert2[1];
                case $dx12:
                    // segments are end-to-end and both need to be flipped
                    list($thisvert1,$thisvert2) = array($thisvert2,$thisvert1);
                    $thislon1 = $thisvert1[0]; $thislat1 = $thisvert1[1];
                    $thislon2 = $thisvert2[0]; $thislat2 = $thisvert2[1];
                    list($prevvert1,$prevvert2) = array($prevvert2,$prevvert1);
                    $prevlon1 = $prevvert1[0]; $prevlat1 = $prevvert1[1];
                    $prevlon2 = $prevvert2[0]; $prevlat2 = $prevvert2[1];
                case $dx22:
                    // current segment end meets previous segment end, flip the current segment
                    list($thisvert1,$thisvert2) = array($thisvert2,$thisvert1);
                    $thislon1 = $thisvert1[0]; $thislat1 = $thisvert1[1];
                    $thislon2 = $thisvert2[0]; $thislat2 = $thisvert2[1];
                case $dx21:
                    // current start is previous end, already fine

            // phase 3: find the azimuth of each, and thus the angle between them
            $thisaz = (180 + rad2deg(atan2(sin(deg2rad($thislon2) - deg2rad($thislon1)) * cos(deg2rad($thislat2)), cos(deg2rad($thislat1)) * sin(deg2rad($thislat2)) - sin(deg2rad($thislat1)) * cos(deg2rad($thislat2)) * cos(deg2rad($thislon2) - deg2rad($thislon1)))) ) % 360;
            $prevaz = (180 + rad2deg(atan2(sin(deg2rad($prevlon2) - deg2rad($prevlon1)) * cos(deg2rad($prevlat2)), cos(deg2rad($prevlat1)) * sin(deg2rad($prevlat2)) - sin(deg2rad($prevlat1)) * cos(deg2rad($prevlat2)) * cos(deg2rad($prevlon2) - deg2rad($prevlon1)))) ) % 360;
            $angle = round($thisaz - $prevaz);
            if ($angle > 180)  $angle = $angle - 360;
            if ($angle < -180) $angle = $angle + 360;
            //printf("%s x %s = %d x %d = %d<br/>\n", $current_step_name, $segment->name, $prevaz, $thisaz, $angle );

            // phase 4: assign a direction word based on that angle
            $turnword = "Turn onto";
            if      ($angle >= -30 and $angle <= 30)   $turnword = "Continue on";
            else if ($angle >= 31  and $angle <= 60)   $turnword = "Take a slight right onto";
            else if ($angle >= 61  and $angle <= 100)  $turnword = "Take a right onto";
            else if ($angle >= 101)                    $turnword = "Take a sharp right onto";
            else if ($angle <= -30 and $angle >= -60)  $turnword = "Take a slight left onto";
            else if ($angle <= -61 and $angle >= -100) $turnword = "Take a left onto";
            else if ($angle <= -101)                   $turnword = "Take a sharp left onto";

            // add the step to the list
            $step = array(
                'stepnumber' => ++$stepnumber,
                'turnword' => $turnword, 'text' => $segment->name,
                'distance' => $current_step_distance, 'duration' => $current_step_duration
            $steps[] = $step;

            // reset the counters for this next step
            $current_step_distance = 0;
            $current_step_duration = 0;
            $current_step_name     = $segment->name;

        // increment the length & duration of the current step, even if that step was just now reset because of a transition
        $current_step_distance += $segment->length;
        $current_step_duration += $segment->seconds;

        // and lastly, if this is the end segment, add the Arrival step so we can indicate the length of travel on this last step
        if ($is_end) {
            $step = array(
                'stepnumber' => ++$stepnumber,
                'turnword' => "Arrive at", 'text' => "your destination",
                'distance' => $current_step_distance, 'duration' => $current_step_duration
            $steps[] = $step;

    // prepend the Start At step, to indicate the name of the street where we start
    array_unshift($steps, array(
        'stepnumber' => null,
        'turnword' => "Start on", 'text' => $segments[0]->name,
        'distance' => null, 'duration' => null
And there you have it: iterating over segments and not only deciding when to start a new step, but also figuring out which direction you would be turning.

Discussion: Which vertices to compare

Phase 1 fetches the first and last vertex of the current linestring and of the previous linestring, and uses these to determine the heading of each path. Technically, this should be the last vertex and the next-to-last vertex, like this:

            $thisvert1 = explode(' ', $thisvert[ sizeof($thisvert)-2 ]);
            $thisvert2 = explode(' ', $thisvert[ sizeof($thisvert)-1 ]);
            $prevvert1 = explode(' ', $prevvert[ sizeof($prevvert)-2 ]);
            $prevvert2 = explode(' ', $prevvert[ sizeof($prevvert)-1 ]);

However, in this case the individual vertices are very high resolution, and the average length of such segments is only 1-3 feet. As such, a single-pixel hand-cramp while generating the data could show the trail as facing northeast when in fact only the last 4 feet are northeast and the rest is due east.

Our use of the first and last does give some surprises: a very, very convoluted segment can have its start and end vertices indicate a southwest azimuth when in fact the last 20 feet of trail leads north. But, using the last segment or two gave such wonky results regularly and I eventually resigned myself that "sometimes slightly off" is better than "usually wrong"

Discussion: segment flipping

When I first read other discussion of generating turning directions, I didn't get the part about flipping segments. But, it made sense later: pgRouting does not present the linestrings with their vertices sorted from the starting point to the ending point of the route. At a transition from one linestring to another, there are 4 combinations of vertex layout: the previous segment's first vertex may be the one closer to the starting point, or it may be the one closer to the end point, and the current segment's first vertex may be the one closer to the starting point or it may be the one closer to the end point.

If we do the azimuth calculation without considering this, we could get an azimuth the exact opposite of correct: the user is heading east on a street but we measured the azimuth from end to start, so got an azimuth of 270 instead of 90.

So, having fetched the vertices in Phase 1, Phase 2 checks the segments' alignment using distance. The proper case is when the previous segment's final vertex ($prevvert2) is closest to the current segment's starting vertex ($thisvert1). The code in Phase 2 does distance calculations, and then inverts the vertices appropriately, so we know that the route properly takes us through $prevvert1, $prevvert2, $thisvert1, and $thisvert2 in that order.

No comments:

Post a Comment