rm(list=ls())
setwd("E:/Dropbox/Research/_DeLoach/Ag_Data_News/")
#############################################################
## Code to create plots in Ag Data News article            ##
##  Title: A Crushing Crush Report           ##
##  Link:  https://asmith.ucdavis.edu/news/crushing-grape-report    ##
#############################################################


# 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,ggpubr,rnassqs,data.table,tigris,tmap,pals)

# The crush report tables are in separate Excel files inside a zipped file at https://www.nass.usda.gov/Statistics_by_State/California/Publications/Specialty_and_Other_Releases/Grapes/Crush/Reports/index.php
# I downloaded and unzipped the 2019 and 2020 data, and posted them to my ftp site to be read by this code.


# Table 2 from 2019 crush report
GET("https://files.asmith.ucdavis.edu/2019erratagcbtb02.xls", write_disk(tf1 <- tempfile(fileext = ".xls")))
region_prod_2019 <- read_excel(tf1,range = cell_rows(3:159)) %>%
  filter(!is.na(`Type and Variety...1`)) %>%
  filter(!is.na(`1`)) %>%
  select(-`Type and Variety...12`) %>%
  mutate(`1`=as.numeric(`1`)) %>%
  mutate(`11`=as.numeric(`11`)) %>%
  rename(`State Total`=`2019 State Total`) %>%
  rename(`Prev State Total`=`2018 State Total`) %>%
  mutate(year=2019)


# Table 2 from 2020 crush report
GET("https://files.asmith.ucdavis.edu/Prelim_2020_gcbtb02.xls", write_disk(tf2 <- tempfile(fileext = ".xls")))
region_prod_2020 <- read_excel(tf2,range = cell_rows(3:156)) %>%
  filter(!is.na(`Type and Variety...1`)) %>%
  filter(!is.na(`1`)) %>%
  select(-`Type and Variety...12`) %>%
  mutate(`1`=as.numeric(`1`)) %>%
  mutate(`11`=as.numeric(`11`)) %>%
  rename(`State Total`=`2020 State Total`) %>%
  rename(`Prev State Total`=`2019 State Total`) %>%
  mutate(year=2020)



# Table 4 from 2019 crush report
GET("https://files.asmith.ucdavis.edu/2019erratagcbtb04.xls", write_disk(tf1a <- tempfile(fileext = ".xls")))
region_purch_2019 <- read_excel(tf1a,range = cell_rows(4:154)) %>%
  filter(!is.na(`Type and Variety...1`)) %>%
  filter(!is.na(`1`)) %>%
  select(-`Type and Variety...12`) %>%
  mutate(`1`=as.numeric(`1`)) %>%
  mutate(`11`=as.numeric(`11`)) %>%
  rename(`State Total`=`2019 State Total`) %>%
  rename(`Prev State Total`=`2018 State Total`) %>%
  mutate(year=2019)


# Table 4 from 2020 crush report
GET("https://files.asmith.ucdavis.edu/Prelim_2020_gcbtb04.xls", write_disk(tf2a <- tempfile(fileext = ".xls")))
region_purch_2020 <- read_excel(tf2a,range = cell_rows(4:151)) %>%
  filter(!is.na(`Type and Variety...1`)) %>%
  filter(!is.na(`1`)) %>%
  select(-`Type and Variety...12`) %>%
  mutate(`1`=as.numeric(`1`)) %>%
  mutate(`11`=as.numeric(`11`)) %>%
  rename(`State Total`=`2020 State Total`) %>%
  rename(`Prev State Total`=`2019 State Total`) %>%
  mutate(year=2020)


# Table 6 from 2019 crush report
GET("https://files.asmith.ucdavis.edu/2019erratagcbtb06.xls", write_disk(tf3 <- tempfile(fileext = ".xls")))
region_val_2019 <- read_excel(tf3,range = cell_rows(4:154)) %>%
  filter(!is.na(`Type and Variety...1`)) %>%
  filter(!is.na(`1`)) %>%
  select(-`Type and Variety...12`) %>%
  mutate(`1`=as.numeric(`1`)) %>%
  mutate(`11`=as.numeric(`11`)) %>%
  rename(`State Total`=`2019 State Total`) %>%
  rename(`Prev State Total`=`2018 State Total`) %>%
  mutate(year=2019) 

