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

#############################################################
## Code to create plots in Ag Data News article            ##
##  Title:A California Avocado Hass Taken Over the W             ##
##  Link:  https://asmith.ucdavis.edu/news/california-avocado-hass-taken-over-world   ##
#############################################################

# 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(tidyverse,readxl,RColorBrewer,httr,blscrapeR,rnassqs,data.table)


########################################################
##  Avocados
########################################################

# Load world production data
df_world <- read_csv("https://files.asmith.ucdavis.edu/FAOSTAT_data_5-3-2021_avo.csv")
countries <- c("Mexico","Dominican Republic","Peru","Colombia","Indonesia","Kenya","Brazil","Haiti","United States of America")

#  Make stacked area plot
world_plot <- df_world %>%
  drop_na(Value) %>%
  filter(Element=="Production") %>%
  select(Area,Year,Value) %>%
  mutate(Area = ifelse(Area %in% countries,Area,"Other")) %>%
  mutate(Area = ifelse(Area=="United States of America","USA",Area))%>%
  group_by(Year,Area) %>%
  summarize(Value=sum(Value)) %>%
  ungroup()%>%
  ggplot(aes(x=Year,y=Value/1000000,fill=factor(Area,levels=c("Mexico","Dominican Republic","Peru","Colombia","Indonesia","Kenya","Brazil","Haiti","Other","USA")))) +
    geom_area()+
    ggtitle(paste0("World Avocado Production (8 largest producers + USA)"))+
    theme_minimal()+
    scale_fill_brewer(palette = "Set3") +
    labs(x="Year",y="Million Metric Tons",fill="Country")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# view and save plot
world_plot
ggsave("world_avo.png")  



#########################################################
#  US Production
#########################################################
# api for NASS (get yours at quickstats.nass.usda.gov/api)
#apikey <- "INSERT YOUR API KEY HERE"
apikey <- "344EEC10-9B53-3725-9F7A-F5F99CC06797"

nassqs_auth(apikey) # authentificate your api key
nassqs_params() # check what variables are there
comm_list <- nassqs_param_values("commodity_desc") # get commodity list


# set years
# you will receive an "bad request - invalid query" error if there are no data for your chosen years.
years <- c(1960:2019)

# set parameters for API calls except for years 
# (we are not allowed to request many records at once, so I iterate API calls by years)
params <-list(source_desc="SURVEY",commodity_desc="AVOCADOS",agg_level_desc=c("NATIONAL","STATE")) 

# create data based on the parameters and years using 'nassqs'
df <- map(years, function(x){params[["year"]] <- x; nassqs(params)}) %>% bind_rows() %>% as.data.table()


# bar plot of production by state
state_plot <- df %>%
  filter(short_desc=="AVOCADOS - PRODUCTION, MEASURED IN TONS") %>%
  filter(!agg_level_desc=="NATIONAL") %>%
  select(state_alpha,year,Value)%>%
  mutate(Value = as.numeric(gsub(",", "", Value))) %>%
  ggplot(aes(x=year,y=(2000/2204)*Value/1000000,fill=factor(state_alpha,levels=c("HI","FL","CA")))) +
    geom_bar(stat="identity",position="stack")+
    ggtitle(paste0("US Avocado Production by State"))+
    theme_minimal()+
    scale_fill_brewer(palette = "Set2") +
    labs(x="Year",y="Million Metric Tons",fill="Country")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

state_plot
ggsave("avo_state.png")  


#############################################
#  supply and use
############################################
GET("https://www.ers.usda.gov/webdocs/DataFiles/54499/FruitYearbookSupplyandUtilization_GTables.xlsx?v=6304.5", write_disk(tf1 <- tempfile(fileext = ".xlsx")))
df_use <- read_xlsx(tf1,sheet="tab-g3") %>%
  drop_na(`...2`) %>%
  rename(year=`Table G-3--Fresh avocados: Supply and utilization, 1980/81 to date`)


# bar plot of supply
supply_plot <- df_use %>%
  rename(Domestic=`...2`) %>%
  rename(Imports=`...3`) %>%
  select(year,Domestic,Imports) %>%
  drop_na(year) %>%
  filter(!year=="Season 1") %>%
  pivot_longer(-year) %>%
  mutate(value=as.numeric(value)/2204) %>%
  mutate(year=as.numeric(str_sub(year,1,4)))%>%
  ggplot(aes(x=year,y=value,fill=factor(name,levels=c("Imports","Domestic")))) +
    geom_area()+
    ggtitle(paste0("Source of Avocados Consumed in US"))+
    theme_minimal()+
    scale_fill_brewer(palette = "Set2") +
    labs(x="Year",y="Million Metric Tons",fill="Source")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

  
supply_plot
ggsave("avo_source.png")



#############################################
#  price
############################################

# get annual CPI from BLS
df_cpi <- setDT(inflation_adjust(2019))[,.(year=as.integer(year),avg_cpi)] # convert to data.table object and convert year as integer.


#  Make price plot
price_plot <- df %>%
  filter(short_desc=="AVOCADOS - PRICE RECEIVED, MEASURED IN $ / TON") %>%
  filter(agg_level_desc=="NATIONAL") %>%
  select(state_alpha,year,Value)%>%
  mutate(Value = as.numeric(gsub(",", "", Value))) %>%
  left_join(df_cpi) %>%
  mutate(real_price=Value/(avg_cpi/tail(avg_cpi,1)))%>%
  filter(year>=1980)%>%
  ggplot(aes(x=year,y=real_price)) +
    geom_line(color="blue",size=1)+
    ggtitle(paste0("US Average Price Received for Avocados ($2019)"))+
    theme_minimal()+
    scale_fill_brewer(palette = "Set2") +
    labs(x="Year",y="$2019/Ton")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

price_plot
ggsave("avo_price.png")  



#  Make acreage plot
acre_plot <- df %>%
  filter(short_desc=="AVOCADOS - ACRES BEARING") %>%
  filter(agg_level_desc=="NATIONAL") %>%
  select(state_alpha,year,Value)%>%
  mutate(Value = as.numeric(gsub(",", "", Value))/1000) %>%
  filter(year>=1980)%>%
  ggplot(aes(x=year,y=Value)) +
    geom_line(color="blue",size=1)+
    ggtitle(paste0("US Avocado Bearing Acreage"))+
    theme_minimal()+
    scale_fill_brewer(palette = "Set2") +
    labs(x="Year",y="000 Acres")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

acre_plot
ggsave("avo_acres.png")