Introduction

The relationship between slope and travel rates is, generally, that steeper slopes, both uphill and downhill, tend to slow one down while traveling on foot. The specific quantitative nature of that relationship has been explored in many studies, including several of our own. Most resulting slope-travel rate functions are asymmetrical or anisotropic, meaning that the same slope uphill and downhill will produce somewhat different travel rates. This makes intuitive sense – traveling on steep slopes is difficult in both directions, but is less physically demanding downhill than uphill. However, the Estimated Ground Evacuation Time layer (EGET) is fundamentally a bi-directional product. Its main function is to provide fire crews with an approximation of how long it will take them to travel to the nearest medical facility from any location on the landcsape. However, it can also be viewed as a useful tool for examining how long it will take a crew to reach a strategic position on the landscape to maximize effective (and safe) fire response. Add to that the quantitative complexity of applying anisotropic cost accumulation on a broad spatial scale and what emerges is a clear need for a symmetrical/isotropic travel rate function. Even though it may not be the truest representation of the fundamental relationship between slope and travel rates, it will represent a useful approximation of that relationship that can be broadly and efficiently applied for the new generation of EGET.

The objective of this document is to provide an overview of a new analysis aimed at generating a novel, symmetrical slope-travel rate function. The best source of data we have, to date, to perform this analysis comes from AllTrails. AllTrails is a website and mobile app geared towards outdoor recreation, where hikers can track their hikes using their phone’s GPS. In 2022, we published a paper that used AllTrails GPS data from nearly 2000 hikes on a diverse range of trails in Utah and California, comprising over 4M GPS points on trails ranging from 2.8 km to 17.9 km from a broad swath of the population. To capture a range in travel rates (from slower to faster people), we used a quantile regression approach, from the 2.5th to the 97.5th percentile at an interval of 2.5 percentiles. The 50th percentile, representing the median travel rate, is the most useful of these functions, in that it should capture the “average” hiker’s slope-travel rate relationship. These functions demonstrated an impressive capacity to capture the range in total hike times for a set of totally independent hikes from different states around the country, with <10% prediction error. However, the functions we introduced were asymmetrical. In this study, we will use the same data and a similar modeling approach to generate a new, symmetrical slope-travel rate function for use in the new generation of EGET.

Analysis

Step 1. Training/test split

First, I will read in the data, and perform an 80/20% training/test split by GPS track (not GPS point, to avoid within-track autocorrelation).

# load libraries
library(quantreg)
## Warning: package 'quantreg' was built under R version 4.2.2
## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
## 
##     backsolve
library(data.table)
## Warning: package 'data.table' was built under R version 4.2.3
library(magrittr)
library(viridis)
## Warning: package 'viridis' was built under R version 4.2.1
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.2.3
# set random seed
set.seed(5757)

# read in the data
df <- fread("S:/ursa/campbell/alltrails/modeling/alltrails_data_all.csv") %>%
  as.data.frame()

# print first few rows of data
head(df)
##                                track.id     hiker.name        datetime.loc
## 1 bells_canyon_trail_to_lower_falls_001 thomas_norbutt 2021-06-01 21:40:16
## 2 bells_canyon_trail_to_lower_falls_001 thomas_norbutt 2021-06-01 20:22:37
## 3 bells_canyon_trail_to_lower_falls_001 thomas_norbutt 2021-06-01 20:11:18
## 4 bells_canyon_trail_to_lower_falls_001 thomas_norbutt 2021-06-01 18:49:09
## 5 bells_canyon_trail_to_lower_falls_001 thomas_norbutt 2021-06-01 21:54:41
## 6 bells_canyon_trail_to_lower_falls_001 thomas_norbutt 2021-06-01 20:11:47
##   utm.e.snp utm.n.snp elev.lidar     slope travel.rate.2d travel.rate.3d
## 1  432606.7   4491412   1656.934 -28.75542      0.2428571      0.2770188
## 2  434611.9   4490306   1986.207 -27.83946      0.2200000      0.2487957
## 3  434767.3   4490248   2049.606 -27.55096      0.2666667      0.3007746
## 4  433410.7   4490499   1780.281 -25.60598      0.3750000      0.4158408
## 5  432568.8   4491524   1595.939 -25.49895      0.2272727      0.2517997
## 6  434767.3   4490247   2049.577 -25.20079      0.3111111      0.3438371
##   travel.rate.2d.med
## 1          0.2314286
## 2          0.2428571
## 3          0.2433333
## 4          0.3111111
## 5          0.3111111
## 6          0.3111111
# create train/test lookup table
df.track <- data.frame(track.id = unique(df$track.id))
df.track$train.test <- "train"
df.track$train.test[sample(nrow(df.track), floor(nrow(df.track) * 0.2))] <- "test"

# join lut to df and split into train/test
df <- merge(df, df.track)
df.train <- df[df$train.test == "train",]
df.test <- df[df$train.test == "test",]

Step 2. Generate symmetrical function

The previous, asymmetrical function took the following form:

\(r=c \left(\frac{1}{\pi b\left(1+\left(\frac{s-a}{b}\right)^{2}\right)}\right)+d+es\)

