Background and Research Questions

The objective of this study is to analyze the tradeoffs between visibility and travel time when generating a least cost path which incorporates a visibility measure.

We have generated landscape rasters from lidar, and LCPs as escape routes based on the landscape rasters and travel rate equations. Now we attempt to analyze the LCPs generated and pursue our objective. In this document, we look at two main research questions:

  1. If I choose an escape route that helps me maintain high visibility, how much more time is it going to take me to get there?
  2. If I choose an escape route that gets me there fastest, how much more exposed to low-visibility environments am I?

These research questions will guide this analysis and provide a basis for additional questions related to our objective.


Set the working directory, packages, and read in data

setwd(wd)
library(sf)
library(terra)
library(dplyr)
library(RColorBrewer)
library(viridis)
library(colorspace)
library(tidyr)
library(ggplot2)
library(gridExtra)

The LCP generation script produces a .csv of attributes and calculated metrics, and a shapefile of all the lines generated with the .csv data as attributes per line. We will read in this data as well as the study area DTM for visualization.

lcps = read.csv(csv_data)
lcp_lines = read_sf(shp_data)
dtm = rast(dtm)


Research Question 1


Question 1 is:

If I choose an escape route that helps me maintain high visibility, how much more time is it going to take me to get there?

A sample answer could be:

“It’s going to take me, on average, in this study area, 2 times longer to arrive at my destination while maintaining high visibility versus choosing the fastest route.”

This section attempts to answer this question. Here are our preliminary steps to analyze this question:

  1. Temporarily ignore rescaling methods

  2. For each start-end point pair, get the time it takes to arrive at destination with the fastest escape route (i.e., the one where visibility bears no weight)

  3. Compare that time to each of the other increasingly-visibility-weighted LCPs

    • assumed output: with increasing visibility weight, increasing time
    • thing we don’t know but need to know: how much does that time increase


Based on our question and sample answer, we are looking to sort through our data to find the following characteristics:

  • visibility weight
  • travel time
  • location

And the analysis we have to perform is going to answer how much of these characteristics are changing as we change scenarios.


First we will look at all of our LCPs generated as points, plotted according to path length and travel time.

lcps %>%
  ggplot(aes(x = pathLength, y = traveltime)) +
  geom_point() +
  labs(x = "Path Length (m)", y = "Travel Time (s)", title = "Path length to travel time; all LCPs")

That’s a lot of points!
Looking into visualizations for each of the steps above, let’s start with ignoring rescaling methods.

lcps_no_rescale = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = pathLength, y = traveltime)) +
  geom_point() +
  labs(x = "Path Length (m)", y = "Travel Time (s)", title = "Path Length to travel time where there is no rescaled visibility")
lcps_no_rescale

This operation performs the following:

  • selects LCPs table as data
  • filters the rescaling factor to equal 0, meaning no rescaled LCPs are included in the table
  • filters out a single v0 value so that we only have one LCP per point pair and weight value
    • it doesn’t matter which, because they are all evaluated at k = 0, i.e. not rescaled
  • generate a plot of travel time and path length to visualize how many LCPs are then used for analysis

Notice in comparing the two scatter plots produced above that the x and y axes are different. We may be able to assume that as we remove the visibility rescaled LCPs, the path length and travel time values are lower; this simple observation follows our sample hypothesis.


Next, let’s filter our data to find the fastest escape route for each start/end point pair, which would be the route where visibility = 0.

lcps_w1_1 = lcps %>%
  group_by(pt1x) %>%
  filter(k == 0, v0 == 1, w1 == 1) %>%
  tibble::rownames_to_column(var="name") %>%
  ggplot(aes(x = name, y = traveltime)) +
  geom_point() +
  labs(x = "LCP Point Pair", y = "Travel Time (s)", title = "Fastest travel time per point pair") +
  scale_x_discrete(breaks = seq(0, 100, by = 10))
lcps_w1_1

