rm(list=ls())
setwd("E:/Dropbox/Research/_DeLoach/Ag_Data_News/")

#############################################################
## Code to create plots in Ag Data News article            ##
##  Title:  New Zealand Agricultural Proposes To Tax Itself, Or Does It?           ##
##  Link:  https://asmith.ucdavis.edu/news/new-zealand-agricultural-sector-proposes-tax-itself    ##
#############################################################


# If a package is installed, it will be loaded. If any are not, the missing package(s) will be installed 
# from CRAN and then loaded.
# pacman::p_load(purrr, tidyverse, RColorBrewer, data.table, rnassqs, blscrapeR, readxl, httr)
pacman::p_load(tidyverse,readxl,RColorBrewer,httr)


## Get emissions data
# Original source:  https://www.climatewatchdata.org/data-explorer/historical-emissions
df_emissions <- read_csv("https://files.asmith.ucdavis.edu/CW_HistoricalEmissions_CAIT.csv") %>%
  pivot_longer(cols="1990":"2019") %>%
  mutate(Year=as.numeric(name))


## Bar plot of emissions by gas

plot_country <- c("NZL")
plot_sector <- c("Agriculture","Energy","Land-Use Change and Forestry","Other")

# Make plot
plot_bar_gas <- df_emissions %>%
  filter(Country %in% plot_country)%>%
  filter(Gas %in% c("CO2","CH4","N2O")) %>%
  filter(Sector %in% c("Agriculture","Bunker Fuels","Energy","Industrial Processes","Land-Use Change and Forestry","Waste")) %>%
  mutate(Sector=ifelse(Sector %in% plot_sector,Sector,"Other")) %>%
  arrange(Year,Gas,Sector)%>%
  group_by(Country,Year,Sector,Gas) %>%
  summarise(Emissions=sum(value)) %>%
  ungroup() %>%
  mutate(Gas=ifelse(Gas=="CO2","Carbon Dioxide",Gas)) %>%
  mutate(Gas=ifelse(Gas=="N2O","Nitrous Oxide",Gas)) %>%
  mutate(Gas=ifelse(Gas=="CH4","Methane",Gas)) %>%
  ggplot(aes(x=Year,y=Emissions,fill=Sector)) +
    #geom_area()+
    geom_bar(position = "stack", stat="identity")+
    facet_grid(rows=vars(Gas)) +
    labs(fill="Sector", x = "Year", y = "Million Metric Tonnes CO2e",caption="Data Source: https://www.climatewatchdata.org/data-explorer/historical-emissions")+
    ggtitle(paste0("GHG Emissions: ", plot_country))+
    theme_minimal()+
    theme(legend.position = "right", 
          plot.title = element_text(hjust = 0.5))+
    scale_fill_brewer(palette="Dark2")    # Set2

# Draw and save plot
plot_bar_gas
ggsave(paste0("GHG_emissions_gas_",plot_country,".png"),bg="white")


## Bar plot of Country emissions by sector
plot_country_sector <- c("USA","NZL")

# Make plot
plot_bar_sector <- df_emissions %>%
  filter(Country %in% plot_country_sector)%>%
  filter(Gas %in% c("CO2","CH4","N2O")) %>%
  filter(Sector %in% c("Agriculture","Bunker Fuels","Energy","Industrial Processes","Land-Use Change and Forestry","Waste")) %>%
  mutate(Sector=ifelse(Sector %in% plot_sector,Sector,"Other")) %>%
  arrange(Year,Gas,Sector)%>%
  group_by(Country,Year,Sector,Gas) %>%
  summarise(Emissions=sum(value)) %>%
  ungroup() %>%
  mutate(Gas=ifelse(Gas=="CO2","Carbon Dioxide",Gas)) %>%
  mutate(Gas=ifelse(Gas=="N2O","Nitrous Oxide",Gas)) %>%
  mutate(Gas=ifelse(Gas=="CH4","Methane",Gas)) %>%
  ggplot(aes(x=Year,y=Emissions,fill=Sector)) +
    #geom_area()+
    geom_bar(position = "stack", stat="identity")+
    facet_grid(rows=vars(Country),scales="free_y") +
    labs(fill="Sector", x = "Year", y = "Million Metric Tonnes CO2e",caption="Data Source: https://www.climatewatchdata.org/data-explorer/historical-emissions")+
    ggtitle("GHG Emissions for New Zealand and USA by Sector")+
    theme_minimal()+
    theme(legend.position = "right", 
          plot.title = element_text(hjust = 0.5))+
    scale_fill_brewer(palette="Dark2")    # Set2

# Draw and save plot
plot_bar_sector
ggsave("GHG_emissions_sector.png",bg="white")


## Bar plot of Country emissions by gas
plot_country_gas <- c("USA","NZL")

# Make plot
plot_bar_gas2 <- df_emissions %>%
  filter(Country %in% plot_country_gas)%>%
  filter(Gas %in% c("CO2","CH4","N2O")) %>%
  filter(Sector %in% c("Agriculture","Bunker Fuels","Energy","Industrial Processes","Land-Use Change and Forestry","Waste")) %>%
  mutate(Sector=ifelse(Sector %in% plot_sector,Sector,"Other")) %>%
  arrange(Year,Gas,Sector)%>%
  group_by(Country,Year,Gas) %>%
  summarise(Emissions=sum(value)) %>%
  ungroup() %>%
  mutate(Gas=ifelse(Gas=="CO2","Carbon Dioxide",Gas)) %>%
  mutate(Gas=ifelse(Gas=="N2O","Nitrous Oxide",Gas)) %>%
  mutate(Gas=ifelse(Gas=="CH4","Methane",Gas)) %>%
  ggplot(aes(x=Year,y=Emissions,fill=Gas)) +
    #geom_area()+
    geom_bar(position = "stack", stat="identity")+
    facet_grid(rows=vars(Country),scales="free_y") +
    labs(fill="Gas", x = "Year", y = "Million Metric Tonnes CO2e",caption="Data Source: https://www.climatewatchdata.org/data-explorer/historical-emissions")+
    ggtitle("GHG Emissions for New Zealand and USA by Gas")+
    theme_minimal()+
    theme(legend.position = "right", 
          plot.title = element_text(hjust = 0.5))+
    scale_fill_brewer(palette="Dark2")    # Set2

# Draw and save plot
plot_bar_gas2
ggsave("GHG_emissions_gas2.png",bg="white")