where \(r\) is travel rate in m/s, \(s\) is slope in degrees, and \(a\), \(b\), \(c\), \(d\), and \(e\) are all model coefficients derived from quantile regression. The terms that make this function asymmetrical are \(a\), which provides an offset that defines where the “peak” of the travel rate function is (the slope of maximum travel rate), and \(e\) which effectively “tilts” the whole curve with a linear term applied to slope. Thus, by removing those terms, we can create an asymmetrical function:

\(r=c \left(\frac{1}{\pi b\left(1+\left(\frac{s}{b}\right)^{2}\right)}\right)+d\)

So, provided that we can find suitable coefficient values for \(b\), \(c\), and \(d\), we should be able to generate a new, symmetrical function. To mimic the previous analysis, I will use a quantile regression approach. However, I will only be generating the 50th percentile. It’s very similar to a more common ordinary least squares regression, but whereas OLS will try to fit to the mean of the sample, quantile will try to fit to the median, which may be better suited for non-normally distributed datasets. The nlrq() (non-linear quantile regression) function requires that you provide it with starting estimates for the function terms, which I will grab from our previous analysis’ median function.

# define starting terms (note that I'm including a and e for now for later use)
a <- -1.4579
b <- 22.0787
c <- 76.3271
d <- 0.0525
e <- -3.2002e-4

# generate symmetrical model
mod <- nlrq(travel.rate.2d.med ~ c*(1/(pi*b*(1+(slope/b)^2)))+d,
            data = df.train, start = list(b = b, c = c, d = d), tau = 0.5)

# get list of coefficients
cf <- coefficients(mod)
b.new <- cf[[1]]
c.new <- cf[[2]]
d.new <- cf[[3]]
print(cf)
##          b          c          d 
## 22.9376861 81.7522684  0.0112471

Step 3. Compare symmetrical and asymmetrical functions

Now that we have our three coefficients, we can use them to visualize the resulting function form and assess the performance of the new model. We can also compare to our previous function. So, let’s do that!

# get travel rates on a range of slopes for the new function
new.fun <- function(slope, b, c, d){
  rate <- c*(1/(pi*b*(1+(slope/b)^2)))+d
  return(rate)
}
slopes <- seq(-90,90)
rates.new <- new.fun(slopes, b.new, c.new, d.new)

# get travel rates on a range of slopes for the old function
old.fun <- function(slope, a, b, c, d, e){
  rate <- c*(1/(pi*b*(1+((slope-a)/b)^2)))+d+e*slope
  return(rate)
}
rates.old <- old.fun(slopes, a, b, c, d, e)

# set up plot
par(mar = c(5,5,1,1), las = 1)

# plot them both out
xlim <- c(-90,90)
ylim <- c(min(c(rates.new, rates.old)), max(c(rates.new, rates.old)))
plot(ylim ~ xlim, type = "n", xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
lines(rates.old ~ slopes, lwd = 3, col = 2)
lines(rates.new ~ slopes, lwd = 3, col = 4)
legend("topright", legend = c("Old (asymmetrical)", "New (symmetrical)"),
       lwd = 3, col = c(2,4))

Visually, the main takeaway should be abundantly clear: there is extremely little difference between the asymmetrical (old) and symmetrical (new) functions. This is great news! That means the resulting travel rate predictions should produce very similar results. Let’s run a quick test on an imaginary trail to see how much the total travel time estimates would vary:

# create an imaginary trail and plot it out
slopes <- rnorm(1000, sd = 10)
dists <- rep(10, 1000)
cumdist <- cumsum(dists)
par(mar = c(5,5,1,1), las = 1)
plot(slopes ~ cumdist, type = "l", xlab = "Distance Along Trail (m)",
     ylab = "Slope (deg)", col = 3)
grid()
abline(h = 0)

…Obviously an unrealistic trail, but captures a realistic range of slopes on a real trail…

# get difference in cumulative travel time predictions
time.old <- cumsum(dists * (1 / old.fun(slopes, a, b, c, d, e)))
time.new <- cumsum(dists * (1 / new.fun(slopes, b.new, c.new, d.new)))
time.diff <- time.new - time.old

# plot the results out
par(mar = c(5,5,2,1), las = 1)
plot(time.old ~ cumdist, type = "l", col = 2, xlab = "Distance Along Trail (m)",
     ylab = "Cumulative Travel Time (s)")
lines(time.new ~ cumdist, col = 4)
grid()
legend("topleft", legend = c("Old (asymmetrical)", "New (symmetrical)"),
       lwd = 1, col = c(2, 4))

In the end, on a topographically diverse trail with a total length of 10km (6.2 mi), the two functions produced a difference of only 17s, or 0.1639028% difference relative to the old function.

Conclusion

We now have a new mathematical function that represents a symmetrical (isotropic, direction-independent) relationship between slope and travel rate. The great news is that this function produces nearly identical results to the asymmetrical function, which we have already demonstrated to be good at estimating hike times among broad populations. Here is the final function form:

\(r=77.61962\left(\frac{1}{\pi\cdot22.4056\left(1+\left(\frac{s}{22.4056}\right)^{2}\right)}\right)+0.04644277\)

We can simplify this to:

\(r=\frac{0.006512s^{2}+80.88869}{0.140214s^{2}+70.389209}\)