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

#############################################################
## Code to create plots in Ag Data News article            ##
##  Title: One Bee for Every 20 Nuts                         ##
##  Link:  https://asmith.ucdavis.edu/news/bees-per-almond       ##
#############################################################


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

#Load agstar data from EPA
GET("https://files.asmith.ucdavis.edu/AlmondShipments.xlsx", write_disk(tf1 <- tempfile(fileext = ".xlsx")))
df <- read_excel(tf1) %>%
  select(-`20/21 vs 19/20`) 

# Subset domestic shipments
domestic <- slice_head(df,n=12) %>%
  pivot_longer(-Domestic) %>%
  mutate(year=as.numeric(str_extract(name, "^.{2}"))+2000) %>%
  mutate(year=ifelse(Domestic %in% c("January","February","March","April","May","June","July"),year+1,year)) %>%
  select(-name) %>%
  mutate(date=mdy(paste(Domestic,"/","01","/",year))) %>%
  arrange(date) %>%
  select(-year,-Domestic) %>%
  mutate(dest="Domestic") %>%
  mutate(value=as.numeric(value))

# Subset exports
export <- slice(df,15:26) %>%
  pivot_longer(-Domestic) %>%
  mutate(year=as.numeric(str_extract(name, "^.{2}"))+2000) %>%
  mutate(year=ifelse(Domestic %in% c("January","February","March","April","May","June","July"),year+1,year)) %>%
  select(-name) %>%
  mutate(date=mdy(paste(Domestic,"/","01","/",year))) %>%
  arrange(date) %>%
  select(-year,-Domestic) %>%
  mutate(dest="Export")%>%
  mutate(value=as.numeric(value))


# Make graph
plot_bar <- rbind(domestic,export) %>%
  ggplot(aes(x=date,y=value,fill=dest)) +
    geom_bar(position="stack",stat="identity") +
    ggtitle(paste0("California Monthly Almond Shipments"))+
    theme_minimal()+
    scale_fill_brewer(palette="Set2")+
    labs(x="Year",y="Million lbs",fill="Destination")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14),axis.text.x = element_text(angle = 0, hjust=1))



plot_bar
ggsave("almond_shipments.png")




########
# api for NASS (get yours at quickstats.nass.usda.gov/api)
apikey <- "INSERT YOUR API KEY HERE"

nassqs_auth(apikey) # authentificate your api key
nassqs_params() # check what variables are there
comm_list <- nassqs_param_values("commodity_desc") # get commodity list
grep("HONEY",comm_list,value=T) # show all commodity names that contain "CARROTS" (we have only one here).



# set years
# you will receive an "bad request - invalid query" error if there are no data for your chosen years.
years <- 1987:2020

# 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)
varh <-"HONEY, BEE COLONIES - INVENTORY, MEASURED IN COLONIES"
params <-list(short_desc = varh, source_desc="SURVEY", agg_level_desc="NATIONAL") 

# create data based on the parameters and years using 'nassqs'
dfb <- map(years, function(x){params[["year"]] <- x; nassqs(params)}) %>% bind_rows() %>% as.data.table() %>%
  mutate(Bees = as.numeric(gsub(",", "", Value))) %>%
  filter(freq_desc=="ANNUAL") %>%
  select(year, Bees) 
  
# 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)
vara <- "ALMONDS, UTILIZED, SHELLED - PRODUCTION, MEASURED IN LB"
params <-list(short_desc = vara, source_desc="SURVEY", agg_level_desc="NATIONAL") 

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

# create data based on the parameters and years using 'nassqs'
dfa <- map(years, function(x){params[["year"]] <- x; nassqs(params)}) %>% bind_rows() %>% as.data.table()  %>%
  mutate(Almonds = as.numeric(gsub(",", "", Value))) %>%
  filter(reference_period_desc=="YEAR") %>%
  select(year, Almonds)



# Make graph
plot_line <- dfb %>%
  left_join(dfa) %>%
  mutate(Bees=Bees/1000000) %>%
  mutate(Almonds=2.5*Almonds/1000000000) %>%
  ggplot(aes(x=year,color=name)) +
    geom_line(aes(y=Bees),size=1,color="blue") +
    geom_line(aes(y=Almonds),size=1,color="red") +
    ggtitle(paste0("Honey Bee Numbers and Almond Production"))+
    theme_minimal()+
    scale_y_continuous(name = "Million Honey Bee Colonies",breaks=seq(0,7,1),limits=c(0,7), sec.axis = sec_axis(trans=~.*.4, name="Billion Pounds of Almonds",breaks=seq(0,3,0.4) )) + 
    labs(x="Year",y="Million lbs",fill="Destination")+
    theme(axis.title.y = element_text(color = "blue", size=14),
        axis.text.y=element_text(colour="blue",size=12),
        axis.title.y.right = element_text(color = "red", size=14),
        axis.text.y.right=element_text(colour="red",size=12),
        axis.text.x=element_text(colour="black",size=12),
        plot.title = element_text(hjust = 0.5,size=18))



plot_line
ggsave("almond_bees.png")