VMT Change During the COVID-19 Pandemic
Nate Bean
10/22/2020
library(tidyverse)
options(dplyr.summarise.inform = FALSE)
library(tidycensus)
library(leaflet)
library(tidyquant)
library(extrafont)
library(fiftystater)
library(readxl)
library(janitor)
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))#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), date_labels = paste0("%b"," ","1st"), limits = as.Date(c(NA, "2020-10-16")),
date_breaks = "1 month")+
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= "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= 15, color = "black", angle= 20, vjust= 0.5),
axis.title.y= element_text(size=20, color="black", margin = margin(r=10)))
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)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")), date_labels = paste0("%b"," ","1st"),
date_breaks = "1 month")+
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 = "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= 15, color = "black", angle= 20, vjust= 0.5),
axis.title.y= element_text(size=18, color="black"),
legend.position = "bottom",
legend.margin = margin(t= -15),
legend.text = element_text(size=18))
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= "Minnesota's 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"))
district_pal <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00')
districts <- read_xlsx("./VMT Analysis/Data/Raw Data/V1.xlsx") %>%
filter(Year==2000,
EntityName != "Statewide/Unattributable") %>%
select(EntityName, ATP) %>%
rename(county= EntityName,
District= ATP)
sl_atp <- streetlight %>%
inner_join(districts, by="county") %>%
group_by(date, District) %>%
summarise(vmt_total= sum(vmt))
ggplot(sl_atp, aes(x=date, y= vmt_total/1000000, color=factor(District)))+
geom_ma(ma_fun=SMA, n=7, size=1, linetype="solid")+
scale_color_manual(values = rev(c(district_pal)))+
scale_x_date(limits = as.Date(c(NA, "2020-10-16")), expand= c(0,1),
date_labels = paste0("%b"," ","1st"), date_breaks = "1 month")+
ylim(0, 150)+
theme_classic()+
labs(
x="",
y="VMT",
color= "",
title = "VMT by Transportation District",
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= 15, color = "black", vjust=0.5, angle=20),
axis.title.y= element_text(size=18, color="black"),
legend.position = "right",
#legend.margin = margin(t= -15),
legend.text = element_text(size=18))
Next, this figure shows the change in per capita VMT by district. There is a pretty clear break between the metro area and the remaining transportation districts around the beginning of July.
#https://mn.gov/admin/demography/data-by-topic/population-data/our-estimates/
library(janitor)
pop <- read_xlsx("./VMT Analysis/Data/Control Data/mn_county_estimates.xlsx") %>%
clean_names() %>%
rename(county= county_name,
population= total_population_2019) %>%
select(county, population)
sl_atp_pop <- streetlight %>%
inner_join(districts, by="county") %>%
inner_join(pop, by="county") %>%
group_by(date, District) %>%
summarise(vmt_total= sum(vmt),
pop_total= sum(population),
vmt_pop_total= vmt_total/pop_total)
ggplot(sl_atp_pop, aes(x=date, y= vmt_pop_total, color=factor(District)))+
geom_ma(ma_fun=SMA, n=7, size=1, linetype="solid")+
scale_color_manual(values = rev(c(district_pal)))+
scale_x_date(limits = as.Date(c(NA, "2020-10-16")), expand= c(0,1), date_labels = paste0("%b"," ","1st"),
date_breaks = "1 month")+
theme_classic()+
labs(
x="",
y="VMT",
color= "",
title = "VMT per capita by Transportation District",
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=15, color = "black", vjust=0.5, angle=20),
axis.title.y= element_text(size=18, color="black"),
legend.position = "right",
#legend.margin = margin(t= -15),
legend.text = element_text(size=18))
Next, we can also reproduce the old chart that compared the metro to Greater MN, but in per capita terms.
sl_metro_pop <- streetlight %>%
inner_join(pop, by="county") %>%
group_by(date, region) %>%
summarise(vmt_total= sum(vmt),
pop_total= sum(population),
vmt_pop_total= vmt_total/pop_total)
ggplot(sl_metro_pop, aes(x=date, y= vmt_pop_total, color=region))+
geom_ma(ma_fun=SMA, n=7, size=2, linetype="solid")+
scale_color_manual(values = c(change_pal))+
scale_x_date(date_labels = paste0("%b."," ","1st"), limits = as.Date(c(NA, "2020-10-16")),
date_breaks = "1 month")+
theme_classic()+
labs(
x="",
y="VMT",
color= "",
title = "Regional VMT per capita",
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= 15, color = "black", vjust=0.5, angle=20),
axis.title.y= element_text(size=18, color="black"),
legend.position = "bottom",
legend.margin = margin(t= -15),
legend.text = element_text(size=18))
Below is a map of race by county, with data coming from the U.S. Census Bureau population estimates for 2019. I defined non-white as those who are non-Hispanic, white alone. This figure can obviously be improved, but I think it’s best to use a continuous color scale, rather than breaks the counties into bins.
# i defined white as "not hispanic, white-alone"
race <- read_csv("./VMT Analysis/Data/Control Data/race_by_county.csv") %>%
clean_names() %>%
filter(
year==12,
agegrp==0) %>%
mutate(nw_pop = (nhwa_male+nhwa_female)/tot_pop) %>% #,
#county= tolower(str_remove(ctyname, " County"))) %>%
select(ctyname, nw_pop)
race_map <- counties %>%
mutate(ctyname= str_remove(NAME, ", Minnesota")) %>%
inner_join(race, by= "ctyname") %>%
mutate(nw_pop= 1- nw_pop)
ggplot(data=race_map)+
geom_sf(aes(fill= nw_pop))+
coord_sf(datum = NA)+
#scale_fill_manual(values = vmt_change, na.value="#D0D0D0", drop=FALSE)+
theme_minimal()+
labs(fill= "% Non-white",
x="",
y="",
title= "Percent non-white by County")+
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"))