Reexamining Levan Travel Rates
Mickey Campbell
2023-06-21
Introduction
In the Fall of 2016, I led a field experiment whereby 31 study subjects were timed as they walked along 22 transects, each 100m in length, near Levan, UT. The transects were strategically placed on the landscape so as to capture a variety of terrain and vegetation conditions. Of particular interest were: (1) the slope of the terrain; (2) the ground surface roughness; and (3) the density of the vegetation. These three structural parameters, derived from airborne lidar data, were compared to travel rates to gain novel insight into the landscape conditions’ effects on movement. We published a paper with the study’s main findings, which can be read here.
Although the study was successful in achieving its objectives, it was not without its limitations. One of the most important limitations is that the relationship between slope and travel rate was assumed to take the form of a negative exponential. Here was the equation from the paper for predicting travel rates:
\(travelrate = \alpha + \beta_{1}*density + \beta_{2}*roughness + \beta_{3}*slope + \beta_{4}*slope^{2}\)
where \(\alpha\) and \(\beta\) are regression coefficients. All of the above are explained in greater detail (units, formulas, etc.) in the previously linked paper. Let’s begin by defining this travel rate function as an R function:
# define slope-travel rate function
a <- 1.662
b1 <- -1.076
b2 <- -9.011
b3 <- -5.191e-3
b4 <- -1.127e-3
travelrate.fun <- function(density,roughness,slope){
travelrate <- a + b1 * density + b2 * roughness + b3 * slope + b4 * slope ^ 2
return(travelrate)
}To demonstrate the problem with the slope element of the function form defined above, let’s assume \(density\) and \(roughness\) are both zero and plot the relationship between slope and travel rates:
# predict travel rates without vegetation or roughness
slopes <- seq(-45,45,0.1)
rates <- travelrate.fun(0,0,slopes)
# plot it out
par(las = 1, mar = c(5,5,1,1))
plot(rates ~ slopes, type = "l", lwd = 2, col = "red",
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
abline(h = 0)
As can be seen above, the function above assumes that travel rates switch to negative at about -40.8 and 36.2. Of course, there are no such thing as negative travel rates, which is one major reason why this function form is problematic. In more recent work (here and here), we have identified the Lorentz function form as being a better representation of the slope-travel rate relationship. That function is defined as:
\(travelrate = c\left(\frac{1}{\pi*b\left(1+(\frac{slope-a}{b})^{2}\right)}\right)+d+e*slope\)
where \(a\), \(b\), \(c\), \(d\), and \(e\) are all coefficients. Through quantile regression modeling of crowdsourced GPS data from AllTrails hikes, we have defined these coefficients for a range of travel rate percentiles, ranging from those who move relatively slow to those who move relatively fast. For the sake of demonstration, we can use the coefficients from the 50th percentile model, which represents the median among a diverse population and should represent the “average” hiker pretty well:
library(rmarkdown)
# define lorentz slope-travel rate function
travelrate.fun <- function(slope,a,b,c,d,e){
travelrate <- c * (1 / (pi * b * (1 + ((slope - a) / b) ^ 2))) + d + e * slope
return(travelrate)
}
# read in table of coefficients and print them out
slope.rate.terms <- read.csv("S:/ursa/campbell/alltrails/modeling/slope_vs_travel_rate_model_percentile_coefficients_med.csv")
paged_table(slope.rate.terms)# get the median terms
a <- slope.rate.terms$a[slope.rate.terms$percentile == 0.5]
b <- slope.rate.terms$b[slope.rate.terms$percentile == 0.5]
c <- slope.rate.terms$c[slope.rate.terms$percentile == 0.5]
d <- slope.rate.terms$d[slope.rate.terms$percentile == 0.5]
e <- slope.rate.terms$e[slope.rate.terms$percentile == 0.5]
# plot out slope-travel rate relationship for median
slopes <- seq(-45,45,0.1)
rates <- travelrate.fun(slopes,a,b,c,d,e)
par(las = 1, mar = c(5,5,1,1))
plot(rates ~ slopes, type = "l", lwd = 2, col = "red",
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
As you can see, not only does the travel rate never reach zero, even at fairly extreme slopes (+/- 45 degrees), but it also has a more logically-consistent form. As you get steeper uphill and downhill from the slope of peak travel rate, rates decrease initially quite sharply, but then begin flatten out on steep slopes.
This brings us to the main objective of this markdown document:
To re-examine the Levan travel rate data, replacing the negative exponential function with the Lorentz function, informed by the coefficients gleaned from our AllTrails data
This will hopefully provide us with a new equation that we can use for Allison’s research examining tradeoffs between travel time and situational awareness for wildland firefighters.
Data Preparation
First things first, I need to load in the travel rate experiment data:
library(stringr)
# read in the travel rate data
df <- read.csv("S:/ursa/campbell/Other/levan/data/experiment_data/final_data_table.csv")
# anonymize study subjects
names <- unique(df$name)
subjs <- paste0("subject ", str_pad(seq(1,length(names)), 2, "left", "0"))
lut <- data.frame(name = names, subject = subjs)
df <- merge(df, lut)
df$name <- NULL
# clean up table and print it out
df <- df[,c("transect", "subject", "rate", "slope", "nrd_0.15_2.75", "roughness")]
colnames(df) <- c("transect", "subject", "rate", "slope", "density", "roughness")
df <- df[order(df$transect, df$subject, df$slope),]
row.names(df) <- seq(1,nrow(df))
paged_table(df)In the original Levan study, we used nonlinear mixed effects regression, which allows you to account for individual-level differences in travel rates. For example, Joe may always travel faster than John, but the relative effects of slope on both Joe and John are the same. Mixed effects regression allows you to enforce that Joe’s individual model will have a higher y-intercept (i.e., a higher baseline travel rate) than John’s model, while modeling the relationship between slope and travel rate. Although this is a nice method in this type of study, for the sake of simplicity, let’s collapse all of the study subjects’ travel rates down to a per-transect/per-direction mean. We have 22 transects in total, and study subjects traveled once in each direction, so we will wind up with 44 data points to drive our new model.
# get mean travel rates per transect and direction
df.mean <- aggregate(x = list(rate.mean = df$rate),
by = list(transect = df$transect,
slope = df$slope,
density = df$density,
roughness = df$roughness),
FUN = mean)
df.mean <- df.mean[order(df.mean$transect, df.mean$slope),]
paged_table(df.mean)Modeling
We now have a dataset ready for modeling. We will use non-linear regression modeling to fit these data to a new function that hopefully captures the relationship between travel rates, slope, vegetation density, and ground surface roughness in a new and more accurate way. Non-linear regression models, however, are quite complex computationally, particularly in situations where you’re trying to model complex function forms like the Lorentz function. Especially with only 44 data points, it may be difficult for the model to “converge” (essentially, find a combination of regression equation terms that are optimized to fit your dataset). To aid in this convergence process, you need to provide a good set of approximate starting terms that are within the ballpark of what the final terms might actually look like. It’s basically like giving the model a “hint” as to where to start when looking for the best-fitting coefficients. The challenge here is that the Levan data are somewhat unique in comparison to “typical” hiking data – this is another limitation of that study. Likely because the transects were so short (100m), study subjects walked them quite fast in comparison to what one might hike over longer distances.
To illustrate this, I’ll plot out the travel rates versus terrain slope for the transects with lowest vegetation density (\(density\) < 0.1) and compare that to the AllTrails travel rate percentile functions:
library(viridis)## Loading required package: viridisLite# subset the data to just the sparsest transects
df.sparse <- df.mean[df.mean$density < 0.1,]
# plot travel rate vs. slope in comparison to AllTrails percentile models
par(las = 1, mar = c(5,5,1,1))
plot(rate.mean ~ slope, data = df.sparse, type = "n",
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
cols <- magma(nrow(slope.rate.terms))
for (i in seq(1,nrow(slope.rate.terms))){
slopes <- seq(-20,20,0.1)
a <- slope.rate.terms$a[i]
b <- slope.rate.terms$b[i]
c <- slope.rate.terms$c[i]
d <- slope.rate.terms$d[i]
e <- slope.rate.terms$e[i]
rates <- travelrate.fun(slopes,a,b,c,d,e)
lines(rates ~ slopes, col = cols[i], lwd = 2)
}
points(rate.mean ~ slope, data = df.sparse, pch = 21, bg = "cyan")
For reference, the highest slope-travel rate model line in light yellow represents the 97.5th travel rate percentile in the AllTrails data. Each line below it sequentially represents a 2.5 percentile decrease from there (95th, 92.5th, etc.). And the points represent the travel rates from the Levan experiment for the sparsest transects. So, clearly, in comparison to a much more broadly-representative group (the AllTrails hiker database), the Levan study subjects were moving fast. Thus, to provide the nonlinear regression model with a useful set of starting coefficients, we will supply it with those 97.5th percentile model. We also need to supply starting coefficients for the two additional model terms: the vegetation density (whose coefficient we will call \(f\)) and the ground surface roughness (\(g\)). For these, we will use the same coefficients as were found in the original paper (\(f = -1.076\) and \(g = -9.011\)). This will bring our new function form to:
\(v = c\left(\frac{1}{\pi b\left(1+(\frac{slope-a}{b})^{2}\right)}\right)+d+e*slope+f*density+g*roughness\)
Let’s give it a whirl!
# define starting coefficients
a.start <- slope.rate.terms$a[slope.rate.terms$percentile == 0.975]
b.start <- slope.rate.terms$b[slope.rate.terms$percentile == 0.975]
c.start <- slope.rate.terms$c[slope.rate.terms$percentile == 0.975]
d.start <- slope.rate.terms$d[slope.rate.terms$percentile == 0.975]
e.start <- slope.rate.terms$e[slope.rate.terms$percentile == 0.975]
f.start <- -1.076
g.start <- -9.011
# build model
mod.full <- nls(rate.mean ~ (c*(1/(pi*b*(1+((slope-a)/b)^2)))+d+e*slope) + f * density + g * roughness,
data = df.mean,
start = list(a = a.start,
b = b.start,
c = c.start,
d = d.start,
e = e.start,
f = f.start,
g = g.start),
control = nls.control(maxiter = 1000))It converged! Which is good. Let’s take a look at the initial versus resulting coefficients:
# compare initial and resulting coefficients
df.coef <- data.frame(coef.init = c(a.start,
b.start,
c.start,
d.start,
e.start,
f.start,
g.start),
coef.rslt = coef(mod.full))
paged_table(df.coef)Fairly similar. To test the performance of this new model, let’s perform a leave-one-out cross-validation. That is, iterate through each point in the dataset, and on each iteration, leave that point aside, build a model using the remaining 43 points, and predict the value for that left-out point. The compiled predicted vs. observed values for each left out point can then be used to generate error metrics such as R2 and RMSE. For each iteration in the cross-validation, we will use the coefficients from the full model above as the starting points.
# get starting coefficients
a.start <- coef(mod.full)[1]
b.start <- coef(mod.full)[2]
c.start <- coef(mod.full)[3]
d.start <- coef(mod.full)[4]
e.start <- coef(mod.full)[5]
f.start <- coef(mod.full)[6]
g.start <- coef(mod.full)[7]
# leave one out cross-validation
for (i in seq(1, nrow(df.mean))){
# split into training and test
df.train <- df.mean[-i,]
df.test <- df.mean[i,]
# build model
mod.fold <- nls(rate.mean ~ (c*(1/(pi*b*(1+((slope-a)/b)^2)))+d+e*slope) + f * density + g * roughness,
data = df.train,
start = list(a = a.start,
b = b.start,
c = c.start,
d = d.start,
e = e.start,
f = f.start,
g = g.start),
control = nls.control(maxiter = 1000))
# apply model
pred <- predict(mod.fold, df.test)
obs <- df.test$rate.mean
# compile results
if (i == 1){
results.df <- data.frame(pred = pred,
obs = obs,
a = coef(mod.fold)[1],
b = coef(mod.fold)[2],
c = coef(mod.fold)[3],
d = coef(mod.fold)[4],
e = coef(mod.fold)[5],
f = coef(mod.fold)[6],
g = coef(mod.fold)[7])
} else {
results.df <- rbind(results.df, data.frame(pred = pred,
obs = obs,
a = coef(mod.fold)[1],
b = coef(mod.fold)[2],
c = coef(mod.fold)[3],
d = coef(mod.fold)[4],
e = coef(mod.fold)[5],
f = coef(mod.fold)[6],
g = coef(mod.fold)[7]))
}
}Modeling
OK, time to check out the results! First, we’ll do a basic predicted vs. observed plot, based on the leave-one-out cross-validation.
# set up plot
par(las = 1, mar = c(5,5,1,1))
# define point size and color
pt.sizes <- 1 + df.mean$density * 4
pt.cols <- colorRampPalette(c("blue", "white", "red"))(nrow(df.mean))
pt.cols <- pt.cols[as.numeric(cut(df.mean$slope, breaks = 44))]
# define x and y
x <- results.df$obs
y <- results.df$pred
# get axis limits
ax.min <- min(c(x, y))
ax.max <- max(c(x, y))
# plot the data
plot(y ~ x, type = "n", xlab = "Observed Travel Rate (m/s)", ylab = "Predicted Travel Rate (m/s)", xlim = c(ax.min, ax.max), ylim = c(ax.min, ax.max))
grid()
lines(c(-100,100),c(-100,100), col = "darkgray", lty = 2)
points(y ~ x, pch = 21, bg = pt.cols, cex = pt.sizes)
# get model performance metrics
mod <- lm(y ~ x)
abline(mod, lwd = 2)
r2 <- round(summary(mod)$adj.r.squared,2)
rrmse <- round(100*(sqrt(mean((y-x) ^ 2)))/mean(x),2)
r2.txt <- bquote(R^2==.(r2))
rrmse.txt <- bquote(rRMSE==.(rrmse)*"%")
a <- round(coef(mod)[2], 2)
b <- round(coef(mod)[1], 2)
if (b < 0){
eq.txt <- bquote(y==.(a)~x-.(abs(b)))
} else {
eq.txt <- bquote(y==.(a)~x+.(abs(b)))
}
legend("topleft",
bty = "n",
legend = c(eq.txt, r2.txt, rrmse.txt),
x.intersp = 0)
# add legend for slope
pt.cols <- colorRampPalette(c("blue", "white", "red"))(nrow(df))
leg.img <- as.raster(matrix(rev(pt.cols), ncol = 1))
rasterImage(image = leg.img,
xleft = par("usr")[1] + 0.9 * (par("usr")[2] - par("usr")[1]),
ybottom = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
xright = par("usr")[1] + 0.95 * (par("usr")[2] - par("usr")[1]),
ytop = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]))
rect(xleft = par("usr")[1] + 0.9 * (par("usr")[2] - par("usr")[1]),
ybottom = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
xright = par("usr")[1] + 0.95 * (par("usr")[2] - par("usr")[1]),
ytop = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]))
text(x = par("usr")[1] + 0.95 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.3 * (par("usr")[4] - par("usr")[3]),
labels = "Slope (deg)",
adj = c(1,0))
text(x = par("usr")[1] + 0.88 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]),
labels = round(max(df$slope)),
adj = c(1,0.5))
text(x = par("usr")[1] + 0.88 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.15 * (par("usr")[4] - par("usr")[3]),
labels = "0",
adj = c(1,0.5))
text(x = par("usr")[1] + 0.88 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
labels = round(min(df$slope)),
adj = c(1,0.5))
# add legend for density
text(x = par("usr")[1] + 0.75 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.3 * (par("usr")[4] - par("usr")[3]),
labels = "Density",
adj = c(1,0))
text(x = par("usr")[1] + 0.68 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]),
labels = str_pad(round(max(df$density),2), 4, "right", "0"),
adj = c(1,0.5))
text(x = par("usr")[1] + 0.68 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.15 * (par("usr")[4] - par("usr")[3]),
labels = str_pad(round(mean(range(df$density)),2), 4, "right", "0"),
adj = c(1,0.5))
text(x = par("usr")[1] + 0.68 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
labels = str_pad(round(min(df$density),2), 4, "right", "0"),
adj = c(1,0.5))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]),
cex = max(pt.sizes))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.2 * (par("usr")[4] - par("usr")[3]),
cex = mean(c(max(pt.sizes), mean(range(pt.sizes)))))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.15 * (par("usr")[4] - par("usr")[3]),
cex = mean(range(pt.sizes)))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.1 * (par("usr")[4] - par("usr")[3]),
cex = mean(c(min(pt.sizes), mean(range(pt.sizes)))))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
cex = min(pt.sizes))
Truly a thing of beauty. I was shocked when I first saw these results and how well they are performing. One of the cool things about this leave-one-out approach is that we can examine the variability in the model coefficients. This can provide us with somewhat of a sense of the confidence in the coefficients for things like vegetation density and ground surface roughness. The slope coefficients are a little bit harder to nail down in terms of their meaning, but for the sake of completeness, I’ll just plot them all out as density distributions.
library(RColorBrewer)
# define function for adding transparency to colors
add.alpha <- function(col, alpha=0.5){
apply(sapply(col, col2rgb)/255, 2,
function(x)
rgb(x[1], x[2], x[3], alpha=alpha))
}
# define colors
cols.fill <- add.alpha(brewer.pal(7,"Set1"))
cols.line <- brewer.pal(7,"Set1")
# plot density distributions
par(las = 1, mfrow = c(3,3), mar = c(5,5,1,1))
for (i in seq(3,ncol(results.df))){
x <- results.df[,i]
d <- density(x)
plot(d$y ~ d$x, type = "n", xlab = colnames(results.df)[i], ylab = "Density")
polygon(d$y ~ d$x, col = cols.fill[i-2], border = NA)
lines(d$y ~ d$x, col = cols.line[i-2], lwd = 2)
abline(v = median(x), lty = 2, lwd = 2)
if (median(x) < mean(par("usr")[1:2])){
legend("topright", legend = paste0("Med = ", round(median(x),3)),
x.intersp = 0, bty = "n")
} else {
legend("topleft", legend = paste0("Med = ", round(median(x), 3)),
x.intersp = 0, bty = "n")
}
}
Conclusion
So there we have it! The new equation for predicting
travel rates as a function of slope, vegetation density, and ground
surface roughness is:
\(v = c\left(\frac{1}{\pi
b\left(1+(\frac{slope-a}{b})^{2}\right)}\right)+d+e*slope+f*density+g*roughness\)
where:
\(a\) = 0.191\(b\) = 10.455\(c\) = 13.269\(d\) = 1.238\(e\) = -0.005\(f\) = -1.05\(g\) = -7.375
Update Upon Further Reflection (6/22/2023)
So, I may have gotten a little bit ahead of myself there… Always a risky move to add a “Conclusion” section to ongoing work, I suppose. I wanted to explore what the resulting predictive equation looks like under different scenarios, and some issues arose. Let’s assume, for a minute, that the results gleaned from the analysis above are valid. Under that assumption, let’s plot out some simple travel rate predictions based on a range of slopes and vegetation densities. For the sake of example, we’ll assume no ground surface roughness.
# define new travel rate function that incorporates slope, vegetation density, and surface roughness
travelrate.fun <- function(slope,density,roughness,a,b,c,d,e,f,g){
travelrate <- c * (1 / (pi * b * (1 + ((slope - a) / b) ^ 2))) + d + e * slope + f * density + g * roughness
return(travelrate)
}
# get coefficients (medians from LOOCV)
a <- median(results.df$a)
b <- median(results.df$b)
c <- median(results.df$c)
d <- median(results.df$d)
e <- median(results.df$e)
f <- median(results.df$f)
g <- median(results.df$g)
# get range of travel rate predictions
min.rate <- 100
max.rate <- -100
for (i in seq(0,1,0.1)){
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
if (min(rates) < min.rate) min.rate <- min(rates)
if (max(rates) > max.rate) max.rate <- max(rates)
}
# plot out the results
cols <- viridis(11)
counter <- 0
par(las = 1, mar = c(5,5,1,8))
plot(x = c(-45,45), y = c(min.rate, max.rate), type = "n",
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
for (i in seq(0,1,0.1)){
counter <- counter + 1
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
lines(rates ~ slopes, lwd = 2, col = cols[counter])
}
par(xpd = NA)
legend(x = par("usr")[2], y = mean(par("usr")[3:4]),
legend = seq(0,1,0.1), title = "Veg Dens", lwd = 2, col = cols,
bty = "n", yjust = 0.5)
A careful (or perhaps not even all that careful) look at these curves
will raise some red flags. The most prominent of those red flags is that
the linear asymmetry term in the slope-travel rate function (\(e\)) has resulted in a pretty strong “tilt”
to these curves, making downhill travel generally quite a bit faster
than uphill travel. I can agree with that in principle, but
quantitatively, the strength of that tilt becomes problematic for the
prediction of downhill travel rates. Let’s dig a little deeper:
# plot travel rates vs. slope for a single veg density (0.5)
slopes <- seq(-45,45,0.1)
densities <- rep(0.5, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
par(las = 1, mar = c(5,5,1,1))
plot(rates ~ slopes, type = "l", col = "blue", lwd = 2,
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
# find minimum downhill travel rate
y.min <- min(rates[slopes < 0])
x.min <- slopes[rates == y.min]
points(x.min, y.min, pch = 21, bg = "red", cex = 1.5)
# find travel rate at -45 degrees
x.45 <- -45
y.45 <- rates[slopes == -45]
points(x.45, y.45, pch = 21, bg = "orange", cex = 1.5)
# find downhill slope at which travel rate equals that of -45 degrees
slopes.sub <- slopes[slopes < 0 & slopes > x.min]
rates.sub <- rates[slopes < 0 & slopes > x.min]
rates.sub.diff <- abs(rates.sub - y.45)
y.equiv <- rates.sub[which.min(rates.sub.diff)]
x.equiv <- slopes.sub[which.min(rates.sub.diff)]
points(x.equiv, y.equiv, pch = 21, bg = "yellow", cex = 1.5)
# add legend
legend("bottomleft",
legend = c(paste0("Slope of Min Rate = ", x.min),
paste0("Steepest Slope = -45"),
paste0("Slope where Rate Equals Steepest Slope = ", x.equiv)),
pch = 21, pt.bg = c("red", "orange", "yellow", pt.cex = 1.5))
So, if we were to move forward with this equation, we would be assuming that people move slowest when hiking downhill at -23 degrees, and that people move at the same speed down a 45 degree slope and a -11.1 degree slope. For obvious reasons, this is problematic. If we were to use it to model least cost paths, the algorithm would choose preferentially to send firefighters down a 45 degree slope over a 23 degree slope. Clearly, this will not work.
So, what do we do?
The way I see it, there are two possible directions forward:
- Approach 1: Reanalyze the Levan data, getting rid of the
asymmetric equation terms (\(a\) and
\(e\)), enforcing a symmetrical
slope-travel rate equation.
- Pros:
- We keep everything within the context of the Levan study
- We plan to develop a symmetrical travel rate function for the Estimated Ground Evacuation Time layer project anyway, so this could be an early test
- The GIS analysis becomes much simpler, as directionality does not need to be taken into account, thus enabling all analyses to stay within the ArcGIS ecosystem
- Cons:
- Perhaps it places too much trust in the Levan data for creating (ideally) generalizable relationships?
- It takes a lot of the nuance that we’ve been able to achieve in our understanding of slope-travel rate relationships from Strava and AllTrails and says “no thanks”
- It winds up, in effect, creating a new formulation of the slope-travel rate relationship, which opens the door to scrutiny
- Pros:
- Approach 2: Use the existing analysis to extract density-
and roughness-based coefficients, but take the slope-travel rate
coefficients from AllTrails study
- Pros:
- The AllTrails slope-travel rate relationships are MUCH more robust than those from Levan for a variety of reasons (sample size, longer hikes, GPS tracked, wider range of slopes, etc.)
- The end result may be a more generalizable equation that is broadly applicable
- We can leverage the percentile nature of the AllTrails relationships to incorporate nuance into the prediction of travel times as a function of slope, density, roughness, and overall faster/slower movement
- Cons:
- It’s a sort of complex mashup between two models, driven by disconnected datasets, which may invite scrutiny
- We have to assume that the density and roughness relationships observed in Levan can be broadly applied
- Pros:
For the sake of completeness (and because I am genuinely enjoying this…) I’ll take a look at both options.
Approach #1
The goal of this approach will be to essentially repeat everything done above, but remove the two asymmetrical terms (\(a\) and \(e\)) from the travel rate equation. I’ll just dump a bunch of code in here and spare you the narrative until the results.
# define starting coefficients
b.start <- slope.rate.terms$b[slope.rate.terms$percentile == 0.975]
c.start <- slope.rate.terms$c[slope.rate.terms$percentile == 0.975]
d.start <- slope.rate.terms$d[slope.rate.terms$percentile == 0.975]
f.start <- -1.076
g.start <- -9.011
# build model
mod.full <- nls(rate.mean ~ (c*(1/(pi*b*(1+((slope)/b)^2)))+d) + f * density + g * roughness,
data = df.mean,
start = list(b = b.start,
c = c.start,
d = d.start,
f = f.start,
g = g.start),
control = nls.control(maxiter = 1000))
# get starting coefficients
b.start <- coef(mod.full)[1]
c.start <- coef(mod.full)[2]
d.start <- coef(mod.full)[3]
f.start <- coef(mod.full)[4]
g.start <- coef(mod.full)[5]
# leave one out cross-validation
for (i in seq(1, nrow(df.mean))){
# split into training and test
df.train <- df.mean[-i,]
df.test <- df.mean[i,]
# build model
mod.fold <- nls(rate.mean ~ (c*(1/(pi*b*(1+((slope)/b)^2)))+d) + f * density + g * roughness,
data = df.train,
start = list(b = b.start,
c = c.start,
d = d.start,
f = f.start,
g = g.start),
control = nls.control(maxiter = 1000))
# apply model
pred <- predict(mod.fold, df.test)
obs <- df.test$rate.mean
# compile results
if (i == 1){
results.df <- data.frame(pred = pred,
obs = obs,
b = coef(mod.fold)[1],
c = coef(mod.fold)[2],
d = coef(mod.fold)[3],
f = coef(mod.fold)[4],
g = coef(mod.fold)[5])
} else {
results.df <- rbind(results.df, data.frame(pred = pred,
obs = obs,
b = coef(mod.fold)[1],
c = coef(mod.fold)[2],
d = coef(mod.fold)[3],
f = coef(mod.fold)[4],
g = coef(mod.fold)[5]))
}
}
# set up plot
par(las = 1, mar = c(5,5,1,1))
# define point size and color
pt.sizes <- 1 + df.mean$density * 4
pt.cols <- colorRampPalette(c("blue", "white", "red"))(nrow(df.mean))
pt.cols <- pt.cols[as.numeric(cut(df.mean$slope, breaks = 44))]
# define x and y
x <- results.df$obs
y <- results.df$pred
# get axis limits
ax.min <- min(c(x, y))
ax.max <- max(c(x, y))
# plot the data
plot(y ~ x, type = "n", xlab = "Observed Travel Rate (m/s)", ylab = "Predicted Travel Rate (m/s)", xlim = c(ax.min, ax.max), ylim = c(ax.min, ax.max))
grid()
lines(c(-100,100),c(-100,100), col = "darkgray", lty = 2)
points(y ~ x, pch = 21, bg = pt.cols, cex = pt.sizes)
# get model performance metrics
mod <- lm(y ~ x)
abline(mod, lwd = 2)
r2 <- round(summary(mod)$adj.r.squared,2)
rrmse <- round(100*(sqrt(mean((y-x) ^ 2)))/mean(x),2)
r2.txt <- bquote(R^2==.(r2))
rrmse.txt <- bquote(rRMSE==.(rrmse)*"%")
a <- round(coef(mod)[2], 2)
b <- round(coef(mod)[1], 2)
if (b < 0){
eq.txt <- bquote(y==.(a)~x-.(abs(b)))
} else {
eq.txt <- bquote(y==.(a)~x+.(abs(b)))
}
legend("topleft",
bty = "n",
legend = c(eq.txt, r2.txt, rrmse.txt),
x.intersp = 0)
# add legend for slope
pt.cols <- colorRampPalette(c("blue", "white", "red"))(nrow(df))
leg.img <- as.raster(matrix(rev(pt.cols), ncol = 1))
rasterImage(image = leg.img,
xleft = par("usr")[1] + 0.9 * (par("usr")[2] - par("usr")[1]),
ybottom = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
xright = par("usr")[1] + 0.95 * (par("usr")[2] - par("usr")[1]),
ytop = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]))
rect(xleft = par("usr")[1] + 0.9 * (par("usr")[2] - par("usr")[1]),
ybottom = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
xright = par("usr")[1] + 0.95 * (par("usr")[2] - par("usr")[1]),
ytop = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]))
text(x = par("usr")[1] + 0.95 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.3 * (par("usr")[4] - par("usr")[3]),
labels = "Slope (deg)",
adj = c(1,0))
text(x = par("usr")[1] + 0.88 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]),
labels = round(max(df$slope)),
adj = c(1,0.5))
text(x = par("usr")[1] + 0.88 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.15 * (par("usr")[4] - par("usr")[3]),
labels = "0",
adj = c(1,0.5))
text(x = par("usr")[1] + 0.88 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
labels = round(min(df$slope)),
adj = c(1,0.5))
# add legend for density
text(x = par("usr")[1] + 0.75 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.3 * (par("usr")[4] - par("usr")[3]),
labels = "Density",
adj = c(1,0))
text(x = par("usr")[1] + 0.68 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]),
labels = str_pad(round(max(df$density),2), 4, "right", "0"),
adj = c(1,0.5))
text(x = par("usr")[1] + 0.68 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.15 * (par("usr")[4] - par("usr")[3]),
labels = str_pad(round(mean(range(df$density)),2), 4, "right", "0"),
adj = c(1,0.5))
text(x = par("usr")[1] + 0.68 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
labels = str_pad(round(min(df$density),2), 4, "right", "0"),
adj = c(1,0.5))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]),
cex = max(pt.sizes))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.2 * (par("usr")[4] - par("usr")[3]),
cex = mean(c(max(pt.sizes), mean(range(pt.sizes)))))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.15 * (par("usr")[4] - par("usr")[3]),
cex = mean(range(pt.sizes)))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.1 * (par("usr")[4] - par("usr")[3]),
cex = mean(c(min(pt.sizes), mean(range(pt.sizes)))))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
cex = min(pt.sizes))
OK, so as we might expect, the predictive capacity of the symmetrical function decreased, but not by a huge margin. An R2 of 0.89 is still quite something. It is quite clear, however, that more reddish dots are above the 1:1 line and more blueish dots are below the 1:1 line. Translation: uphill travel rates are being overpredicted and downhill travel rates are being underpredicted. Totally to be expected when we’re ignoring travel direction.
# define function for adding transparency to colors
add.alpha <- function(col, alpha=0.5){
apply(sapply(col, col2rgb)/255, 2,
function(x)
rgb(x[1], x[2], x[3], alpha=alpha))
}
# define colors
cols.fill <- add.alpha(brewer.pal(5,"Set1"))
cols.line <- brewer.pal(5,"Set1")
# plot density distributions
par(las = 1, mfrow = c(2,3), mar = c(5,5,1,1))
for (i in seq(3,ncol(results.df))){
x <- results.df[,i]
d <- density(x)
plot(d$y ~ d$x, type = "n", xlab = colnames(results.df)[i], ylab = "Density")
polygon(d$y ~ d$x, col = cols.fill[i-2], border = NA)
lines(d$y ~ d$x, col = cols.line[i-2], lwd = 2)
abline(v = median(x), lty = 2, lwd = 2)
if (median(x) < mean(par("usr")[1:2])){
legend("topright", legend = paste0("Med = ", round(median(x),3)),
x.intersp = 0, bty = "n")
} else {
legend("topleft", legend = paste0("Med = ", round(median(x), 3)),
x.intersp = 0, bty = "n")
}
}
To me, the striking thing here is just how consistent the vegetation density and ground surface roughness coefficients are in comparison to the previous, asymmetrical function. That says to me that these relationships are quite robust.
# define new symmetrical travel rate function that incorporates slope, vegetation density, and surface roughness
travelrate.fun <- function(slope,density,roughness,b,c,d,f,g){
travelrate <- c * (1 / (pi * b * (1 + ((slope) / b) ^ 2))) + d + f * density + g * roughness
return(travelrate)
}
# get coefficients (medians from LOOCV)
b <- median(results.df$b)
c <- median(results.df$c)
d <- median(results.df$d)
f <- median(results.df$f)
g <- median(results.df$g)
# get range of travel rate predictions
min.rate <- 100
max.rate <- -100
for (i in seq(0,1,0.1)){
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, b, c, d, f, g)
if (min(rates) < min.rate) min.rate <- min(rates)
if (max(rates) > max.rate) max.rate <- max(rates)
}
# plot out the results
cols <- viridis(11)
counter <- 0
par(las = 1, mar = c(5,5,1,8))
plot(x = c(-45,45), y = c(min.rate, max.rate), type = "n",
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
for (i in seq(0,1,0.1)){
counter <- counter + 1
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, b, c, d, f, g)
lines(rates ~ slopes, lwd = 2, col = cols[counter])
}
par(xpd = NA)
legend(x = par("usr")[2], y = mean(par("usr")[3:4]),
legend = seq(0,1,0.1), title = "Veg Dens", lwd = 2, col = cols,
bty = "n", yjust = 0.5)
Clearly, the symmetrical equation form resolves the strange downhill travel rate issues. And, importantly, the curves look logically consistent. Are they “correct”? Well, they fit the Levan data quite nicely. But, as we know, the Levan data aren’t terribly representative of broad populations.
Approach #2
As a reminder, this approach will involve extracting the density and roughness coefficients from the Levan data and merging them with the slope-travel rate relationships from the AllTrails data. To produce a generalizable relationship, I’ll use the 50th percential from the AllTrails models. Once again, code dump with some results summaries to follow.
# define new travel rate function that incorporates slope, vegetation density, and surface roughness
travelrate.fun <- function(slope,density,roughness,a,b,c,d,e,f,g){
travelrate <- c * (1 / (pi * b * (1 + ((slope - a) / b) ^ 2))) + d + e * slope + f * density + g * roughness
return(travelrate)
}
# define coefficients
a <- slope.rate.terms$a[slope.rate.terms$percentile == 0.5]
b <- slope.rate.terms$b[slope.rate.terms$percentile == 0.5]
c <- slope.rate.terms$c[slope.rate.terms$percentile == 0.5]
d <- slope.rate.terms$d[slope.rate.terms$percentile == 0.5]
e <- slope.rate.terms$e[slope.rate.terms$percentile == 0.5]
f <- median(results.df$f)
g <- median(results.df$g)
# get range of travel rate predictions
min.rate <- 100
max.rate <- -100
for (i in seq(0,1,0.1)){
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
if (min(rates) < min.rate) min.rate <- min(rates)
if (max(rates) > max.rate) max.rate <- max(rates)
}
# plot out the results
cols <- viridis(11)
counter <- 0
par(las = 1, mar = c(5,5,1,8))
plot(x = c(-45,45), y = c(min.rate, max.rate), type = "n",
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
for (i in seq(0,1,0.1)){
counter <- counter + 1
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
lines(rates ~ slopes, lwd = 2, col = cols[counter])
}
par(xpd = NA)
legend(x = par("usr")[2], y = mean(par("usr")[3:4]),
legend = seq(0,1,0.1), title = "Veg Dens", lwd = 2, col = cols,
bty = "n", yjust = 0.5)
Interesting… And potentially problematic. Why? Because of the negative travel rate problem. This says that at vegetation densities less than 20%, all slopes from -45 - +45 can be readily traversed. However, once you get to 30% or more, some portion of the curve dips into the negative travel rate zone. To take a glass half full look at this, can we just treat negative travel rates as impassible conditions? Here’s how that would look:
library(plot.matrix)
# create matrices of slopes, densities, and roughnesses
slopes <- matrix(rep(seq(-45,45), length(seq(-45,45))),
byrow = T, nrow = length(seq(-45,45)))
densities <- matrix(rep(seq(1,0,length.out = length(seq(-45,45))),
length(seq(-45,45))),
byrow = F, nrow = length(seq(-45,45)))
roughnesses <- matrix(rep(0, length(seq(-45,45)) ^ 2),
nrow = length(seq(-45,45)))
# use these matrices to create a travel rate prediction matrix
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
# define all negatives as impassible (0) and all positives as passible (1)
rates[rates <= 0] <- 0
rates[rates > 0] <- 1
# plot it out
par(las = 1, mar = c(5,5,1,8))
plot(rates, border = NA, col = c("red", "blue"),
key = NULL, axis.col = NULL, axis.row = NULL,
xlab = "Slope (deg)", ylab = "Density", main = NA)
axis(1, at = seq(1,91, length.out = 7),
labels = seq(-45,45, length.out = 7))
axis(2, at = seq(1,91, length.out = 6),
labels = seq(0,1, length.out = 6))
for (x in seq(1,91,length.out = 7)){
abline(v = x, lty = 3, col = "lightgray")
}
for (y in seq(1,91,length.out = 6)){
abline(h = y, lty = 3, col = "lightgray")
}
box()
par(xpd = NA)
legend(x = par("usr")[2], y = mean(par("usr")[3:4]),
legend = c("Passible", "Impassible"), fill = c("blue", "red"),
bty = "n", yjust = 0.5)
Looking at it this way, it’s not entirely an unreasonable set of thresholds… 100% vegetation density from 0.15 - 2.75 m basically doesn’t exist. That would require a situation where zero pulse energy can penetrate through an extremely dense understory to reach the ground surface. Highly improbable unless you’re in some tropical environment. And I kind of like the idea that, slope and vegetation density interact in this manner. On low slopes I can likely pass through even the densest vegetation. But, on steep slopes, it doesn’t take much vegetation for it to become practically impassible. Of course, all of this is extrapolation, since we don’t have any experimental data on the steepest slopes or in the densest vegetation. But, I’m not sure that we’ll ever have such data, even with new Hovermap-driven travel rate experiments. So, at some point, we have to do some extrapolation, and this does not look unreasonable to me.
Vegetation Density and Ground Surface Roughness as Multiplicative Factors
Taking a big step back from all of this, perhaps the underlying issue is how I’ve been treating vegetation density and ground surface roughness. Up to this point, everything has been additive (or, I guess, subtractive). That is, the difference in travel rates between traveling through vegetation with densities of 0 and 0.1 is the same as the difference in travel rates between traveling through vegetation with densities of 0.9 and 1 across all slopes. If this relationship was multiplicative (or, I guess, divisive), then perhaps we could avoid this negative travel rate problem altogether. Below, I’ll try to simulate what I’m talking about. Ignoring ground surface roughness for the moment, I’ll use the AllTrails median function as the basis of calculating slope-based travel rate, and then formulate an imaginary vegetation density-based multiplicative travel rate impedance factor.
# define new travel rate function that incorporates slope, vegetation density, and surface roughness
travelrate.fun <- function(slope,density,roughness,a,b,c,d,e,f){
travelrate <- (c * (1 / (pi * b * (1 + ((slope - a) / b) ^ 2))) + d + e * slope) / (f * density + 1)
return(travelrate)
}
# define coefficients
a <- slope.rate.terms$a[slope.rate.terms$percentile == 0.5]
b <- slope.rate.terms$b[slope.rate.terms$percentile == 0.5]
c <- slope.rate.terms$c[slope.rate.terms$percentile == 0.5]
d <- slope.rate.terms$d[slope.rate.terms$percentile == 0.5]
e <- slope.rate.terms$e[slope.rate.terms$percentile == 0.5]
f <- 10
# get range of travel rate predictions
min.rate <- 100
max.rate <- -100
for (i in seq(0,1,0.1)){
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f)
if (min(rates) < min.rate) min.rate <- min(rates)
if (max(rates) > max.rate) max.rate <- max(rates)
}
# plot out the results
cols <- viridis(11)
counter <- 0
par(las = 1, mar = c(5,5,1,8))
plot(x = c(-45,45), y = c(min.rate, max.rate), type = "n",
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
for (i in seq(0,1,0.1)){
counter <- counter + 1
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f)
lines(rates ~ slopes, lwd = 2, col = cols[counter])
}
par(xpd = NA)
legend(x = par("usr")[2], y = mean(par("usr")[3:4]),
legend = seq(0,1,0.1), title = "Veg Dens", lwd = 2, col = cols,
bty = "n", yjust = 0.5) 
Now this is quite an intriguing concept. The interpretation here is that, with high vegetation density, you move slowly no matter what the slope. That is, the effects of slope become muted with increasing vegetation density (which I buy!). And, no matter what the specific coefficient \(f\) is, it never goes below zero, as long as the original slope-travel rate relationship never goes below zero. Let’s see if we can make this work… I’ll start by using the asymmetrical function form, despite its shortcomings, because it fits the data the best.
# define starting coefficients
a.start <- slope.rate.terms$a[slope.rate.terms$percentile == 0.975]
b.start <- slope.rate.terms$b[slope.rate.terms$percentile == 0.975]
c.start <- slope.rate.terms$c[slope.rate.terms$percentile == 0.975]
d.start <- slope.rate.terms$d[slope.rate.terms$percentile == 0.975]
e.start <- slope.rate.terms$e[slope.rate.terms$percentile == 0.975]
f.start <- 5
g.start <- 0
# build model
mod.full <- nls(rate.mean ~ (c*(1/(pi*b*(1+((slope-a)/b)^2)))+d+e*slope) / (f * density + g * roughness + 1),
data = df.mean,
start = list(a = a.start,
b = b.start,
c = c.start,
d = d.start,
e = e.start,
f = f.start,
g = g.start),
control = nls.control(maxiter = 1000))
# get starting coefficients
a.start <- coef(mod.full)[1]
b.start <- coef(mod.full)[2]
c.start <- coef(mod.full)[3]
d.start <- coef(mod.full)[4]
e.start <- coef(mod.full)[5]
f.start <- coef(mod.full)[6]
g.start <- coef(mod.full)[7]
# leave one out cross-validation
for (i in seq(1, nrow(df.mean))){
# split into training and test
df.train <- df.mean[-i,]
df.test <- df.mean[i,]
# build model
mod.fold <- nls(rate.mean ~ (c*(1/(pi*b*(1+((slope-a)/b)^2)))+d+e*slope) / (f * density + g * roughness + 1),
data = df.train,
start = list(a = a.start,
b = b.start,
c = c.start,
d = d.start,
e = e.start,
f = f.start,
g = g.start),
control = nls.control(maxiter = 1000,
minFactor = 1e-10,
warnOnly = T))
# apply model
pred <- predict(mod.fold, df.test)
obs <- df.test$rate.mean
# compile results
if (i == 1){
results.df <- data.frame(pred = pred,
obs = obs,
a = coef(mod.fold)[1],
b = coef(mod.fold)[2],
c = coef(mod.fold)[3],
d = coef(mod.fold)[4],
e = coef(mod.fold)[5],
f = coef(mod.fold)[6],
g = coef(mod.fold)[7])
} else {
results.df <- rbind(results.df, data.frame(pred = pred,
obs = obs,
a = coef(mod.fold)[1],
b = coef(mod.fold)[2],
c = coef(mod.fold)[3],
d = coef(mod.fold)[4],
e = coef(mod.fold)[5],
f = coef(mod.fold)[6],
g = coef(mod.fold)[7]))
}
}## Warning in nls(rate.mean ~ (c * (1/(pi * b * (1 + ((slope - a)/b)^2))) + : step
## factor 5.82077e-11 reduced below 'minFactor' of 1e-10# set up plot
par(las = 1, mar = c(5,5,1,1))
# define point size and color
pt.sizes <- 1 + df.mean$density * 4
pt.cols <- colorRampPalette(c("blue", "white", "red"))(nrow(df.mean))
pt.cols <- pt.cols[as.numeric(cut(df.mean$slope, breaks = 44))]
# define x and y
x <- results.df$obs
y <- results.df$pred
# get axis limits
ax.min <- min(c(x, y))
ax.max <- max(c(x, y))
# plot the data
plot(y ~ x, type = "n", xlab = "Observed Travel Rate (m/s)", ylab = "Predicted Travel Rate (m/s)", xlim = c(ax.min, ax.max), ylim = c(ax.min, ax.max))
grid()
lines(c(-100,100),c(-100,100), col = "darkgray", lty = 2)
points(y ~ x, pch = 21, bg = pt.cols, cex = pt.sizes)
# get model performance metrics
mod <- lm(y ~ x)
abline(mod, lwd = 2)
r2 <- round(summary(mod)$adj.r.squared,2)
rrmse <- round(100*(sqrt(mean((y-x) ^ 2)))/mean(x),2)
r2.txt <- bquote(R^2==.(r2))
rrmse.txt <- bquote(rRMSE==.(rrmse)*"%")
a <- round(coef(mod)[2], 2)
b <- round(coef(mod)[1], 2)
if (b < 0){
eq.txt <- bquote(y==.(a)~x-.(abs(b)))
} else {
eq.txt <- bquote(y==.(a)~x+.(abs(b)))
}
legend("topleft",
bty = "n",
legend = c(eq.txt, r2.txt, rrmse.txt),
x.intersp = 0)
# add legend for slope
pt.cols <- colorRampPalette(c("blue", "white", "red"))(nrow(df))
leg.img <- as.raster(matrix(rev(pt.cols), ncol = 1))
rasterImage(image = leg.img,
xleft = par("usr")[1] + 0.9 * (par("usr")[2] - par("usr")[1]),
ybottom = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
xright = par("usr")[1] + 0.95 * (par("usr")[2] - par("usr")[1]),
ytop = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]))
rect(xleft = par("usr")[1] + 0.9 * (par("usr")[2] - par("usr")[1]),
ybottom = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
xright = par("usr")[1] + 0.95 * (par("usr")[2] - par("usr")[1]),
ytop = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]))
text(x = par("usr")[1] + 0.95 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.3 * (par("usr")[4] - par("usr")[3]),
labels = "Slope (deg)",
adj = c(1,0))
text(x = par("usr")[1] + 0.88 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]),
labels = round(max(df$slope)),
adj = c(1,0.5))
text(x = par("usr")[1] + 0.88 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.15 * (par("usr")[4] - par("usr")[3]),
labels = "0",
adj = c(1,0.5))
text(x = par("usr")[1] + 0.88 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
labels = round(min(df$slope)),
adj = c(1,0.5))
# add legend for density
text(x = par("usr")[1] + 0.75 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.3 * (par("usr")[4] - par("usr")[3]),
labels = "Density",
adj = c(1,0))
text(x = par("usr")[1] + 0.68 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]),
labels = str_pad(round(max(df$density),2), 4, "right", "0"),
adj = c(1,0.5))
text(x = par("usr")[1] + 0.68 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.15 * (par("usr")[4] - par("usr")[3]),
labels = str_pad(round(mean(range(df$density)),2), 4, "right", "0"),
adj = c(1,0.5))
text(x = par("usr")[1] + 0.68 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
labels = str_pad(round(min(df$density),2), 4, "right", "0"),
adj = c(1,0.5))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.25 * (par("usr")[4] - par("usr")[3]),
cex = max(pt.sizes))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.2 * (par("usr")[4] - par("usr")[3]),
cex = mean(c(max(pt.sizes), mean(range(pt.sizes)))))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.15 * (par("usr")[4] - par("usr")[3]),
cex = mean(range(pt.sizes)))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.1 * (par("usr")[4] - par("usr")[3]),
cex = mean(c(min(pt.sizes), mean(range(pt.sizes)))))
points(x = par("usr")[1] + 0.715 * (par("usr")[2] - par("usr")[1]),
y = par("usr")[3] + 0.05 * (par("usr")[4] - par("usr")[3]),
cex = min(pt.sizes))
OK, so I know we’ve seen a lot of these types of plots so far, but this one is the best yet! Really impressive stuff. Now, a word of caution about these results… On one of the iterations of the leave-one-out cross-validation, the model failed to converge. Meaning, it didn’t technically find the “optimal” coefficients to fit. However, you can still use the last set of coefficients the algorithm tested, albeit perhaps violating some statistical principles. Nevertheless, the full model fit just fine, and 43/44 of the iterations did as well, so I’m not too concerned about it.
# define function for adding transparency to colors
add.alpha <- function(col, alpha=0.5){
apply(sapply(col, col2rgb)/255, 2,
function(x)
rgb(x[1], x[2], x[3], alpha=alpha))
}
# define colors
cols.fill <- add.alpha(brewer.pal(7,"Set1"))
cols.line <- brewer.pal(7,"Set1")
# plot density distributions
par(las = 1, mfrow = c(3,3), mar = c(5,5,1,1))
for (i in seq(3,ncol(results.df))){
x <- results.df[,i]
d <- density(x)
plot(d$y ~ d$x, type = "n", xlab = colnames(results.df)[i], ylab = "Density")
polygon(d$y ~ d$x, col = cols.fill[i-2], border = NA)
lines(d$y ~ d$x, col = cols.line[i-2], lwd = 2)
abline(v = median(x), lty = 2, lwd = 2)
if (median(x) < mean(par("usr")[1:2])){
legend("topright", legend = paste0("Med = ", round(median(x),3)),
x.intersp = 0, bty = "n")
} else {
legend("topleft", legend = paste0("Med = ", round(median(x), 3)),
x.intersp = 0, bty = "n")
}
}
Very interesting… So, the new function form changes all of the coefficients quite dramatically. I guess that’s to be expected, but I was sort of surprised to see how much the slope-travel rate coefficients (\(a\) - \(e\)) changed.
# define new travel rate function that incorporates slope, vegetation density, and surface roughness
travelrate.fun <- function(slope,density,roughness,a,b,c,d,e,f,g){
travelrate <- (c * (1 / (pi * b * (1 + ((slope - a) / b) ^ 2))) + d + e * slope)/ (f * density + g * roughness + 1)
return(travelrate)
}
# get coefficients (medians from LOOCV)
a <- median(results.df$a)
b <- median(results.df$b)
c <- median(results.df$c)
d <- median(results.df$d)
e <- median(results.df$e)
f <- median(results.df$f)
g <- median(results.df$g)
# get range of travel rate predictions
min.rate <- 100
max.rate <- -100
for (i in seq(0,1,0.1)){
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
if (min(rates) < min.rate) min.rate <- min(rates)
if (max(rates) > max.rate) max.rate <- max(rates)
}
# plot out the results
cols <- viridis(11)
counter <- 0
par(las = 1, mar = c(5,5,1,8))
plot(x = c(-45,45), y = c(min.rate, max.rate), type = "n",
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
for (i in seq(0,1,0.1)){
counter <- counter + 1
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
lines(rates ~ slopes, lwd = 2, col = cols[counter])
}
par(xpd = NA)
legend(x = par("usr")[2], y = mean(par("usr")[3:4]),
legend = seq(0,1,0.1), title = "Veg Dens", lwd = 2, col = cols,
bty = "n", yjust = 0.5)
Whoa! Cool! So, even though we’re using the asymmetric function, this multiplicative density-roughness function form has removed the weird steep downhill travel rate uptick.
OK, so that was the asymmetric function… Since it resolved that downhill travel problem, I’m not sure it’s worth exploring how this multiplicative function form applies to the symmetric function form (Approach #1 above). But I am interested to see how it performs with the hybrid AllTrails-Levan model (Approach 2 above).
# define new travel rate function that incorporates slope, vegetation density, and surface roughness
travelrate.fun <- function(slope,density,roughness,a,b,c,d,e,f,g){
travelrate <- (c * (1 / (pi * b * (1 + ((slope - a) / b) ^ 2))) + d + e * slope) / (f * density + g * roughness + 1)
return(travelrate)
}
# define coefficients
a <- slope.rate.terms$a[slope.rate.terms$percentile == 0.5]
b <- slope.rate.terms$b[slope.rate.terms$percentile == 0.5]
c <- slope.rate.terms$c[slope.rate.terms$percentile == 0.5]
d <- slope.rate.terms$d[slope.rate.terms$percentile == 0.5]
e <- slope.rate.terms$e[slope.rate.terms$percentile == 0.5]
f <- median(results.df$f)
g <- median(results.df$g)
# get range of travel rate predictions
min.rate <- 100
max.rate <- -100
for (i in seq(0,1,0.1)){
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
if (min(rates) < min.rate) min.rate <- min(rates)
if (max(rates) > max.rate) max.rate <- max(rates)
}
# plot out the results
cols <- viridis(11)
counter <- 0
par(las = 1, mar = c(5,5,1,8))
plot(x = c(-45,45), y = c(min.rate, max.rate), type = "n",
xlab = "Slope (deg)", ylab = "Travel Rate (m/s)")
grid()
for (i in seq(0,1,0.1)){
counter <- counter + 1
slopes <- seq(-45,45,0.1)
densities <- rep(i, length(slopes))
roughnesses <- rep(0, length(slopes))
rates <- travelrate.fun(slopes, densities, roughnesses, a, b, c, d, e, f, g)
lines(rates ~ slopes, lwd = 2, col = cols[counter])
}
par(xpd = NA)
legend(x = par("usr")[2], y = mean(par("usr")[3:4]),
legend = seq(0,1,0.1), title = "Veg Dens", lwd = 2, col = cols,
bty = "n", yjust = 0.5)
Fascinating. In contrast to the linear approach, which resulted in some pretty restrictive thresholds above which movement (positive travel rates) was possible. The multiplicative approach is probably a bit too permissive. For example, is it really reasonable to expect that someone could travel ~0.5 m/s through 100% vegetation density? I really don’t think so. But, is it better to present a set of equations that doesn’t result in the somewhat hard-to-explain negative travel rates? Or is it better to present a set of equations that “look” better, though perhaps provide an overly-optimistic view of realistic travel in a wildland environment? Something to keep in mind is that the ground surface roughness has been assumed to be zero throughout all of these simulations, so that would certainly slow things down a bit too.
New Conclusions (6/22/2023)
There are seemingly four paths forward:
- Linear Density-Roughness Effects 1a. Using Levan-driven slope-symmetric function 1b. Using Levan-AllTrails hybrid function
- Multiplicative Density-Roughness Effects 2a. Using Levan-driven slope-asymmetric function 2b. Using Levan-AllTrails hybrid function
My gut says either 1b or 2b, with a slight preference for the latter. But, if you’ve made it this far, I’m interested to get your thoughts as well.
Phew… That was a lot.