library(kableExtra)
library(readxl)
library(janitor)
library(tidyverse)
library(lubridate)
library(scales)
library(viridis)
library(cowplot)Data Assignment #2 (due November 2, 2025 by 11:59pm)
This assignment addresses two related subjects: the fiscal impact of oil and gas in Alberta and the pipelines we use to transport our oil and gas.
Canada is an energy exporting country, and our energy exports are an important source of income to governments, in particular in Alberta.
Let’s start with the pipelines that transport some of those exports. The Canada Energy Regulator produces pipeline profiles for the largest oil and gas pipelines in Canada, and also tracks other export-related data. The pipeline throughput and capacity data are available here. Crude oil export data sets are available here, while gas trade data are available here.
The US EIA also tracks crude oil imports from Canada and elsewhere, and often provides much more detail than the Canadian export data. Unfortunately, their new Application Programming Interface (API) is very cumbersome and so I’ve skipped teaching you that step and downloaded some data from their API for you to use for this assignment.
There are five deliverables for this assignment, and they build on skills we have already learned:
- a graph of throughput by grade and capacity on the Enbridge mainline into Sarnia, ON;
- a graph of throughput by grade and destination on the TransMountain pipeline;
- a graph of throughput and capacity on the TC Energy Canadian Mainline;
- a graph of US oil imports by origin (Canada vs Rest of the World) and refining region (US PADDs);
- a graph of Alberta resource royalty revenue over time.
To execute my versions of the deliverables that you see below, I’ve used the following packages. Use this as a guide to set up your document:
You may also find it useful to have a look at the Functions Demo before you start this assignment.
Deliverable 1: The Enbridge Mainline
The Enbridge Mainline moves more crude oil than all other systems combined. It is also unique in that it allows Western Canadian crude to move to Ontario via the United States.
For this graph, I expect you to download data from the pipeline profiles and produce a graph in R that shows thoughput into Sarnia by grade as well as capacity at that point. You’ll need to filter the data to make it useful for your graph, and you’ll need to combine an area plot and a line plot.
You should produce something like (but not necessarily identical to) this:
enb_data<-get_pipe_data()%>%filter(key_point=="Into-Sarnia",direction_of_flow=="east") %>% mutate(
product=str_to_title(product),
product=gsub("Ngl", "NGL",product)
)
ggplot(enb_data) +
geom_area(aes(date,throughput_1000_m3_d,group = product,fill=product),position = "stack",color="black",linewidth=.25) +
geom_line(aes(date,available_capacity_1000_m3_d,color="Available Capacity"),linewidth=.85, lty="21") +
scale_fill_manual("",values=viridis(n=4,alpha=1,begin=1,end=0,option = "B",direction=1))+
scale_x_date(name=NULL,date_breaks = "2 years", date_labels = "%b\n%Y",expand=c(0,0)) +
scale_color_manual("",values=c("black"))+
scale_y_continuous(expand = c(0, 0),
sec.axis = sec_axis( trans=~.*1/.16, name="Shipments (Monthly, Thousands of Barrels per Day)")) +
guides(fill = guide_legend(order = 1) )+
theme_classic() +
theme(axis.line.x = element_line(color = "gray"),
axis.line.y = element_line(color = "gray"),
legend.key.width=unit(2,"line"),
legend.position = "bottom",
legend.box = "horizontal",
legend.text = element_text(size = 12),
plot.title = element_text(hjust=0.5,size = 14))+
labs(y="Shipments (Monthly, Thousands of Cubic Metres per Day)",x="Date",
title=paste("Enbridge Canadian Mainline Shipments by Product into Sarnia",sep=""),
caption="Source: CER Data for Enbridge Mainline, graph by Andrew Leach.")
Hints:
To add the second y axis, use
scale_y_continuous(expand = c(0, 0),sec.axis = sec_axis( trans=~.*1/.16, name="Shipments (Monthly, Thousands of Barrels per Day)"))which translates the units using the multiplier in trans=~;
You’ll probably want to download a csv and look at it. You’ll also likely want to do some fixing of cases. You can do this in R or in your CSV (but it’s easier in R);
If you’re getting a weird looking graph, you’re probably graphing by year. Check your data to see what you really should be using.
To add a capacity line, just add an extra geom_line object to your graph. This let’s me show you a neat trick too - if you name the colour in a ggplot aesthetic, this will carry through to your legends.
geom_line(aes(date,available_capacity_1000_m3_d,color="Available Capacity"),linewidth=.85, lty="21") + scale_color_manual("",values=c("black"))
Wonder what happened in 2016? Look here.
Deliverable 2: the Trans-Mountain Pipeline
The Trans-Mountain pipeline has been around since the 1950s and is the only major crude oil pipeline to cross the continental divide. Trans-Mountain is also unique in that it ships both crude oil and refined products in batches (other pipelines ship in batches too, just not with refined products in the mix). The following description is taken from the Trans-Mountain website”
Transmountain also has three delivery points: exports to the US via Sumas, deliveries to domestic refining in Burnaby, and a port delivery point at Westridge.
To see what gets delivered where on TransMountain, you should produce something like (but not necessarily identical to) this for TransMountain:
tm_data<-get_pipe_data(pipe_name="Trans-Mountain")%>% filter(key_point!="system")%>%mutate(
trade_type=str_to_title(trade_type) ,
product=str_to_title(product),
pair=paste(key_point," (",trade_type,")",sep=""))
ggplot(tm_data,aes(date,throughput_1000_m3_d,group = product,fill=product)) +
geom_area(position = "stack",,color="black",linewidth=.25) +
scale_fill_manual("",values=viridis(n=4,alpha=1,begin=.7,end=0,option = "E",direction=-1))+
scale_x_date(name=NULL,date_breaks = "2 years", date_labels = "%b\n%Y",expand=c(0,0)) +
scale_y_continuous(expand = c(0, 0),
sec.axis = sec_axis( trans=~.*1/.16, name="Shipments (Monthly, Thousands of Barrels per Day)")) +
guides(alpha = guide_legend(override.aes = list(fill=viridis(n=3,alpha=1,begin=.8,end=0,option = "E",direction=-1)[1]),order = 10) ,
fill = guide_legend(order = 1) )+
facet_wrap(~pair,ncol = 1, scales = "free_y")+
theme_classic() +
theme(panel.border = element_blank(),
panel.grid = element_blank(),
panel.grid.major.y = element_line(color = "gray",linetype="dotted"),
axis.line.x = element_line(color = "gray"),
axis.line.y = element_line(color = "gray"),
axis.text = element_text(size = 12),
axis.text.x = element_text(margin = margin(t = 10)),
axis.title = element_text(size = 12),
#axis.label.x = element_text(size=20,vjust=+5),
plot.subtitle = element_text(size = 12,hjust=0.5),
plot.caption = element_text(face="italic",size = 12,hjust=0),
legend.key.width=unit(2,"line"),
legend.position = "bottom",
legend.box = "vertical",
#legend.direction = "horizontal",
#legend.box = "horizontal",
legend.text = element_text(size = 12),
plot.title = element_text(hjust=0.5,size = 14))+
labs(y="Shipments (Monthly, Thousands of Cubic Metres per Day)",x="Date",
title=paste("Trans-Mountain Pipeline Shipments by Product and Destination",sep=""),
caption="Source: CER Data for Trans-Mountain, graph by Andrew Leach.")
Hint: use
facet_wrap(~key_point,ncol = 1, scales = "free_y")to allow different y-axes on each of the plots, and to stack the plots vertically if you like. You are not required to present them this way, but I do want to see a three-facet plot
And, take note of what’s happened since the expansion was opened and what we’re moving on the pipeline now vs. before.
Deliverable 3: The TCPL Mainline
The first pipeline controversy in Canadian political history was not about an oil pipeline, but about a gas pipeline. The Pipeline Debate of 1956 eventually led to the construction of the Canadian mainline. The Mainline also became the subject of one of the most important regulatory decisions in recent history, as the NEB refused to allow the pipeline to enter into a downward spiral and capped tolls and directed TransCanada to find new business to improve the viability of its asset. Both tolls and throughputs are available on the pipeline profile site.
To see how flows have changed on the Mainline over time at three points - the Prairies gate, the entry to the Northern Ontario Line and the start of the Eastern Triangle- you should make a graph that looks like this:
options(timeout=300)
if (!file.exists("mainline_data.rdata") ||
difftime(Sys.time(), file.info("mainline_data.rdata")$mtime, units = "days") > 10)
{
mainline_data<-get_pipe_data(pipe_name="tcpl-mainline")%>%
filter(direction_of_flow=="east")
save(mainline_data,file = "mainline_data.rdata")
}
load(file = "mainline_data.rdata")
ggplot(
#test<-
mainline_data %>%
filter(key_point %in% c("Prairies","Northern Ontario Line","Eastern Triangle - NOL Receipts"))%>%
group_by(key_point)%>%
#mutate(throughput_1000_m3_d=mean(throughput_1000_m3_d))%>%
#mutate(throughput_1000_m3_d=zoo::rollmean(throughput_1000_m3_d,30,align = "right",na.pad = TRUE))%>%
ungroup()%>%
filter(!is.na(throughput_1000_m3_d))%>%
mutate(key_point=as_factor(key_point),
key_point=fct_relevel(key_point,"Prairies","Northern Ontario Line"))%>%
rename(throughput=throughput_1000_m3_d,
capacity=capacity_1000_m3_d)) +
geom_line(aes(date,throughput/1000,color="Throughput"),linewidth=0.85) +
geom_line(aes(date,capacity/1000,color="Available Capacity"),linewidth=1.25, lty="21") +
scale_x_date(name=NULL,date_breaks = "1 year", date_labels = "%b\n%Y",expand=c(0,0)) +
scale_color_manual("",values=c("black","firebrick"))+
scale_y_continuous(expand = c(0, 0)) +
guides(alpha = guide_legend(override.aes = list(fill=viridis(n=3,alpha=1,begin=.8,end=0,option = "E",direction=-1)[1]),order = 10) ,
fill = guide_legend(order = 1) )+
theme_classic() +
expand_limits(y=0)+
facet_grid(rows=vars(key_point),scales="free_y")+
theme(panel.border = element_blank(),
panel.grid = element_blank(),
panel.grid.major.y = element_line(color = "gray",linetype="dotted"),
axis.line.x = element_line(color = "gray"),
axis.line.y = element_line(color = "gray"),
axis.text = element_text(size = 12),
#axis.label.x = element_text(size=20,vjust=+5),
plot.subtitle = element_text(size = 12,hjust=0.5),
#plot.caption = element_text(face="italic",size = 12,hjust=0),
legend.key.width=unit(2,"line"),
legend.position = "bottom",
legend.box = "horizontal",
#legend.direction = "horizontal",
#legend.box = "horizontal",
plot.caption=element_blank(),
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
plot.title = element_text(hjust=0.5,size = 14))+
labs(y="Throughput or Capacity (Millions of Cubic Metres per Day)",x="Date",
title=paste("TCPL Canadian Mainline Throughput and Capacity (Prairies and Northern Ontario)",sep=""),
caption="Source: CER Data for TCPL Canadian Mainline, graph by Andrew Leach.")
Hints:
The Mainline data files are large, so depending on the speed of your conneciton, you may get timeout errors. If you do, just add the command options(timeout=300) to your code before you call the download command. You likely also want to make sure you’re not downloading them every time you run your code, so use an if(!file.exists..) in your code as we’ve done before;
If you don’t like the look of the provided daily data, you can use a monthly average throughput, by adding something like this to your code:
group_by(month,year,key_point)%>%
mutate(throughput_1000_m3_d=mean(throughput_1000_m3_d))%>%
ungroup()You can also try something a bit more interesting and use a rolling mean. Using the
zoolibrary, you could do something like this:
group_by(key_point)%>%
#mutate(throughput_1000_m3_d=mean(throughput_1000_m3_d))%>%
mutate(throughput_1000_m3_d=zoo::rollmean(throughput_1000_m3_d,30,align = "right",na.pad = TRUE))%>%
ungroup()%>%
filter(!is.na(throughput_1000_m3_d))You should specify in your subtitle or somewhere on your graph if and how you’ve gone from the provided daily data to a smoother series. Like this:
ggplot(
#test<-
mainline_data %>%
filter(key_point %in% c("Prairies","Northern Ontario Line","Eastern Triangle - NOL Receipts"))%>%
group_by(key_point)%>%
#mutate(throughput_1000_m3_d=mean(throughput_1000_m3_d))%>%
mutate(throughput_1000_m3_d=zoo::rollmean(throughput_1000_m3_d,30,align = "right",na.pad = TRUE))%>%
ungroup()%>%
filter(!is.na(throughput_1000_m3_d))%>%
mutate(key_point=as_factor(key_point),
key_point=fct_relevel(key_point,"Prairies","Northern Ontario Line"))%>%
rename(throughput=throughput_1000_m3_d,
capacity=capacity_1000_m3_d)) +
geom_line(aes(date,throughput/1000,color="Throughput"),linewidth=0.85) +
geom_line(aes(date,capacity/1000,color="Available Capacity"),linewidth=1.25, lty="21") +
scale_fill_manual("",values=viridis(n=4,alpha=1,begin=.7,end=0,option = "E",direction=-1))+
scale_x_date(name=NULL,date_breaks = "1 year", date_labels = "%b\n%Y",expand=c(0,0)) +
scale_color_manual("",values=c("black","firebrick"))+
scale_y_continuous(expand = c(0, 0)) +
guides(alpha = guide_legend(override.aes = list(fill=viridis(n=3,alpha=1,begin=.8,end=0,option = "E",direction=-1)[1]),order = 10) ,
fill = guide_legend(order = 1) )+
theme_classic() +
expand_limits(y=0)+
facet_grid(rows=vars(key_point),scales="free_y")+
theme(panel.border = element_blank(),
panel.grid = element_blank(),
panel.grid.major.y = element_line(color = "gray",linetype="dotted"),
axis.line.x = element_line(color = "gray"),
axis.line.y = element_line(color = "gray"),
axis.text = element_text(size = 12),
#axis.label.x = element_text(size=20,vjust=+5),
plot.subtitle = element_text(size = 12,hjust=0.5),
#plot.caption = element_text(face="italic",size = 12,hjust=0),
legend.key.width=unit(2,"line"),
legend.position = "bottom",
legend.box = "horizontal",
#legend.direction = "horizontal",
#legend.box = "horizontal",
plot.caption=element_blank(),
#axis.text.x = element_blank(),
#axis.title.x = element_blank(),
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
plot.title = element_text(hjust=0.5,size = 14))+
labs(y="Throughput or Capacity (Millions of Cubic Metres per Day)",x="Date",
title=paste("TCPL Canadian Mainline Throughput and Capacity (Prairies and Northern Ontario)",sep=""),
subtitle="Trailing 30-day average throughput and monthly available capacity",
caption="Source: CER Data for TCPL Canadian Mainline, graph by Andrew Leach.")
Context: When Energy East was proposed in 2014, there was a lot less gas being moved on the Mainline. One aspect of the Energy East project was converting one pipe on the gas mainline to oil service. Do you think that option is still available today? Why or why not?
Deliverable 4: US imports of Canadian crude
We haven’t talked a lot about the US market, but it’s worth having a sense of how the US relies on imports, where that oil comes from, and where it is refined. Using data that I’ve provided for you here, you should be able to reproduce a graph that looks like this one for use imports by grade and refining region. Note that refining regions give you a better sense of where demand lies, rather than a focus on the import port or pipeline border crossing point.
assignment_data_conf<-read_csv("assignment_2_data.csv")%>%
mutate(destination_name=as_factor(destination_name))
import_plot<-
ggplot(
assignment_data_conf%>%
#mutate(origin_name=as_factor(origin_name))%>%
#mutate(origin_name=fct_relevel(origin_name,"ROW"))%>%
#mutate(origin_name=fct_recode(origin_name,"Rest of the World"="ROW")
mutate(origin_name=gsub("ROW","Rest of the World",origin_name)
))+
geom_area(aes(period,quantity,group=origin_name,fill=origin_name),color="black",linewidth=0.25,position="stack")+
scale_x_date(expand=c(0,0))+
scale_y_continuous(breaks=pretty_breaks(), expand=c(0,0))+
#look here, adding in the facet wrap and the labeler
facet_grid(cols=vars(destination_name),rows=vars(grade_name))+
expand_limits(y=c(0,140))+
#scale_fill_manual("",values=c("darkblue","red"))+
scale_fill_manual("",values=c("red","darkblue"))+
theme_classic() +
theme(panel.border = element_blank(),
panel.grid = element_blank(),
panel.grid.major.y = element_line(color = "gray",linetype="dotted"),
axis.line.x = element_line(color = "gray"),
axis.line.y = element_line(color = "gray"),
axis.text = element_text(size = 12),
axis.text.x = element_text(margin = margin(t = 10)),
axis.title = element_text(size = 12),
#axis.label.x = element_text(size=20,vjust=+5),
plot.subtitle = element_text(size = 12,hjust=0.5),
plot.caption = element_text(face="italic",size = 12,hjust=0),
legend.key.width=unit(2,"line"),
legend.position = "bottom",
legend.box = "vertical",
#legend.direction = "horizontal",
#legend.box = "horizontal",
legend.text = element_text(size = 12),
plot.title = element_text(hjust=0.5,size = 14))+
theme(text = element_text(size=16))+
theme(legend.position = "bottom", legend.box = "vertical",
legend.margin=margin(t = 0,b=0, unit='cm'))+
#t, r, b, l (To remember order, think trouble)
theme(plot.margin = unit(c(1,1,0.2,1), "cm"))+
guides(color = guide_legend(keywidth = unit(1.6,"cm"),nrow = 1),
linetype = guide_legend(keywidth = unit(1.6,"cm"),nrow = 1))+
labs(y="Imports (thousands of barrels per day)",x="",
title="US Crude Imports from Canada and the Rest of the World (ROW)",
subtitle="Imports by Refining PADD and Grade",
caption = paste("Data via US Energy Information Administration, current to ",format(max(assignment_data$period),"%B %Y"),". Graph by Andrew Leach.",sep="")
)+
NULL
import_plot_factors<-
ggplot(
assignment_data_conf%>%
mutate(origin_name=as_factor(origin_name))%>%
mutate(origin_name=fct_relevel(origin_name,"ROW"))%>%
mutate(origin_name=fct_recode(origin_name,"Rest of the World"="ROW")
#mutate(origin_name=gsub("ROW","Rest of the World",origin_name)
))+
geom_area(aes(period,quantity,group=origin_name,fill=origin_name),color="black",linewidth=0.25,position="stack")+
scale_x_date(expand=c(0,0))+
scale_y_continuous(breaks=pretty_breaks(), expand=c(0,0))+
#look here, adding in the facet wrap and the labeler
facet_grid(cols=vars(destination_name),rows=vars(grade_name))+
expand_limits(y=c(0,140))+
scale_fill_manual("",values=c("darkblue","red"))+
#scale_fill_manual("",values=c("red","darkblue"))+
theme_classic() +
theme(panel.border = element_blank(),
panel.grid = element_blank(),
panel.grid.major.y = element_line(color = "gray",linetype="dotted"),
axis.line.x = element_line(color = "gray"),
axis.line.y = element_line(color = "gray"),
axis.text = element_text(size = 12),
axis.text.x = element_text(margin = margin(t = 10)),
axis.title = element_text(size = 12),
#axis.label.x = element_text(size=20,vjust=+5),
plot.subtitle = element_text(size = 12,hjust=0.5),
plot.caption = element_text(face="italic",size = 12,hjust=0),
legend.key.width=unit(2,"line"),
legend.position = "bottom",
legend.box = "vertical",
#legend.direction = "horizontal",
#legend.box = "horizontal",
legend.text = element_text(size = 12),
plot.title = element_text(hjust=0.5,size = 14))+
theme(text = element_text(size=16))+
theme(legend.position = "bottom", legend.box = "vertical",
legend.margin=margin(t = 0,b=0, unit='cm'))+
#t, r, b, l (To remember order, think trouble)
theme(plot.margin = unit(c(1,1,0.2,1), "cm"))+
guides(color = guide_legend(keywidth = unit(1.6,"cm"),nrow = 1),
linetype = guide_legend(keywidth = unit(1.6,"cm"),nrow = 1))+
labs(y="Imports (thousands of barrels per day)",x="",
title="US Crude Imports from Canada and the Rest of the World (ROW)",
subtitle="Imports by Refining PADD and Grade",
caption = paste("Data via US Energy Information Administration, current to ",format(max(assignment_data$period),"%B %Y"),". Graph by Andrew Leach.",sep="")
)+
NULL
ggsave("../../slides/oil_transp/import_plot.png",plot = import_plot_factors,dpi=400,width=14,height=7,bg="white")
Hints:
If you want to have the PADDs in the same order as my plot, use this after you read the data:
mutate(destination_name=as_factor(destination_name)
If you want to have regions in the reverse order:
mutate(origin_name=as_factor(origin_name))%>% mutate(origin_name=fct_relevel(origin_name,"ROW"))%>% mutate(origin_name=fct_recode(origin_name,"Rest of the World"="ROW")
Think about your facet_grid elements (rows and cols), your grouping of data, and your fill aesthetic.’
Deliverable 5: Alberta Resource Royalties
Alberta’s budget is always a big economic news day, and resource royalties are always a big part of the story. In the last few fiscal years, resource royalties have hit record levels and the forecast of future royalties will drive government spending and savings decisions. To have some perspective on these decisions, let’s make a plot of Alberta resource royalties by product over time. Use the newest data here to create a graph similar to the following:
options(timeout=300)
#if(!file.exists("royalties.xlsx"))
download.file("https://open.alberta.ca/dataset/382b7a1e-9c34-47c7-9531-38e67ca5441d/resource/1373ea7a-ae5d-4ef7-acf7-7a966adb21d5/download/em-alberta-historical-royalty-revenue-2024-2025.xlsx",destfile ="royalties.xlsx",mode="wb" )
royalties<-read_excel("royalties.xlsx",range="B26:BE35")%>%clean_names()%>%slice(-1)%>%slice(-c(7,8))%>%
rename(product=1)%>%
pivot_longer(-product,names_to = "year",values_to = "royalties")%>%
mutate(year = as.numeric(gsub("x","",str_split(year, "_", simplify = TRUE)[ , 1])))
ggplot(royalties%>%mutate(product=as_factor(product)))+
geom_area(aes(year,royalties/1000,fill=product),linewidth=.25) +
scale_fill_manual("",values=colors_ua10())+
#scale_x_discrete(name=NULL,expand=c(0,0)) +
scale_color_manual("",values=c("black"))+
scale_y_continuous(expand = c(0, 0)) +
guides(alpha = guide_legend(override.aes = list(fill=viridis(n=3,alpha=1,begin=.8,end=0,option = "E",direction=-1)[1]),order = 10) ,
fill = guide_legend(order = 1) )+
theme_classic() +
theme(panel.border = element_blank(),
panel.grid = element_blank(),
panel.grid.major.y = element_line(color = "gray",linetype="dotted"),
axis.line.x = element_line(color = "gray"),
axis.line.y = element_line(color = "gray"),
axis.text = element_text(size = 12),
#axis.label.x = element_text(size=20,vjust=+5),
plot.subtitle = element_text(size = 12,hjust=0.5),
#plot.caption = element_text(face="italic",size = 12,hjust=0),
legend.key.width=unit(2,"line"),
legend.position = "bottom",
legend.box = "horizontal",
#legend.direction = "horizontal",
#legend.box = "horizontal",
plot.caption=element_blank(),
#axis.text.x = element_blank(),
#axis.title.x = element_blank(),
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
plot.title = element_text(hjust=0.5,size = 14))+
labs(y="Royalty Revenue (CA$ Billions)",x="Year",
title=paste("Alberta Royalty Revenue",sep=""),
caption="Source: Government of Alberta data, graph by Andrew Leach.")
Hints:
To turn the fiscal years into single years, I used this command:
mutate(year = as.numeric(gsub("x","",str_split(year, "_", simplify = TRUE)[ , 1]))). This splits the string at the “”, and then replaces the x with a blank. So, for example, if your year entry was ”x1970_71”, you run this command and in first grabs the first part of the outcome of splitting the strong at the ”” (x1970) and then uses gsub to replace the x with a blank.
Other than altering how the fiscal years are read in, you should not need any other new commands. You do not need to replicate the U of A colours.
Context: notice how much more important natural gas royalties used to be in Alberta!
RMD File and HTML/PDF Preparation
As before, use the basic RMD file to complete this (and future) assignments, just rename it assignment 2. You only need to submit the HTML.