# Table 6 from 2020 crush report
GET("https://files.asmith.ucdavis.edu/Prelim_2020_gcbtb06.xls", write_disk(tf4 <- tempfile(fileext = ".xls")))
region_val_2020 <- read_excel(tf4,range = cell_rows(4:151)) %>%
  filter(!is.na(`Type and Variety...1`)) %>%
  filter(!is.na(`1`)) %>%
  select(-`Type and Variety...12`) %>%
  mutate(`1`=as.numeric(`1`)) %>%
  mutate(`11`=as.numeric(`11`)) %>%
  rename(`State Total`=`2020 State Total`) %>%
  rename(`Prev State Total`=`2019 State Total`) %>%
  mutate(year=2020)


# Time series from crush reports (pasted into Excel from unnumbered table on page 2 of 2020 report)
GET("https://files.asmith.ucdavis.edu/crush_2020.xlsx", write_disk(tf5 <- tempfile(fileext = ".xlsx")))
crush_2020 <- read_excel(tf5)



###################################
### Make crush bar chart
###################################
plot_bar <- crush_2020 %>%
  pivot_longer(-year) %>%
  filter(str_detect(name, "tons")) %>%
  mutate(name=gsub(" type -tons", "", name)) %>%
  mutate(name=gsub(" types -tons", "", name))%>%
  filter(!name=="All"&!name=="Total wine") %>% 
  mutate(name=factor(name,levels=c("Red wine","White wine","Table","Raisin"))) %>%
  ggplot(aes(x=year,y=value,fill=name)) +
    geom_bar(position="stack", stat="identity") +
    theme_minimal()+
    ggtitle(paste0("California Grape Crush"))+
    labs(x="Year",y="Thousand Tons",fill="Variety")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw plot
plot_bar
ggsave("crush_bar.png")


###################################
### Make time series price plot
###################################
plot_ts <- crush_2020 %>%
  pivot_longer(-year) %>%
  filter(str_detect(name, "price")) %>%
  mutate(name=gsub(" type -price", "", name)) %>%
  mutate(name=gsub(" types -price", "", name)) %>%
  filter(!name=="Total wine") %>% 
  mutate(name=factor(name,levels=c("Red wine","White wine","Table","Raisin","All"))) %>%
  ggplot(aes(x=year,y=value,color=name)) +
    geom_line(size=1) +
    theme_minimal()+
    ggtitle(paste0("California Grape Crush Prices"))+
    labs(x="Year",y="$/Tons",color="Variety")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw plot
plot_ts
ggsave("crush_timeseries.png")


##################################
### Maps
##################################
# get county boundary and merge with growing regions
cb <- counties(cb=T) %>% # for state map, use `states(cb=T)`
  filter(STATEFP %in% c("06")) %>%
  mutate(fips = as.integer(GEOID)) %>%
  mutate(region=ifelse(NAME %in% c("Mendocino"),1,NA)) %>%
  mutate(region=ifelse(NAME %in% c("Lake"),2,region)) %>%
  mutate(region=ifelse(NAME %in% c("Marin","Sonoma"),3,region)) %>%
  mutate(region=ifelse(NAME %in% c("Napa"),4,region)) %>%
  mutate(region=ifelse(NAME %in% c("Solano"),5,region)) %>%
  mutate(region=ifelse(NAME %in% c("San Francisco","Alameda","Santa Clara","San Mateo","Santa Cruz", "Contra Costa"),6,region)) %>%
  mutate(region=ifelse(NAME %in% c("Monterey","San Benito"),7,region)) %>%
  mutate(region=ifelse(NAME %in% c("San Luis Obispo","Ventura","Santa Barbara"),8,region)) %>%
  mutate(region=ifelse(NAME %in% c("Siskiyou","Shasta","Lassen","Humboldt","Butte","Modoc","Glenn","Plumas","Sierra","Colusa","Yuba","Trinity","Del Norte","Tehama","Sutter"),9,region)) %>%
  mutate(region=ifelse(NAME %in% c("Nevada","Tuolumne","Placer","Calaveras","Amador","El Dorado","Mariposa"),10,region)) %>%
  mutate(region=ifelse(NAME %in% c("Sacramento","San Joaquin"),11,region)) %>%
  mutate(region=ifelse(NAME %in% c("Stanislaus","Merced"),12,region)) %>%
  mutate(region=ifelse(NAME %in% c("Fresno","Mono","Inyo","Alpine","Madera"),13,region)) %>%
  mutate(region=ifelse(NAME %in% c("Kern","Tulare","Kings"),14,region)) %>%
  mutate(region=ifelse(NAME %in% c("Los Angeles","San Bernardino"),15,region)) %>%
  mutate(region=ifelse(NAME %in% c("San Diego","Orange","Riverside","Imperial"),16,region)) %>%
  mutate(region=ifelse(NAME %in% c("Yolo"),17,region)) 