This operation performs the following:

  • selects LCPs table as data
  • groups the data by the starting point x-coordinate
    • this means the point pairs are grouped, considering they all have discrete start points
  • filter out rescaled LCPs
  • filter points to find the LCP where w1 (speed) is 1, so visibility has no weight in the LCP generation
  • renames pt1x to ascending numerical value for plotting purposes
  • plot renamed point pair as x, and fastest escape route as y

Finally, compare the other visibility-weighted LCPs travel times.

lcps_min_travel = lcps %>%
  filter(k == 0, v0 == 1, w1 == 1) %>%
  select(pt1x, traveltime)

lcps = left_join(lcps, lcps_min_travel, by = "pt1x")

lcps = mutate(lcps, traveltime_mult = traveltime.x / traveltime.y)

tt_mult_plot = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = (as.factor(1-w1)), y = traveltime_mult)) +
  geom_boxplot() +
  labs(x = "Visibility Weights", y = "Travel Time Multiplier")
tt_mult_plot

This operation performs the following:

  • creates a table of just minimum travel time specified by 0 visibility weight, and pt1x as the identifier
  • join the minimum travel time values to the main data table with pt1x as the key
  • calculate a travel time multiplier as the actual travel time over the minimum travel time for each row of our data set
  • plot the results with visibility weights as our x value, and the travel time multiplier as y.

Our results match our expected hypothesis; as we increase the weight or importance of visibility for generating our LCPs, travel time increases. In this relationship, it appears to have an exponential effect.

According to Phil, either;

  1. Visibility is already moderately good on least cost paths (i.e. the path doesn’t change much as the weight of visibility is increased), or
  2. Visibility can be substantially improved with minor impacts on least cost path length (increasing the weight of visibility only makes the LCP a little longer until you get to very high weights).


2/13/2024 –

After yesterday’s meeting, we want to look at a “spaghetti plot” to show raw travel rate values rather than a multiplier. Hopefully what we will see are the absolute values of travel times which aims to provide clarity to the relationship between travel rate and visibility weight.

This plot should have:

  • x-axis: visibility weight
  • y-axis: travel time
  • 100 lines: one line per point pair
spaghetti1 = lcps %>%
  filter(k == 0, v0 == 0) %>%
  ggplot(aes(as.factor(1-w1), traveltime.x, group = pt1x, color = pt1x)) +
  geom_point(show.legend = F) +
  geom_line(show.legend = F) +
  labs(x = "Visibility Weights", y="Travel Time", title="Travel time value \n as visibility weight increases, \n per point pair", color = "Point Pair") +
  theme(legend.position = "left", legend.key.size = unit(0.1, 'cm'), legend.box="horizontal", plot.title=element_text(size=10, hjust=0.5), axis.title.x=element_text(size=9),axis.title.y=element_text(size=9)) +
  scale_color_distiller(palette = "Spectral")
spaghetti1

Now I’ll look at the two visualizations next to each other:

grid.arrange(spaghetti1, tt_mult_plot, ncol = 2)

Extra plot with the travel time multiplier as lines instead of box plots for funsies:

spaghetti_tt_mult = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = as.factor(1-w1), y = traveltime_mult, group = pt1x, color=pt1x)) +
  geom_point(show.legend = F) +
  geom_line(show.legend = F) +
  labs(x = "Visibility Weights", y="Travel Time Multiplier", title = "Travel time multiplier \n as visibility weight increases, \n per point pair", color="Point Pair") + theme(legend.position = "right", legend.key.size = unit(0.1, 'cm'), plot.title=element_text(size=10,hjust=0.5), axis.title.x=element_text(size=9),axis.title.y=element_text(size=9)) +
  scale_color_distiller(palette = "Spectral")
grid.arrange(spaghetti1, spaghetti_tt_mult, ncol = 2)


Research Question 2


Question 2 is:

If I choose an escape route that gets me there fastest, how much more exposed to low-visibility environments am I?

A sample answer could be:

“When I choose the fastest path, my average visibility decreases by half, compared to equally weighting visibility and travel rate.”

