library(tidyverse)
options(dplyr.summarise.inform = FALSE)
library(tidycensus)
library(leaflet)
library(tidyquant)
library(extrafont)
library(fiftystater)
library(readxl)

epsg2163 <- leafletCRS(
  crsClass = "L.Proj.CRS",
  code = "EPSG:2163",
  proj4def = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs",
  resolutions = 2^(16:7))

metro_counties <- c("Hennepin", "Ramsey", "Dakota", "Carver", "Scott", "washington", "Anoka")

streetlight <- read_csv("./VMT Analysis/Data/streetlight_data_all.csv") %>%
  filter(state== "Minnesota") %>%
  mutate(date= paste(year, month, day, sep = "-"),
         date= as.Date(date),
         region = case_when(
           county %in% metro_counties ~ "Metro",
           TRUE ~ "Greater MN"
         )) %>%
  select(-c("year", "month", "day"))

sl_MNtotal <- read_csv("./VMT Analysis/Data/streetlight_data_mn.csv") %>%
  mutate(date= paste(year, month, day, sep = "-"),
         date= as.Date(date))

COVID-19 has had an impact on when and how Minnesotans travel due to voluntary social distancing, state, and changes in employment. It is important to understand changes in vehicle miles traveled (VMT) in the state for many reasons including as an indicator for social distancing, its impact on state revenues generated through vehicle travel, and its impact on the condition of transportation infrastructure. To estimate changes in VMT across the state, this paper uses data generated by Streetlight Data. A thorough explanation of the methodology used to estimate VMT changes is available through this white paper.

Figure 1 shows the statewide trend in VMT from between March 1 and October 16, 2020. VMT declined drastically in late April, and has slowly recovered afterward. By October 16, VMT was just above 80 percent of what it had been at the beginning of March. VMT has declined declined since the beginning of August, before which it had consistently increased since mid-April. There could be several explanations for the end to the consistent VMT gains made in previous months. It may be that the remaining decline is largely due to those working from home in jobs that will be remote in the long-term. It could also track with the general increase in COVID-19 cases over this time frame. The data shows a consistent rise and fall pattern related to the difference in travel patterns throughout the week. A moving average is used to smooth out the variation in the trends caused by these differences in travel patterns.

#vmt change in MN with important events
#orders from https://mn.gov/covid19/for-minnesotans/stay-safe-mn/stay-safe-plan.jsp