# Prep production data
reg_1920_prod_all <- rbind(region_prod_2019,region_prod_2020) %>%
  filter(`Type and Variety...1`=="TOTAL ALL VARIETIES") %>%
  select(-`Type and Variety...1`,-`State Total`,-`Prev State Total`) %>%
  pivot_longer(-year)%>%
  pivot_wider(names_from=year,values_from=value)%>%
  mutate(growth=100*(`2020`/`2019`-1)) %>%
  #mutate(growth=(`2020`-`2019`)) %>%
  mutate(region=as.integer(name)) 

# Make map
plot_prod_map <- cb %>%
  left_join(reg_1920_prod_all) %>%
  tm_shape()+ 
    tm_polygons("growth",id="fips", title="% Change", palette="RdYlGn", style="cont", n=8, midpoint=0)+ # see below for other palettes
    tm_legend(position = c("right", "top"), 
            frame = TRUE,
            #legend.outside=T,
            bg.color="white")+
    tm_layout(main.title = "Percent Change in Tons of Grapes Crushed: 2019 to 2020",
            main.title.size=1.2,
            main.title.position = "center") +
    tm_text("NAME", scale=.5) 

# Draw Map
plot_prod_map
tmap_save(plot_prod_map,"CA_prod_map.png")


# Prep returns data
reg_1920_val_all <- rbind(region_val_2019,region_val_2020) %>%
  filter(`Type and Variety...1`=="TOTAL ALL VARIETIES") %>%
  select(-`Type and Variety...1`,-`State Total`,-`Prev State Total`) %>%
  pivot_longer(-year) %>%
  pivot_wider(names_from=year,values_from=value) %>%
  mutate(growth=100*(`2020`/`2019`-1)) %>%
  mutate(region=as.integer(name)) 


# Make map
plot_ret_map <- cb %>%
  left_join(reg_1920_val_all) %>%
  tm_shape()+ 
  tm_polygons("growth",id="fips", title="% Change", palette="RdYlGn", style="cont", n=8, midpoint=0)+ # see below for other palettes
  tm_legend(position = c("right", "top"), 
            frame = TRUE,
            #legend.outside=T,
            bg.color="white")+
  tm_layout(main.title = "Percent Change in Grower Returns per Ton: 2019 to 2020",
            main.title.size=1.2,
            main.title.position = "center") +
  tm_text("NAME", scale=.5) 

# Draw map
plot_ret_map
tmap_save(plot_ret_map,"CA_crushreturns_map.png")



#################################################
#   percent purchased
#################################################

# Create plot
plot_purch <- rbind(region_purch_2019,region_purch_2020) %>%
  filter(`Type and Variety...1`=="TOTAL ALL VARIETIES") %>%
  select(-`Type and Variety...1`,-`State Total`,-`Prev State Total`) %>%
  pivot_longer(-year)%>%
  pivot_wider(names_from=year,values_from=value) %>%
  rename(`2019pu`=`2019`,`2020pu`=`2020`) %>%
  left_join(reg_1920_prod_all) %>%
  mutate(prop_pu_2019=100*`2019pu`/`2019`,prop_pu_2020=100*`2020pu`/`2020`) %>%
  select(region,prop_pu_2019,prop_pu_2020) %>%
  pivot_longer(-region) %>%
  mutate(name=ifelse(name=="prop_pu_2019","2019","2020")) %>%
  ggplot(aes(x=factor(region,levels=seq(1:17)),y=value,fill=name))+
    geom_bar(stat="identity", position = "dodge",width=.5)+
    theme_minimal()+
    ggtitle(paste0("California Percent of Grapes Sold for Crushing"))+
    labs(x="Region",y="Percent",fill="Year")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# draw plot
plot_purch
ggsave("crush_purchases.png")