##########################################################
#       Tequila and Mezcal: What about agave?
#   Author: Laura Alcocer, lalcocer@ucdavis.edu
#   Last updated: 10/07/2021
#   Purpose: Generate figures for blog article for ARE231
###########################################################


rm(list=ls())

#Set wd and load packages
setwd("ENTER YOUR WORKING DIRECTORY HERE")

# 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,ggplot2,RColorBrewer,rgdal,rgeos,raster,maptools,ggmap,httr,downloader,ggpubr)


#################################
#   Agave production in Mexico
#################################
GET("https://files.asmith.ucdavis.edu/agave_production.csv", write_disk(tf1 <- tempfile(fileext = ".csv")))
prodmx <- read.csv(tf1, header=T) 

ggplot(prodmx, aes(x=year, y=production/1000)) +
  geom_line() +
  ylab("Production (thousand tonnes)") + 
  xlab("Year") +
  ggtitle("Agave Production in Mexico") +
  theme_bw() 
ggsave("agave_production_mexico.png")

ggplot(prodmx, aes(x=year, y=harvested)) + 
  geom_line() + ylab("Area Harvested (Ha)") + 
  xlab("Year") + 
  ggtitle("Area Harvested for Agave in Mexico") + 
  theme_bw()
ggsave("agave_harvested_mexico.png")



###################################
#   Production broken down by type
###################################
GET("https://files.asmith.ucdavis.edu/agave_tequilero_production.csv", write_disk(tf2 <- tempfile(fileext = ".csv")))
tequila <- read.csv(tf2)
GET("https://files.asmith.ucdavis.edu/agave_mezcalero_production.csv", write_disk(tf3 <- tempfile(fileext = ".csv")))
mezcal <- read.csv(tf3) 


prod.type <- tequila %>% 
  rbind(mezcal, prodmx) %>% 
  base::subset(., year >= 1999 & crop != "Agave") 

ggplot(prod.type, aes(x=year, y=production/1000, fill=factor(crop,levels=c("Agave mezcalero","Agave tequilero")))) + 
  geom_area(alpha=0.8) +
  scale_fill_brewer(palette = "YlGnBu") + 
  labs(x="Year",y="Production (thousand tonnes)",fill="crops") +
  theme_bw() +
  ggtitle("Agave Production by Type")
ggsave("agave_production_bytype.png")


#################################
#   Tequila and Mezcal Exports
#################################
GET("https://files.asmith.ucdavis.edu/Indicadores20211003144430.csv", write_disk(tf4 <- tempfile(fileext = ".csv")))
exports <- read.csv(tf4) %>% 
  mutate(year = substr(Date,1,4)) %>%
  group_by(year) %>% 
  summarise(total = sum(Exports)) 

ggplot(exports, aes(x=as.integer(year), y=total)) + 
  geom_line() + 
  xlab("Year") + 
  ylab("Thousand USD") +
  ggtitle("Tequila and Mezcal Exports") + 
  theme_bw()
ggsave("tequila_exports.png")


##############################
#   Tequila production 
##############################
GET("https://files.asmith.ucdavis.edu/EstProduccionTotalTequilaAnual.csv", write_disk(tf5 <- tempfile(fileext = ".csv")))
tequilaprod <- read.csv(tf5, skip=3) %>% 
  rename(type = chart1_SeriesGroup1_label, year = chart1_SeriesGroup1_chart1_CategoryGroup1_label,
         production = chart1_SeriesGroup1_chart1_CategoryGroup1_Value_DataValue0) %>%
  mutate(production = as.numeric(production), type = as.factor(type)) %>% subset(., (type != "Total"))

ggplot(tequilaprod, aes(x=as.integer(year), y=production, fill=type)) + 
  geom_area(alpha=0.9) +
  scale_fill_brewer(palette = "YlGnBu") +  
  ylab("Production (Millions of liters)") + 
  xlab("Year") + 
  theme_bw() +
  ggtitle("Tequila Production in Volume")
ggsave("tequila_production_volume.png")

##################################
#   Differences in area harvested 
##################################
#Load in production data per tequila producing state in 2020 and 1999
GET("https://files.asmith.ucdavis.edu/production_2020.csv", write_disk(tf6 <- tempfile(fileext = ".csv")))
  prod20 <- read.csv(tf6, skip=9) %>% 
  rename(state = X.1, production = X.2, planted = Sembrada, harvested = Cosechada) %>% 
  dplyr::select(state, production, harvested) 

GET("https://files.asmith.ucdavis.edu/production_1999.csv", write_disk(tf7 <- tempfile(fileext = ".csv")))
  prod99 <- read.csv(tf7, skip=9) %>% 
  rename(state = X.1, production = X.2, planted = Sembrada, harvested = Cosechada) %>% 
  dplyr::select(state, production, harvested) 

  
#Load in shapefile for adm1 level in Mexico
download("https://files.asmith.ucdavis.edu/gadm36_MEX_1.zip",dest="./shapefile/gadm36_MEX_1.zip",mode="wb")
unzip(zipfile="shapefile/gadm36_MEX_1.zip",exdir="shapefile")  
shpfile <- readOGR(dsn = "shapefile", layer = "gadm36_MEX_1", stringsAsFactors = FALSE )
shpfile <- spTransform(shpfile, CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")) %>% gBuffer(byid=TRUE, width=0)
shp.plot <- fortify(shpfile, region="NAME_1")

map.df <- shp.plot %>% left_join(prod20, by=c("id" = "state")) 

p <- ggplot(data = map.df) + geom_polygon(aes(x=long, y=lat, group=group, fill=harvested)) +
  labs(title = "Agave Harvested in 2020") + 
  theme(plot.title = element_text(hjust = 0.5, size=14),
        plot.subtitle = element_text(hjust = 0.5, size=10),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.background = element_blank()) + 
        labs(fill = "Area Harvested (Ha)") 

map.df <- shp.plot %>% left_join(prod99, by=c("id" = "state")) 
map.df$production[is.na(map.df$production)] <- 0
map.df$harvested[map.df$harvested == 0] <- NA

q <- ggplot(data = map.df) + geom_polygon(aes(x=long, y=lat, group=group, fill=harvested)) +
  labs(title = "Agave Harvested in 1999") + 
  theme(plot.title = element_text(hjust = 0.5, size=14),
        plot.subtitle = element_text(hjust = 0.5, size=10),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.background = element_blank()) +
        labs(fill = "Area Harvested (Ha)") 

ggarrange(q,p,ncol=2,nrow=1)
ggsave("agave_harvested_map.png")