ggplot(sl_MNtotal, aes(x=date, y= vmt/1000000))+
  geom_bar(stat="identity",color= "white", fill="#1f78b4", alpha=0.3)+
  geom_ma(ma_fun=SMA, n=7, color="#1f78b4", size=2, linetype="solid")+
  scale_y_continuous(breaks= c(0, 100, 200, 300, 400), expand= c(0,2))+
  scale_x_date(expand= c(0,1))+
  geom_segment(aes(x=as.Date(c("2020-03-05")), y=242, xend=as.Date(c("2020-03-05")), yend=260), size=1.25, 
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-03-26")), y=69, xend=as.Date(c("2020-03-26")), yend=130), size=1.25,
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-05-18")), y=97, xend=as.Date(c("2020-05-18")), yend=115), size=1.25,
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-06-01")), y=120, xend=as.Date(c("2020-06-01")), yend=138), size=1.25,
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-06-10")), y=136, xend=as.Date(c("2020-06-10")), yend=151), size=1.25,
               color="black")+
  annotate("label", x=as.Date(sl_MNtotal$date[10]), y=274, label="1st Case", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  annotate("label", x=as.Date(sl_MNtotal$date[27]), y=115, label="Stay-at-Home \nOrder", family= "Roboto Condensed",
           size=5, fill='white', alpha= 0.1, label.size= NA, hjust=0)+
  annotate("label", x=as.Date(sl_MNtotal$date[73]), y=128, label="Phase 1", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  annotate("label", x=as.Date(sl_MNtotal$date[87]), y=150, label="Phase 2", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  annotate("label", x=as.Date(sl_MNtotal$date[102]), y=168, label="Phase 3", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  #scale_x_date(limits = as.Date(c(NA, "2020-07-03")))+
  theme_classic()+
  labs(x="", y="VMT (Millions)",
       title= "Figure 1. Minnesota VMT During COVID-19",
       caption = "Minnesota's reopening has been broken into phases. Phase 1 ended limits on non-essential \ntravel and allowed some business reopenings. Following phases lifted further restrictions.")+
  theme(text= element_text(family= "Roboto Condensed"),
        panel.grid.major = element_line(colour="light grey"),
        plot.title = element_text(size = 18),
        plot.caption = element_text(size=11),
        panel.grid.major.x = element_blank(),
        axis.text.y = element_text(size=18, color = "black"),
        axis.text.x = element_text(size= 18, color = "black"),
        axis.title.y= element_text(size=20, color="black", margin = margin(r=10)))

This figure also shows the dates of the Governor’s stay-at-home order and reopening phases 1-3. Phase 4 is planned to be the final phase of reopening. Figure 1 does not provide clear evidence that the orders impacted travel behavior, although a more sophisticated analysis could be used to explore this further. VMT dropped significantly prior to the Stay-at-Home order, giving evidence that the order itself was not the primary cause of the decline. After flattening out daily variation, this figure provides some visual evidence that the recovery rate of daily VMT may have briefly increased after reopening Phase 1 began . This effect would be understandable because it was the first period travel was allowed for non-essential activities. However, the increase in the VMT recovery only lasted for about a week.

In addition to smoothing out weekly variation, it’s very important to seasonally adjust VMT data because average VMT changes significantly by month.

adjustments <- read_csv("./VMT Analysis/Data/Raw Data/seasonal_adjustments.csv")

test2 <-sl_MNtotal %>%
  inner_join(adjustments, by="month") %>%
  mutate(vmt_adj= vmt * avg_adjust) %>%
  ggplot(aes(x=date, y= vmt/1000000))+
  geom_bar(stat="identity",color= "white", fill="#1f78b4", alpha=0.3)+
  geom_ma(ma_fun=SMA, n=7, color="#1f78b4", size=2, linetype="solid")+
  scale_y_continuous(breaks= c(0, 100, 200, 300, 400), expand= c(0,2))+
  scale_x_date(expand= c(0,1))+
  geom_segment(aes(x=as.Date(c("2020-03-05")), y=242, xend=as.Date(c("2020-03-05")), yend=260), size=1.25, 
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-03-26")), y=69, xend=as.Date(c("2020-03-26")), yend=130), size=1.25,
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-05-18")), y=97, xend=as.Date(c("2020-05-18")), yend=115), size=1.25,
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-06-01")), y=120, xend=as.Date(c("2020-06-01")), yend=138), size=1.25,
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-06-10")), y=136, xend=as.Date(c("2020-06-10")), yend=151), size=1.25,
               color="black")+
  annotate("label", x=as.Date(sl_MNtotal$date[10]), y=274, label="1st Case", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  annotate("label", x=as.Date(sl_MNtotal$date[27]), y=115, label="Stay-at-Home \nOrder", family= "Roboto Condensed",
           size=5, fill='white', alpha= 0.1, label.size= NA, hjust=0)+
  annotate("label", x=as.Date(sl_MNtotal$date[73]), y=128, label="Phase 1", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  annotate("label", x=as.Date(sl_MNtotal$date[87]), y=150, label="Phase 2", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  annotate("label", x=as.Date(sl_MNtotal$date[102]), y=168, label="Phase 3", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  #scale_x_date(limits = as.Date(c(NA, "2020-07-03")))+
  theme_classic()+
  labs(x="", y="VMT (Millions)",
       title= "Figure 1. Minnesota VMT During COVID-19",
       caption = "Minnesota's reopening has been broken into phases. Phase 1 ended limits on non-essential \ntravel and allowed some business reopenings. Following phases lifted further restrictions.")+
  theme(text= element_text(family= "Roboto Condensed"),
        panel.grid.major = element_line(colour="light grey"),
        plot.title = element_text(size = 18),
        plot.caption = element_text(size=11),
        panel.grid.major.x = element_blank(),
        axis.text.y = element_text(size=18, color = "black"),
        axis.text.x = element_text(size= 18, color = "black"),
        axis.title.y= element_text(size=20, color="black", margin = margin(r=10)))

Figure 2 shows the percent change of Minnesota’s VMT during the pandemic using average VMT in the first week of March as the baseline. It also shows the national change in VMT. Minnesota dropped slightly more than the nation as a whole, to about 80 percent less than the baseline, and has recovered at a similar rate.

change_pal <- c("#1f78b4","#b2df8a")

change_both <- read_csv("./VMT Analysis/Data/streetlight_data_all.csv") %>%
  mutate(date= paste(year, month, day, sep = "-"),
         date= as.Date(date)) %>%
  select(-c("year", "month", "day")) %>%
  mutate(natl_base= sum(vmt[date == "2020-03-01" | date == "2020-03-02" |date == "2020-03-03" |date == "2020-03-04" |date == "2020-03-05" |date == "2020-03-06" |date == "2020-03-07"])/7,
        mn_base= 246645371) %>% #calculated above
  group_by(date) %>%
  summarise(vmt_natl= sum(vmt),
            natl_base= mean(natl_base),
            mn_base= mean(mn_base),
            natl_change= (vmt_natl/natl_base-1)*100) %>%
  inner_join(sl_MNtotal, by= "date") %>%
  mutate(mn_change= (vmt/mn_base-1)*100)

change_mn <- change_both %>%
  select(date, mn_base, mn_change) %>%
  mutate(area= "Minnesota") %>%
  rename(
    base= mn_base,
    perc_change= mn_change)

change_natl <- change_both %>%
  select(date, natl_base, natl_change) %>%
  mutate(area= "U.S.") %>%
  rename(
    base= natl_base,
    perc_change= natl_change)

change_final <- rbind(change_mn, change_natl)

ggplot(change_final, aes(x=date, y= perc_change, color= area))+
  geom_ma(ma_fun=SMA, n=7, size=2, linetype="solid")+
  scale_color_manual(values = c(change_pal))+
  scale_x_date(limits = as.Date(c(NA, "2020-10-16")))+
  #ylim(0, 150)+
  theme_classic()+
  labs(
    x="",
    y="% Change in VMT",
    color= "",
    title = "Figure 2. Minnesota vs. National VMT Change",
    subtitle = "") +
  theme(text= element_text(family= "Roboto Condensed"),
        plot.title = element_text(size=18),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(colour="light grey"),
        panel.grid.major.x = element_blank(),
        axis.text.y = element_text(size=18, color = "black"),
        axis.text.x = element_text(size= 18, color = "black"),
        axis.title.y= element_text(size=18, color="black"),
        legend.position = "bottom",
        legend.margin = margin(t= -15),
        legend.text = element_text(size=18))

It is also important to explore how the pattern of VMT change varies across the state. One key divide in Minnesota is between the Twin Cities and the remainder of the state. Figure 3 shows the trend in VMT for both sections of the state. The figure shows that VMT fell slightly less outside of the metro area, and has recovered at a faster rate until as well. This difference has become larger during the fall, despite the fact that new COVID-19 infections in Minnesota are disproportionately occurring outside of the metro area.

metro2 <- c("Dakota", "Carver", "Scott", "washington", "Anoka")

sl_metro <- streetlight %>%
  group_by(date, region) %>%
  summarise(vmt_total= sum(vmt))

ggplot(sl_metro, aes(x=date, y= vmt_total/1000000, color=region))+
  geom_ma(ma_fun=SMA, n=7, size=2, linetype="solid")+
  scale_color_manual(values = c(change_pal))+
  scale_x_date(limits = as.Date(c(NA, "2020-10-16")))+
  geom_segment(aes(x=as.Date(c("2020-03-05")), y=115, xend=as.Date(c("2020-03-05")), yend=125), size=1.25, 
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-03-26")), y=39, xend=as.Date(c("2020-03-26")), yend=60), size=1.25,
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-05-18")), y=58, xend=as.Date(c("2020-05-18")), yend=68), size=1.25,
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-06-01")), y=71, xend=as.Date(c("2020-06-01")), yend=81), size=1.25,
               color="black")+
  geom_segment(aes(x=as.Date(c("2020-06-10")), y=79, xend=as.Date(c("2020-06-10")), yend=89), size=1.25,
               color="black")+
  annotate("label", x=as.Date(sl_MNtotal$date[4]), y=132, label="1st Case", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  annotate("label", x=as.Date(sl_MNtotal$date[26]), y=57, label="Stay-at-Home \nOrder", family= "Roboto Condensed",
           size=5, fill='white', alpha= 0.1, label.size= NA, hjust=0)+
  annotate("label", x=as.Date(sl_MNtotal$date[75]), y=73, label="Phase 1", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  annotate("label", x=as.Date(sl_MNtotal$date[89]), y=86, label="Phase 2", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  annotate("label", x=as.Date(sl_MNtotal$date[102]), y=95, label="Phase 3", family= "Roboto Condensed", size=5,
           fill='white', alpha= 0.4, label.size= NA)+
  ylim(0, 150)+
  theme_classic()+
  labs(
    x="",
    y="VMT (Millions)",
    color= "",
    title = "Figure 3. Regional VMT Change",
    subtitle = "")+
  theme(text= element_text(family= "Roboto Condensed"),
        plot.title = element_text(size=18),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(colour="light grey"),
        panel.grid.major.x = element_blank(),
        axis.text.y = element_text(size=18, color = "black"),
        axis.text.x = element_text(size= 18, color = "black"),
        axis.title.y= element_text(size=18, color="black"),
        legend.position = "bottom",
        legend.margin = margin(t= -15),
        legend.text = element_text(size=18))

The change in VMT can also be analyzed at the county level. Figure 4 shows how the pandemic impacted VMT across the state on April 1st and October 15th. The former date was the tail end of the period where statewide VMT dropped the most, while October 15th is very recent. By mid-October, most counties were within 30 percentage points or above their VMT from the first week of March. Counties in the metro area were among those with the largest continuing declines.

counties_NA <- streetlight %>% #find a way to make this more elegant
  add_row(county_id= 27077, county= "Lake of the Woods", date= as.Date(c("2020-03-01"))) %>%
  add_row(county_id= 27077, county= "Lake of the Woods", date= as.Date(c("2020-04-01"))) %>%
  add_row(county_id= 27077, county= "Lake of the Woods", date= as.Date(c("2020-10-15"))) %>% 
  add_row(county_id= 27031, county= "Cook", date= as.Date(c("2020-03-01"))) %>%
  add_row(county_id= 27031, county= "Cook", date= as.Date(c("2020-04-01"))) %>%
  add_row(county_id= 27031, county= "Cook", date= as.Date(c("2020-10-15"))) %>%
  add_row(county_id= 27155, county= "Traverse", date= as.Date(c("2020-03-01"))) %>%
  add_row(county_id= 27155, county= "Traverse", date= as.Date(c("2020-04-01"))) %>%
  add_row(county_id= 27155, county= "Traverse", date= as.Date(c("2020-10-15"))) %>%
  add_row(county_id= 27087, county= "Mahnomen", date= as.Date(c("2020-03-01"))) %>%
  add_row(county_id= 27087, county= "Mahnomen", date= as.Date(c("2020-04-01"))) %>%
  add_row(county_id= 27087, county= "Mahnomen", date= as.Date(c("2020-10-15")))

sl_counties <- counties %>%
  mutate(county_id= as.numeric(GEOID)) %>%
  select(county_id, geometry) %>%
  full_join(counties_NA, by="county_id") %>%
  mutate(change_bin= case_when(
           mar_perc_change < -75 ~ "< -75",
           mar_perc_change >= -75 & mar_perc_change < -45 ~ "-75 - -46",
           mar_perc_change >= -45 & mar_perc_change < -30 ~ "-45 - -31",
           mar_perc_change >= -30 & mar_perc_change < 15 ~ "-30 - 16",
           mar_perc_change >= -15 & mar_perc_change <= 0 ~ "-15 - 0",
           mar_perc_change > 0  ~ "> 0"))

#quantile(sl_counties$mar_perc_change, probs= seq(0,1,0.2), na.rm = TRUE)

sl_counties$change_bin <- factor(sl_counties$change_bin, levels = c("< -75","-75 - -46","-45 - -31","-30 - 16","-15 - 0","> 0"))

vmt_change= c("#e31a1c","#fc4e2a", "#fd8d3c", "#feb24c", "#fed976", "#FCFC85")

date_labels <- c("2020-03-01"= "March 1st, 2020", "2020-04-01"= "April 1st, 2020", "2020-05-01"= "May 1st, 2020", "2020-06-01"= "June 1st, 2020", "2020-07-01"= "July 1st, 2020", "2020-10-15"= "October 15th, 2020")

sl_map <- sl_counties %>%
  filter(date== "2020-04-01"| date== "2020-10-15")

ggplot(data=sl_map)+
  geom_sf(aes(fill= change_bin))+
  coord_sf(datum = NA)+
  scale_fill_manual(values = vmt_change, na.value="#D0D0D0", drop=FALSE)+
  theme_minimal()+
  labs(fill= "% VMT Change",
       x="",
       y="",
       title= "Figure 4. Minnesota VMT by County")+
  facet_wrap(~date, ncol = 2, labeller = labeller(date= date_labels)) +
  theme(text= element_text(family= "Roboto Condensed"),
        plot.title = element_text(size=18),
        strip.text = element_text(size=16),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        legend.title = element_text(size=16),
        legend.position = c(0.96, 0.37),
        legend.text = element_text(size=16),
        legend.key.height = unit(0.6, "cm"),
        legend.key.width = unit(0.6, "cm"),
        legend.spacing.x = unit(0.1, "cm"))

The VMT data provided by Streetlight are an effective way to explore county-level travel patterns in Minnesota. There are some additional sources of mobility data that can also be useful, particularly the data provided through Apple’s Mobility Reports. This data is generated by aggregating requests through the Apple Maps app. Google makes similar data available. It is evaluated using an index, where a value of 100 the volume on January 13, 2020 (the baseline value). This data breaks down requests into three categories: driving, walking, and transit. It also has a couple weaknesses, compared to the Streetlight data.It will be biased because it tracks mapping requests rather than actual mobility. Given that people don’t typically use maps to drive to places they’re familiar with (work, neighborhood trips, etc.) it seems likely this data is sampling a biased portion of overall trips. It also has no information about trip duration, which is captured by VMT data.

Figure 5 shows how the volume of Apple mapping requests for driving have changed over the course of the pandemic in Minnesota. Volume dropped to about 50 percent of the baseline request at the beginning of April, and rebounded to about 150 percent by July. However, by September request began to decline again. It is not clear how much of this change is due to recovery from COVID-19 compared to seasonal changes in mapping requests. As a percentage, the change in mapping requests dropped less and recovered faster than the change in VMT. Importantly, the rises and falls in the mapping request data also track with the Streetlight VMT data, suggesting it is a fairly reliable measure of travel behavior.

read_csv("./VMT Analysis/Data/apple.csv") %>%
  filter(region=="Minnesota",
         transportation_type=="driving") %>%
  ggplot(aes(x=date, y= index))+
  geom_line(size=1.5, color= "#1f78b4", alpha=0.3)+
  geom_ma(ma_fun=SMA, n=7, size=1.4, linetype="solid", color="#1f78b4")+
  ylim(0, 200)+
  theme_classic()+
  labs(
    x="",
    y="Requests",
    title= "Figure 5. Change in Apple Mapping Requests (Driving)")+
  theme(text= element_text(family= "Roboto Condensed"),
        plot.title = element_text(size=18),
        panel.grid.major = element_line(colour="light grey"),
        panel.grid.major.x = element_blank(),
        axis.text.y = element_text(size=18, color = "black"),
        axis.text.x = element_text(size= 18, color = "black"),
        axis.title.y= element_text(size=20, color="black"))

The mapping request data data can also shows how mapping requests have changed for transit and walking. Figure 6 shows these changes. Transit requests declined more significantly than transit or walking requests and have not increased since June. As of mid-October, they are still at about half of the baseline. Walking requests also dropped initially, but recovered strongly over the summer. Walking requests have declined once again in the fall months. It seems likely a significant amount of the variation in walking is tied to seasonal changes.

read_csv("./VMT Analysis/Data/apple.csv") %>%
  filter(region=="Minnesota",
         transportation_type != "driving") %>%
  ggplot(aes(x=date, y= index))+
  geom_line(size=1.2, color= "#1f78b4", alpha=0.2)+
  geom_ma(ma_fun=SMA, n=7, size=1.4, linetype="solid", color="#1f78b4")+
  theme_classic()+
  labs(
    x="",
    y="Apple Maps Requests",
    title= "Figure 6. Change in Transit and Walking Requests") +
  xlab("")+
  ylab("Apple Maps Requests")+
  facet_wrap(~transportation_type, ncol=2)+
  theme(text= element_text(family= "Roboto Condensed"),
        strip.text = element_text(size=16),
        plot.title = element_text(size=18),
        panel.grid.major = element_line(colour="light grey"),
        panel.grid.major.x = element_blank(),
        axis.text.y = element_text(size=18, color = "black"),
        axis.text.x = element_text(size= 18, color = "black"),
        axis.title.y= element_text(size=20, color="black"))