quit()
reticulate::repl_python()
import subprocess
import sys
subprocess.check_call([sys.executable, "-m", "pip", "install",
"git+https://github.com/USEPA/flowsa.git@ghgi-recursive-yaml#"
"egg=flowsa"])
Y
packages = c("httr","RJSONIO","tidyr","RCurl","ggplot2","SAScii","maps",
"RColorBrewer","dplyr","matrixStats","bea.R","mapproj","devtools")
# packages used by sage but not currently on the system
new.packages = packages[!(packages %in% installed.packages()[,"Package"])]
# install the required packages not already installed
if (length(new.packages))
install.packages(new.packages, repos = "http://cran.rstudio.com/")
# required: gdxrrw
library("devtools")
install_github("GAMS-dev/gdxrrw/gdxrrw@v1.0.10")
library(devtools)
devtools::install_github("wmurphyrd/fiftystater"))
devtools::install_github("wmurphyrd/fiftystater")
wd = "~/OneDrive - Environmental Protection Agency (EPA)/epa_elasticity_databank/econ_forum_files/"
wd
rm(list = ls()[!(ls() %in% keep.items)])
list_of_packages <- c('dplyr', 'tidyr', 'stringr', 'arrow')
new_packages <- list_of_packages[!(list_of_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages, repos = "http://cran.rstudio.com/")
lapply(list_of_packages, library, character.only = TRUE)
INPUT = "C:\Users\ASchreib\OneDrive - Environmental Protection Agency (EPA)\CGETeam\data"
fl <- read_parquet(file.path(INPUT,
'emissions',
'flowsa',
'2024-08-16',
'FlowBySector',
'GHG_national_2016_m2_v2.0.3_a743c0b.parquet'))
INPUT = "C:/Users/ASchreib/OneDrive - Environmental Protection Agency (EPA)/CGETeam/data"
fls <- read_parquet(file.path(INPUT,
'emissions',
'flowsa',
'2024-08-15',
'FlowBySector',
'GHG_state_2016_m1_v2.0.3_a743c0b.parquet'))
View(fls)
fl <- read_parquet(file.path(INPUT,
'emissions',
'flowsa',
'2024-08-16',
'FlowBySector',
'GHG_national_2016_m2_v2.0.3_a743c0b.parquet'))
View(fl)
unique(fls$Flowable)
unique(fls$Year)
fls_co2 = subset(fls, Flowable %in% "Carbon dioxide") %>%
select(SectorProducedBy,Location,ActivityProducedBy,MetaSources,FlowAmount)
View(fls_co2)
sum(fls_co2$FlowAmount)
sum(fls_co2$FlowAmount)/1e12
sum(fls_co2$FlowAmount)/1e9
unique(fls_co2$MetaSources)
unique(fls_co2$SectorProducedBy)
fls_ele = subset(fls_co2) %>% filter(SectorProducedBy %in% "2211")
View(fls_ele)
fl_co2 = subset(fl, Flowable %in% "Carbon dioxide") %>%
select(SectorProducedBy,Location,ActivityProducedBy,MetaSources,FlowAmount)
View(fl_co2)
sum(fl_co2$FlowAmount)/1e9
unique(fl_co2$MetaSources)
unique(fl_co2$SectorProducedBy)
fl_ele = subset(fl_co2) %>% filter(SectorProducedBy %in% "2211")
View(fl_ele)
fl_co2 = subset(fl, Flowable %in% "Carbon dioxide") %>%
select(SectorProducedBy,Location,ActivityProducedBy,MetaSources,SourceName,FlowAmount)
fl_ele = subset(fl_co2) %>% filter(SectorProducedBy %in% "2211")
unique(fls_co2$MetaSources)
unique(fl_co2$SourceName)
write.csv(unique(fls_co2$MetaSources),file="git/test/emissions/state_metasources.csv", row.names=FALSE)
mapsource = read.csv("git/test/emissions/state_metasources.csv")
head(mapsource)
mapsource = read.csv("git/test/emissions/state_metasources.csv") %>%
rename("MetaSources"="metasources")
fls_co2 = fls_co2 %>% left_join(mapsource)
head(fls_co2)
fls_nat_co2 = fls_co2 %>% group_by(SectorProducedBy,map) %>%
summarize(pollution = sum(FlowAmount, na.rm=TRUE))
View(fls_nat_co2)
fls_map = fls_nat_co2 %>%
group_by(map) %>%
summarize(pollution = sum(pollution)) %>%
ggplot(aes(x=map, y=pollution), fill = "dodgerblue3") +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="% Change", x="Sectors") +
theme_bw() +
coord_flip() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "Inconsolata"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none")
list_of_packages <- c('dplyr', 'tidyr', 'stringr', 'arrow','ggplot')
new_packages <- list_of_packages[!(list_of_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages, repos = "http://cran.rstudio.com/")
lapply(list_of_packages, library, character.only = TRUE)
list_of_packages <- c('dplyr', 'tidyr', 'stringr', 'arrow','ggplot2')
new_packages <- list_of_packages[!(list_of_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages, repos = "http://cran.rstudio.com/")
lapply(list_of_packages, library, character.only = TRUE)
fls_map = fls_nat_co2 %>%
group_by(map) %>%
summarize(pollution = sum(pollution)) %>%
ggplot(aes(x=map, y=pollution), fill = "dodgerblue3") +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="% Change", x="Sectors") +
theme_bw() +
coord_flip() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "Inconsolata"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none")
fls_map
fls_map = fls_nat_co2 %>%
group_by(map) %>%
summarize(pollution = sum(pollution)) %>%
ggplot(aes(x=map, y=pollution), fill = "dodgerblue3") +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="% Change", x="Sectors") +
theme_bw() +
coord_flip() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none")
fls_map
fls_map
fls_map = fls_nat_co2 %>%
group_by(map) %>%
summarize(pollution = sum(pollution)) %>%
ggplot(aes(x=map, y=pollution, fill = "dodgerblue3")) +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="% Change", x="Emissions Category") +
theme_bw() +
coord_flip() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none")
fls_map
fls_map = fls_nat_co2 %>%
group_by(map) %>%
summarize(pollution = sum(pollution)) %>%
ggplot(aes(x=map, y=pollution/1e9, fill = "dodgerblue3")) +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="Million metric tonnes", x="Emissions Category") +
theme_bw() +
coord_flip() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none")
fls_map
fls_map = fls_nat_co2 %>%
group_by(map) %>%
summarize(pollution = sum(pollution)) %>%
ggplot(aes(x=map, y=pollution/1e9), fill = "dodgerblue3") +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="Million metric tonnes", x="Emissions Category") +
theme_bw() +
coord_flip() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none")
fls_map
fls_map = fls_nat_co2 %>%
group_by(map) %>%
summarize(pollution = sum(pollution)) %>%
ggplot(aes(x=map, y=pollution/1e9)) +
geom_bar(stat="identity", position=position_dodge(), color="black", fill = "dodgerblue3") +
labs(y="Million metric tonnes", x="Emissions Category") +
theme_bw() +
coord_flip() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none")
fls_map
fls_plot = fls_nat_co2 %>%
ggplot(aes(x=SectorProducedBy, y=pollution/1e9, fill=map)) +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="Million metric tonnes", x="NAICS codes", fill="Emissions/ncategory") +
theme_bw() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "right")
fls_plot
fls_plot = fls_nat_co2 %>%
ggplot(aes(x=SectorProducedBy, y=pollution/1e9, fill=map)) +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="Million metric tonnes", x="NAICS codes", fill="Emissions/ncategory") +
coord_flip() +
theme_bw() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "right")
fls_plot
fls_plot = fls_nat_co2 %>%
ggplot(aes(x=SectorProducedBy, y=pollution/1e9, fill=map)) +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="Million metric tonnes", x="NAICS codes", fill="Emissions\ncategory") +
coord_flip() +
theme_bw() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "right")
fls_plot
ggsave("git/test/emissions/summary_state_emissions.png",fls_plot,width=10,height=20)
write.csv(unique(fl_co2$MetaSources),file="git/test/emissions/national_metasources.csv", row.names=FALSE)
mapsource_nat = read.csv("git/test/emissions/national_metasources.csv")
mapsource_nat
fl_co2 = fl_co2 %>% left_join(mapsource_nat)
head(fl_co2)
fl_map = fl_co2 %>% group_by(map) %>%
summarize(pollution = sum(FlowAmount, na.rm=TRUE))
head(fl_map)
fl_map = fl_co2 %>% group_by(map) %>%
summarize(pollution = sum(FlowAmount, na.rm=TRUE)) %>%
mutate(flows = "national")
fls_map = fls_nat_co2 %>%
group_by(map) %>%
summarize(pollution = sum(pollution)) %>%
ungroup() %>%
mutate(flows = "state") %>%
rbind(., fl_map) %>%
ggplot(aes(x=map, y=pollution/1e9, fill = flows)) +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="Million metric tonnes", x="Emissions Category", fill="dataset") +
theme_bw() +
coord_flip() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "none")
fls_map
fls_map = fls_nat_co2 %>%
group_by(map) %>%
summarize(pollution = sum(pollution)) %>%
ungroup() %>%
mutate(flows = "state") %>%
rbind(., fl_map) %>%
ggplot(aes(x=map, y=pollution/1e9, fill = flows)) +
geom_bar(stat="identity", position=position_dodge(), color="black") +
labs(y="Million metric tonnes", x="Emissions Category", fill="Dataset") +
theme_bw() +
coord_flip() +
geom_hline(yintercept=0) +
theme(text = element_text(color="black",family = "mono"),
axis.text.y = element_text(color="black", hjust=0),
axis.title.y = element_text(vjust = +3),
axis.title.x = element_text(vjust = -0.75),
panel.border = element_rect(colour = "black"),
axis.line = element_line(colour = "black"),
legend.position = "right")
fls_map
ggsave("git/test/emissions/state_v_national_emissions.png",fls_map,width=8,height=6)
unique(fls$SourceName)
fls <- read_parquet(file.path(INPUT,
'emissions',
'flowsa',
'2024-08-15',
'FlowBySector',
'GHG_state_2016_m1_v2.0.3_a743c0b.parquet'))
unique(fls$SourceName)
getwd()
INPUT = file.path("OneDrive - Environmental Protection Agency (EPA)","CGETeam","data")
INPUT
INPUT = file.path(getwd(),"OneDrive - Environmental Protection Agency (EPA)","CGETeam","data")
INPUT
fs <- read_parquet(file.path(INPUT,
'emissions',
'flowsa',
'2024-08-15',
'FlowBySector',
'GHG_state_2016_m1_v2.0.3_a743c0b.parquet'))
#check if a package is installed, and if not, install it
list_of_packages <- c('dplyr', 'tidyr', 'stringr', 'arrow')
new_packages <- list_of_packages[!(list_of_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages, repos = "http://cran.rstudio.com/")
lapply(list_of_packages, library, character.only = TRUE)
INPUT = file.path(getwd(),"OneDrive - Environmental Protection Agency (EPA)","CGETeam","data")
######
fs <- read_parquet(file.path(INPUT,
'emissions',
'flowsa',
'2024-08-15',
'FlowBySector',
'GHG_state_2016_m1_v2.0.3_a743c0b.parquet'))
head(fs)
View(fs)
unique(fs$MetaSources)
fs.direct = subset(fs, MetaSources %in% "EPA_StateGHGI.direct")
head(fs.direct)
unique(fs.direct$ActivityProducedBy)
wd <- "~/git/test/loe_replication_files"
setwd(wd)
# directory to store model results
output.dir = file.path("output","largeopeneconomy")
# create output directory if not exist
if (!dir.exists(output.dir)) dir.create(output.dir)
# type of image format in which to save plots
plot.type = "png"
# dpi of saved plots
plot.dpi = 150
# add sage pallette for plot colors
sage.palette = c("#9C9F84","#A97D5D","#F7DCB4","#5C755E","#6D929B","#B7AFA3",
"#FF9966","#DEA259","#9A3334")
# ------------------------------------------------------------------------------
# load libraries
# ------------------------------------------------------------------------------
# load in the R utilities for sage
source(file.path("utilities","R_utilities.R"))
# reading gdx files package
if (!("gdx" %in% installed.packages()[,"Package"])) {
if (!("remotes" %in% installed.packages()[,"Package"])) {
install.packages("remotes")
}
library("remotes")
remotes::install_github("GAMS-dev/gdxrrw-miro")
options(repos = c(CRAN = "@CRAN@", pik = "https://rse.pik-potsdam.de/r/packages"))
install.packages("gdx")
}
# declare needed packages
list.of.packages <- c("ggplot2","dplyr","tidyr","ggrepel","stringr",
"xtable","forcats","ggforce","gdx","latex2exp")
# install packages if needed
new.packages <-
list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages))
install.packages(new.packages, repos = "http://cran.rstudio.com/")
# initialize packages
lapply(list.of.packages, library, character.only = TRUE)
# ------------------------------------------------------------------------------
# set modeling dimensions
# ------------------------------------------------------------------------------
# model file
model_file = "sage_loe.gms"
# sectors shocked
shock_sector = c("ele","bom","chm","prm","ref")
# range of shocks considered (look at $1 billion shock)
shock_size = c("1")
# large open economy trade assumptions
loe_trade = c("ref","soe","high_growth","low_growth")
# regional aggregation (default, national)
region_aggregation = "default"
# -----------------------------------------------------------------------------
# plot SOE framework for imports and exports
# -----------------------------------------------------------------------------
# SOE framework for exports -- perfectly elastic export demand curve (demand
# for US exports is perfectly elastic). Model has a supply curve for exports.
ex_sup <- function(x){x}
ex_dem_loe <- function(x){-x + 2}
png(filename=file.path(output.dir,"soe_exports_row.png"), width=600, height=600)
par(mar = c(4, 4, 0.1, 0.1))
plot(ex_sup, 0, 2, xlab="Quantity", ylab="Export Price", lwd=2, cex.lab=1.3, xaxt = "n", yaxt = "n")
plot(ex_dem_loe, 0, 2, lwd=2, lty=2, add=TRUE)
segments(0,1,2,1, lwd=2)
legend(x=1.1, y=1.8, TeX(r"(Export Supply (modeled))"), box.col="transparent", bg="white", cex=1.3)
legend(x=1.05, y=1.2, TeX(r"(Export Demand: SOE ($\epsilon \rightarrow -\infty$))"), box.col="transparent", bg="transparent", cex=1.3)
legend(x=1.1, y=0.3, TeX(r"(Export Demand: LOE ($\epsilon < 0$))"), box.col="transparent", bg="white", cex=1.3)
dev.off()
# SOE framework for imports -- perfectly elastic import supply curve (supply
# of world imports is perfectly elastic). Model has a demand curve for imports.
imp_dem <- function(x){-x + 2}
imp_sup_loe <- function(x){x}
png(filename=file.path(output.dir,"soe_imports_row.png"), width=600, height=600)
par(mar = c(4, 4, 0.1, 0.1))
plot(imp_dem, 0, 2, xlab="Quantity", ylab="Import Price", lwd=2, cex.lab=1.3, xaxt = "n", yaxt = "n")
plot(imp_sup_loe, 0, 2, lwd=2, lty=2, add=TRUE)
segments(0,1,2,1, lwd=2)
legend(x=1.1, y=0.3, TeX(r"(Import Demand (modeled))"), box.col="transparent", bg="white", cex=1.3)
legend(x=1.1, y=1.2, TeX(r"(Import Supply: SOE ($\eta \rightarrow \infty))"), box.col="transparent", bg="transparent", cex=1.3)
legend(x=1.1, y=1.8, TeX(r"(Import Supply: LOE ($\eta > 0))"), box.col="transparent", bg="white", cex=1.3)
dev.off()
setwd(file.path(wd,"gtapingams","largeopeneconomy"))
# run launching routine
system("gams run.gms")
getwd()
system("gams run.gms")
system("gams run.gms")
system("gams run.gms")
system("gams run.gms")
