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

#############################################################
## Code to create plots in Ag Data News article            ##
##  Title:       ##
##  Link:  https://asmith.ucdavis.edu/news/    ##
#############################################################


# 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,lubridate,scales)


###################################################################################
##  FAO Stat data
###################################################################################

# I didn't have time to figure out the FAOstat API, so working with downloaded csv files

## Get and clean production and export data for Oils and Oilseeds
# Original source:  http://www.fao.org/faostat/en/#data/QCL
prod_oil <- read_csv("https://files.asmith.ucdavis.edu/FAOSTAT_data_4-21-2022_prod.csv") %>%
  select(Year,Area,Element,Item,Unit,Value) 

# Original source:  http://www.fao.org/faostat/en/#data/FBS
df_oil <- read_csv("https://files.asmith.ucdavis.edu/FAOSTAT_data_4-21-2022_balance.csv") %>%
  select(Year,Area,Element,Item,Unit,Value) 

# Prices: World Bank Pink Sheet
GET("https://thedocs.worldbank.org/en/doc/5d903e848db1d1b83e0ec8f744e55570-0350012021/related/CMO-Historical-Data-Monthly.xlsx", write_disk(tf1 <- tempfile(fileext = ".xlsx")))
df_prices <- read_xlsx(tf1,sheet="Monthly Prices",skip=6) %>%
  rename(Date=`...1`)


# Commodities/countries to plot
comms_crop <- c("Soyabeans","Sunflower seed","Rape and Mustardseed","Palm kernels")
comms_oil1 <- c("Soyabean Oil","Groundnut Oil","Sunflowerseed Oil","Rape and Mustard Oil",
                "Cottonseed Oil","Palmkernel Oil","Palm Oil","Coconut Oil","Sesameseed Oil",
                "Olive Oil","Ricebran Oil","Maize Germ Oil","Oilcrops Oil, Other")
comms_oil <- c("Sunflower Oil","Palm Oil","Rapeseed Oil","Soybean Oil","Other")
countries <- c("Ukraine","Russia", "USA","Canada","Indonesia","Malaysia","China","Brazil","Other")
countries_oil <- c("Ukraine","Russia","USA","Indonesia","China")



# Sunflower production plot
plot_sunflower_prod <- prod_oil %>%
  filter(!Area=="China") %>%
  mutate(Area=ifelse(Area == "China, mainland","China",Area)) %>%
  filter(Element %in% c("Production")) %>%
  mutate(Element=ifelse(Element=="Export Quantity","Exports",Element)) %>%
  filter(Item %in% c("Sunflower seed")) %>%
  mutate(Area=ifelse(Area == "United States of America","USA",Area)) %>%
  mutate(Area=ifelse(Area == "Russian Federation","Russia",Area)) %>%
  mutate(Area=ifelse(Area %in% c("Ukraine","Russia","USSR","USA","Argentina","China","Romania","Bulgaria","Turkey","France"),Area,"Other")) %>%
  group_by(Year,Element,Item,Area) %>%
  summarise(Value=sum(Value/1000000,na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x=Year,y=Value,fill=factor(Area,levels=c("Ukraine","Russia","USSR","USA","Argentina","China","Romania","Bulgaria","Turkey","France","Other")))) +
    geom_area(size=1)+
    labs(fill="Country", x = "Year", y = "Million Tonnes",caption="https://agdatanews.substack.com \n Source: http://www.fao.org/faostat/en/#data/QCL") +
    ggtitle(paste0("Sunflower Production by Country"))+
    theme_minimal()+
    theme(legend.position = "right", 
          plot.title = element_text(hjust = 0.5,size=16),
          text = element_text(size=14))+
    scale_fill_brewer(palette="Set3")    # Set2

# Draw and save plot
plot_sunflower_prod
ggsave(paste0("world_sunflower_prod.png"),bg="white")