So the two variables we are going to be plotting here are travel time and average visibility. Steps to answer this question are as follows:

  1. Temporarily ignore rescaling methods

  2. For each start-end point pair, get the average visibility with the most-visibility-optimized route (i.e., the one where speed bears no weight)

  3. Compare this (maximum) average visibility to other decreasingly-visibility-weighted LCPs

  4. Compare percent high visibility to other decreasingly-visibility-weighted LCPs

Notice these steps are very similar to our first question. Basically, where question 1 is: higher visibility, how much longer?, question 2 is: faster, how much less visible?


Okay anyways, let’s look at some plots of the average visibility of our LCP data. See my second plot’s code in this document for how we set up the plots displayed below.

First, we’ll compare average visibility to travel time and path length.

Not a lot of relationship to see there…


Let’s also look at all three variables together, with average visibility as our point colors

And here’s the same plot with speed weights as our point sizes


Next we will find the maximum average visibility value per point pair, and plot those results.

lcps_w1_0 = lcps %>%
  group_by(pt1x) %>%
  filter(k == 0, v0 == 1, w1 == 0) %>%
  tibble::rownames_to_column(var='name') %>%
  ggplot(aes(x=name,y=avgVis)) +
  geom_point() +
  labs(x = "LCP Point Pair", y = "Average Visibility", title = "Maximum average visibility per point pair") +
  scale_x_discrete(breaks = seq(0,100,by=10))
lcps_w1_0


Now we can calculate the average visibility multiplier, and compare this value to our speed weights.

lcps_max_vis = lcps %>%
  filter(k == 0, v0 == 1, w1 == 0) %>%
  select(pt1x, avgVis)

lcps = left_join(lcps, lcps_max_vis, by='pt1x')

lcps = mutate(lcps, avgvis_mult = avgVis.x/avgVis.y)

avg_vis_mult_plot = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = (as.factor(w1)), y = avgvis_mult)) +
  geom_boxplot() +
  labs(x = "Speed Weight", y = "Average Visibility Multiplier")
avg_vis_mult_plot


2/13/2024 –

Again, repeated spaghetti plots for average visibility visualizations

spaghetti2 = lcps %>%
  filter(k == 0, v0 == 0) %>%
  ggplot(aes(as.factor(w1), avgVis.x, group = pt1x, color = pt1x)) +
  geom_point(show.legend = F) +
  geom_line(show.legend = F) +
  labs(x = "Speed Weights", y="Average Visibility", title="Average visibility value \n as speed weight increases, \n per point pair", color = "Point Pair") +
  theme(plot.title=element_text(size=10, hjust=0.5), axis.title.x=element_text(size=9),axis.title.y=element_text(size=9)) +
  scale_color_distiller(palette = "Paired")
spaghetti2

Side-by-side with box plot, then side-by-side with average visibility multiplier spaghetti plot

grid.arrange(spaghetti2, avg_vis_mult_plot, ncol = 2)

spaghetti_avg_vis_mult = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = as.factor(w1), y = avgvis_mult, group = pt1x, color=pt1x)) +
  geom_point(show.legend = F) +
  geom_line(show.legend = F) +
  labs(x = "Speed Weights", y="Average Visibility Multiplier", title = "Average visibility multiplier \n as speed weight increases, \n per point pair", color="Point Pair") + theme(plot.title=element_text(size=10,hjust=0.5), axis.title.x=element_text(size=9),axis.title.y=element_text(size=9)) +
  scale_color_distiller(palette = "Paired")
grid.arrange(spaghetti1, spaghetti_avg_vis_mult, ncol = 2)



Moving on from average visibility to our metric of percent high visibility per LCP, we will create the same plots as the ones we made for average visibility. The percent high visibility metric is evaluated as the percent of LCPs with visibility values greater than 0.25. Our percent high visibility multiplier will be computed the same way as our previous two analyses.

First, I’ll repeat the path length/travel time scatter plot, with percent high visibility as the point color

