ESTIMATING HIKING TIME ALONG A SEGMENT - TOBLER


The simplest way to calculate travel time over your route segments, is to make an assumption about walking speed, e.g. a constant 3 miles per hour. However, this client wanted to account for slope a la Tobler's hiking function, since this is substantially hilly terrain and a constant walking speed is not necessarily realistic.

This process requires 3D data, containing elevation. If your data is only 2D, get a DEM and load up ArcGIS and produce 3D data, aka LINESTRING Z. Don't worry, I'll wait...

And now it's fairly easy:

-- Supplement each trail segment with a length. The example here transforms the geometry to
-- a State Plane SRS with feet as units.
ALTER TABLE trails ADD COLUMN length_feet float;
UPDATE trails SET length_feet = ST_LENGTH(ST_TRANSFORM(the_geom,3734));

-- Now the Tobler calculation. We fetch the elevation of the first and last
-- vertices of the linestring, divide by length, and that's the slope.
-- Remember, we split our segments into 100-foot lengths so this is mostly accurate.
-- If you keep yours as mile-long trails your accuracy may vary.
-- The 0.911344 converts from Tobler's own units (kilometers per hour)
-- to our preferred units (feet per second, to match the length in feet)
ALTER TABLE routing_trails ADD COLUMN tobler_fps FLOAT;
UPDATE routing_trails SET tobler_fps =
0.911344 * 6 * EXP(
-3.5 * ABS(
    (
      ST_Z(ST_StartPoint(ST_TRANSFORM(the_geom,3734)))
      -
      ST_Z(ST_EndPoint(ST_TRANSFORM(the_geom,3734)))
    )/length )
);

-- Now that each segment has its speed in feet per second and its length in feet,
-- the duration in seconds is simple. In our case, we calculate for three different
-- travel modes and we assume that bicyclists are being safe and only doing 3X as fast as walking.
UPDATE routing_trails SET duration_hike = length_feet / tobler_fps;
UPDATE routing_trails SET duration_bike = duration_hike * 0.33;
UPDATE routing_trails SET duration_horse = duration_hike * 0.66;

Discussion

I don't really know how accurate Tobler's function is, nor the degree to which it makes a difference in the accuracy of route duration estimates. Still, our original approach is to presume 3 miles per hour (4.4 feet per second) on foot, 3X that for bicycles, and 1.5X that for equestrian. Taking the slope into account may miss a lot of factors, but it can't be any worse!

TURNING DIRECTIONS

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

SELECT
ST_AsText(ST_Transform(route.the_geom,4326)) AS wkt,
ST_Length(route.the_geom) AS length,
hiking_trails.name, hiking_trails.duration
FROM
hiking_trails,
(SELECT gid, the_geom FROM astar_sp_delta_directed('hiking_trails',?,?,2640,true,true)) AS route
WHERE
hiking_trails.gid=route.gid

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];
                    break;
                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];
                    break;
                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];
                    break;
                case $dx21:
                    // current start is previous end, already fine
                    break;
            }

            // 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.

Leaflet 0.4, Tried and Failed !

This big app I've been writing lately, has been in Leaflet 0.3 I was urged to upgrade to 0.4 for a few minor things: panning intertia, the Retina support, the new Control.Scale which is prettier than mine, and general bugfixes.

But it didn't work out, despite my careful reading and note-taking of the Breaking Changes.

On my Android (2.3.6 on a LG Spectrum), I had serious problems.

• Inertia didn't work at all.
• Animations didn't happen anymore; apparently this is intentional but they worked in 0.3.
• After loading, the map would become obscured by something, a DIV I couldn't identify, and would only reappear after I click the zoom out button, though tapping would perform click-queries and show popups over the blankness.
• The window.scrollTo() trick to hide the address bar, caused the page to flash annoyingly for some 15 seconds.
• And then it crashed the browser completely.
  

So yeah, at least for this non-trivial mobile app, Leaflet 0.3 is where we'll stay.

Unfortunately, Android was the last of the platforms I tested, and the app did work perfectly on Desktop, and on the iPad 4. So in the meantime I did discover two features which may not have been entirely clear from the "breaking changes" guide:

• Layers can now be given an explicit "zIndex" property in the constructor options. Excellent! Less excellent, is that you must do this or else the layers stack randomly; in my case, the basemaps came in on top completely obscuring the overlay layers. I know, it's not the documented behavior and they should have stacked in the order they were added, but I could make the randomness happen reliably here.
• Markers can no longer have "draggable" set unless "clickable" is also set. Not a problem, but it did cost me some 15 minutes. True to their word, clickable:false means NO mouse events, and I can see where this is a win for efficiency.