# Sunflower production and export plot
plot_sunflower <- df_oil %>%
  filter(!Area=="China") %>%
  mutate(Area=ifelse(Area == "China, mainland","China",Area)) %>%
  filter(Element %in% c("Production","Export Quantity")) %>%
  mutate(Element=ifelse(Element=="Export Quantity","Exports",Element)) %>%
  filter(Item %in% c("Sunflower seed","Sunflowerseed Oil")) %>%
  mutate(Area=ifelse(Area == "United States of America","USA",Area)) %>%
  mutate(Area=ifelse(Area == "Russian Federation","Russia",Area)) %>%
  mutate(Area=ifelse(Area %in% c("Ukraine","Russia","USA","Argentina","China","Romania","Bulgaria","Turkey","France"),Area,"Other")) %>%
  group_by(Year,Element,Item,Area) %>%
  summarise(Value=sum(Value/1000,na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x=Year,y=Value,fill=factor(Area,levels=c("Ukraine","Russia","USA","Argentina","China","Romania","Bulgaria","Turkey","France","Other")))) +
    geom_area(size=1)+
    labs(fill="Country", x = "Year", y = "Million Tonnes",caption="https://agdatanews.substack.com \n Source: http://www.fao.org/faostat/en/#data/FBS") +
    facet_grid(cols=vars(factor(Element,levels=c("Production","Exports"))),rows=vars(Item),scales="free_y") +
    ggtitle(paste0("Sunflower Production and Exports"))+
    scale_x_continuous(breaks=c(2010:2019))+
    theme_minimal()+
    theme(legend.position = "right", 
          plot.title = element_text(hjust = 0.5,size=16),
          text = element_text(size=14))+
    scale_fill_brewer(palette="Set3")    # Set2

# Draw and save plot
plot_sunflower
ggsave(paste0("world_sunflower.png"),bg="white")





# World oil production plot
plot_world_oil <- df_oil %>%
  filter(!Area=="China") %>%
  mutate(Area=ifelse(Area == "China, mainland","China",Area)) %>%
  filter(Element=="Production") %>%
  filter(Item %in% comms_oil1) %>%
  mutate(Item=ifelse(Item == "Sunflowerseed Oil","Sunflower Oil",Item)) %>%  
  mutate(Item=ifelse(Item == "Rape and Mustard Oil","Rapeseed Oil",Item)) %>%  
  mutate(Item=ifelse(Item == "Soyabean Oil","Soybean Oil",Item)) %>% 
  mutate(Item=ifelse(Item %in% comms_oil,Item,"Other")) %>%
  group_by(Year,Item) %>%
  summarise(Value=sum(Value/1000,na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x=Year,y=Value,fill=factor(Item,levels=comms_oil))) +
  geom_area(size=1)+
  labs(fill="Commodity", x = "Year", y = "Million Tonnes",caption="https://agdatanews.substack.com \n Source: http://www.fao.org/faostat/en/#data/FBS") +
  ggtitle(paste0("World Vegetable Oil Production"))+
  scale_x_continuous(breaks=c(2010:2019))+
  theme_minimal()+
  theme(legend.position = "right", 
        plot.title = element_text(hjust = 0.5,size=16),
        text = element_text(size=14))+
  scale_fill_brewer(palette="Set3")    # Set2

# Draw and save plot
plot_world_oil
ggsave(paste0("oil_world_prod.png"),bg="white")




# Oil crop production plot
plot_oil_bar <- df_oil %>%
  filter(!Area=="China") %>%
  mutate(Area=ifelse(Area == "China, mainland","China",Area)) %>%
  filter(Element=="Production") %>%
  filter(Year==2019) %>%
  filter(Item %in% comms_oil1) %>%
  mutate(Area=ifelse(Area == "United States of America","USA",Area)) %>%
  mutate(Area=ifelse(Area == "Russian Federation","Russia",Area)) %>%  
  mutate(Item=ifelse(Item == "Rape and Mustardseed","Rapeseed",Item)) %>%  
  mutate(Item=ifelse(Item == "Soyabeans","Soybeans",Item)) %>% 
  mutate(Item=ifelse(Item == "Sunflowerseed Oil","Sunflower Oil",Item)) %>%  
  mutate(Item=ifelse(Item == "Rape and Mustard Oil","Rapeseed Oil",Item)) %>%  
  mutate(Item=ifelse(Item == "Soyabean Oil","Soybean Oil",Item)) %>% 
  mutate(Area=ifelse(Area %in% countries,Area,"Other")) %>%
  mutate(Item=ifelse(Item %in% comms_oil,Item,"Other")) %>%
  group_by(Year,Item,Area) %>%
  summarise(Value=sum(Value/1000,na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x=factor(Item,levels=comms_oil),y=Value,fill=factor(Area,levels=countries))) +
    geom_bar(position="stack", stat="identity")+
    labs(fill="Country", x = "Year", y = "Million Tonnes",caption="https://agdatanews.substack.com \n Source: http://www.fao.org/faostat/en/#data/FBS") +
    ggtitle(paste0("Vegetable Oil Production"))+
    theme_minimal()+
    theme(legend.position = "right", 
          plot.title = element_text(hjust = 0.5,size=16),
          text = element_text(size=14))+
    scale_fill_brewer(palette="Set3")    # Set2