length_time_highvis = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = pathLength, y=traveltime.x)) +
  geom_point(aes(color = percentVisT25)) +
  labs(x = "Path Length (m)", y = "Travel Time (s)", title = "Travel time and path length with no rescaled LCPs") +
  scale_color_gradient(low = "darkgreen", high = "darkseagreen1")
length_time_highvis

Next we will create the box plots with the high visibility multiplier and speed weights.

lcps_max_highvis = lcps %>%
  filter(k == 0, v0 == 1, w1 == 0) %>%
  select(pt1x, percentVisT25)

lcps = left_join(lcps, lcps_max_highvis, by='pt1x')

lcps = mutate(lcps, highvis_mult = percentVisT25.x/percentVisT25.y)

high_vis_mult_plot = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = (as.factor(w1)), y = highvis_mult)) +
  geom_boxplot() +
  labs(x = "Speed Weight", y = "Percent High Visibility Multiplier")
high_vis_mult_plot
## Warning: Removed 407 rows containing non-finite values (`stat_boxplot()`).

2/13/2024 –

One more time, repeated spaghetti plots for high visibility visualizations

spaghetti3 = lcps %>%
  filter(k == 0, v0 == 0) %>%
  ggplot(aes(as.factor(w1), percentVisT25.x, group = pt1x, color = pt1x)) +
  geom_point(show.legend = F) +
  geom_line(show.legend = F) +
  labs(x = "Speed Weights", y="Percent High Visibility (>0.25)", title="Percent High visibility value \n as speed weight increases, \n per point pair", color = "Point Pair") +
  theme(plot.title=element_text(size=10, hjust=0.5), axis.title.x=element_text(size=9),axis.title.y=element_text(size=9)) +
  scale_color_distiller(palette = "Paired")
spaghetti3

Side-by-side with box plot, then side-by-side with high visibility multiplier spaghetti plot

grid.arrange(spaghetti3, high_vis_mult_plot, ncol = 2)
## Warning: Removed 407 rows containing non-finite values (`stat_boxplot()`).

spaghetti_high_vis_mult = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = as.factor(w1), y = highvis_mult, group = pt1x, color=pt1x)) +
  geom_point(show.legend = F) +
  geom_line(show.legend = F) +
  labs(x = "Speed Weights", y="Percent High Visibility Multiplier", title = "High visibility multiplier \n as speed weight increases, \n per point pair", color="Point Pair") + theme(plot.title=element_text(size=10,hjust=0.5), axis.title.x=element_text(size=9),axis.title.y=element_text(size=9)) +
  scale_color_distiller(palette = "Paired")
grid.arrange(spaghetti1, spaghetti_high_vis_mult, ncol = 2)
## Warning: Removed 406 rows containing missing values (`geom_point()`).
## Warning: Removed 406 rows containing missing values (`geom_line()`).


Visibility thresholds

After generating new data, we are looking at different thresholds of visibility with the same analysis of the high visibility assessment above. Other thresholds were 5, 10, and 50.


THRESHOLD = 5

length_time_visT5 = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = pathLength, y=traveltime.x)) +
  geom_point(aes(color = percentVisT5)) +
  labs(x = "Path Length (m)", y = "Travel Time (s)", title = "Travel time and path length with no rescaled LCPs") +
  scale_color_gradient(low = "darkgreen", high = "darkseagreen1")
length_time_visT5

lcps_visT5 = lcps %>%
  filter(k == 0, v0 == 1, w1 == 0) %>%
  select(pt1x, percentVisT5)

lcps = left_join(lcps, lcps_visT5, by='pt1x')

lcps = mutate(lcps, vis_t5_mult = percentVisT5.x/percentVisT5.y)

lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = (as.factor(w1)), y = vis_t5_mult)) +
  geom_boxplot() +
  labs(x = "Speed Weight", y = "Vis threshold = 5")


THRESHOLD = 10

length_time_visT10 = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = pathLength, y=traveltime.x)) +
  geom_point(aes(color = percentVisT10)) +
  labs(x = "Path Length (m)", y = "Travel Time (s)", title = "Travel time and path length with no rescaled LCPs") +
  scale_color_gradient(low = "darkgreen", high = "darkseagreen1")
length_time_visT10

lcps_visT10 = lcps %>%
  filter(k == 0, v0 == 1, w1 == 0) %>%
  select(pt1x, percentVisT10)

lcps = left_join(lcps, lcps_visT10, by='pt1x')

lcps = mutate(lcps, vis_t10_mult = percentVisT10.x/percentVisT10.y)

lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = (as.factor(w1)), y = vis_t10_mult)) +
  geom_boxplot() +
  labs(x = "Speed Weight", y = "Vis threshold = 10")
## Warning: Removed 44 rows containing non-finite values (`stat_boxplot()`).


THRESHOLD = 25

Note this is the same as high visibility above

length_time_highvis

high_vis_mult_plot
## Warning: Removed 407 rows containing non-finite values (`stat_boxplot()`).


THRESHOLD = 50

length_time_visT50 = lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = pathLength, y=traveltime.x)) +
  geom_point(aes(color = percentVisT50)) +
  labs(x = "Path Length (m)", y = "Travel Time (s)", title = "Travel time and path length with no rescaled LCPs") +
  scale_color_gradient(low = "darkgreen", high = "darkseagreen1")
length_time_visT50

lcps_visT50 = lcps %>%
  filter(k == 0, v0 == 1, w1 == 0) %>%
  select(pt1x, percentVisT50)

lcps = left_join(lcps, lcps_visT50, by='pt1x')

lcps = mutate(lcps, vis_t50_mult = percentVisT50.x/percentVisT50.y)

lcps %>%
  filter(k == 0, v0 == 1) %>%
  ggplot(aes(x = (as.factor(w1)), y = vis_t50_mult)) +
  geom_boxplot() +
  labs(x = "Speed Weight", y = "Vis threshold = 50")
## Warning: Removed 924 rows containing non-finite values (`stat_boxplot()`).


Rescaled Visibility

Applying above analyses to rescaled visibility LCPs

Question 1.

v0_labs = c("Logarithmic", "Logistic Growth", "Exponential")
names(v0_labs) = c("0", "0.5", "1")
k_labs = c("No Rescale", "k = 5", "k = 10", "k = 20")
names(k_labs) = c("0", "5", "10", "20")

tt_mult_plot_rescale = lcps %>%
  ggplot(aes(x = (as.factor(1-w1)), y = traveltime_mult)) +
  geom_boxplot() +
  labs(x = "Visibility Weights", y = "Travel Time Multiplier", title = "Travel time multiplier per Visibility weight, \n rescaled visibility variables") +
  facet_grid(v0 ~ k, labeller = labeller(v0 = v0_labs, k = k_labs)) +
  theme(plot.title=element_text(size=12,hjust=0.5))

tt_mult_plot_rescale

Question 2.

avg_vis_mult_plot_rescale = lcps %>%
  ggplot(aes(x = (as.factor(w1)), y = avgvis_mult)) +
  geom_boxplot() +
  labs(x = "Speed Weights", y = "Average Visibility Multiplier", title = "Average visibility multiplier per speed weight, \n rescaled visibility variables") +
  facet_grid(v0 ~ k, labeller = labeller(v0 = v0_labs, k = k_labs)) +
  theme(plot.title=element_text(size=12,hjust=0.5))

avg_vis_mult_plot_rescale

high_vis_mult_plot_rescale = lcps %>%
  ggplot(aes(x = (as.factor(w1)), y = highvis_mult)) +
  geom_boxplot() +
  labs(x = "Speed Weights", y = "Percent High Visibility Multiplier", title = "Percent High Visibility (>0.25) multiplier per speed weight, \n rescaled visibility variables") +
  facet_grid(v0 ~ k, labeller = labeller(v0 = v0_labs, k = k_labs)) +
  theme(plot.title=element_text(size=12,hjust=0.5))

high_vis_mult_plot_rescale
## Warning: Removed 4884 rows containing non-finite values (`stat_boxplot()`).