# Draw and save plot
plot_oil_bar
ggsave(paste0("oil_bar_plot.png"),bg="white")






# World oil production plot
plot_prices <- df_prices %>%
  select(Date,PALM_OIL,SOYBEAN_OIL,SUNFLOWER_OIL,RAPESEED_OIL) %>%
  rename(`Palm Oil`=PALM_OIL) %>%
  rename(`Soybean Oil`=SOYBEAN_OIL) %>%
  rename(`Sunflower Oil`=SUNFLOWER_OIL) %>%
  rename(`Rapeseed Oil`=RAPESEED_OIL) %>%
  mutate(`Sunflower Oil`=as.numeric(`Sunflower Oil`)) %>%
  mutate(`Rapeseed Oil`=as.numeric(`Rapeseed Oil`)) %>%
  mutate(Date=ym(Date)) %>%
  filter(Date>=as.Date("2000-01-01")) %>%
  pivot_longer(-Date) %>%
  ggplot(aes(x=Date,y=value,color=name)) +
  geom_line(size=1)+
  labs(color="Commodity", x = "Year", y = "$US/Tonne",caption="https://agdatanews.substack.com \n Source: https://www.worldbank.org/en/research/commodity-markets") +
  ggtitle(paste0("Vegetable Oil Prices"))+
  scale_x_date(labels = date_format("%b-%Y"))+
  theme_minimal()+
  theme(legend.position = "right", 
        plot.title = element_text(hjust = 0.5,size=16),
        text = element_text(size=14))+
  scale_fill_brewer(palette="Set3")    # Set2

# Draw and save plot
plot_prices
ggsave(paste0("oil_prices.png"),bg="white")





# World oil production plot
plot_prices_2020 <- df_prices %>%
  select(Date,PALM_OIL,SOYBEAN_OIL,SUNFLOWER_OIL,RAPESEED_OIL) %>%
  rename(`Palm Oil`=PALM_OIL) %>%
  rename(`Soybean Oil`=SOYBEAN_OIL) %>%
  rename(`Sunflower Oil`=SUNFLOWER_OIL) %>%
  rename(`Rapeseed Oil`=RAPESEED_OIL) %>%
  mutate(`Sunflower Oil`=as.numeric(`Sunflower Oil`)) %>%
  mutate(`Rapeseed Oil`=as.numeric(`Rapeseed Oil`)) %>%
  mutate(Date=ym(Date)) %>%
  filter(Date>=as.Date("2020-01-01")) %>%
  pivot_longer(-Date) %>%
  ggplot(aes(x=Date,y=value,color=name)) +
    geom_line(size=1)+
    labs(color="Commodity", x = "Year", y = "$US/Tonne",caption="https://agdatanews.substack.com \n Source: https://www.worldbank.org/en/research/commodity-markets") +
    ggtitle(paste0("Vegetable Oil Prices"))+
      scale_x_date(labels = date_format("%b-%Y"))+
      theme_minimal()+
      theme(legend.position = "right", 
            plot.title = element_text(hjust = 0.5,size=16),
            text = element_text(size=14))+
      scale_fill_brewer(palette="Set3")    # Set2

# Draw and save plot
plot_prices_2020
ggsave(paste0("oil_prices_2020.png"),bg="white")