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/local    ##
#############################################################


# 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,scales,tmap,pals,tigris)

########################################
#   GLOBAL EXPORTS/IMPORTS PRODUCTION
########################################

## Get and clean production and trade data
# Original source:  http://www.fao.org/faostat/en/#data/QV
GET("https://files.asmith.ucdavis.edu/FAOSTAT_prod_data_10-19-2021.csv", write_disk(tf1 <- tempfile(fileext = ".csv")))
df_prod <- read_csv(tf1) %>%
  select(Year,Area,Element,Item,Unit,Value) 

# Original source:  http://www.fao.org/faostat/en/#data/TCL
GET("https://files.asmith.ucdavis.edu/FAOSTAT_trade_data_10-19-2021.csv", write_disk(tf2 <- tempfile(fileext = ".csv")))
df_trade <- read_csv(tf2) %>%
  select(Year,Area,Element,Item,Unit,Value) 


df <- bind_rows(df_prod,df_trade) %>%
  mutate(Item=ifelse(Item=="Vegetables and Fruit Primary","Fruit and Vegetables",Item)) %>%
  mutate(Item=ifelse(Item=="Agricultural Products","Agriculture",Item)) %>%
  mutate(Item=ifelse(Item=="Cereals, Total","Cereals",Item)) %>%
  mutate(Item=ifelse(Item=="Food Excluding Fish","Food",Item)) %>%
  filter(Item %in% c("Food","Fruit and Vegetables","Cereals","Agriculture")) %>%
  pivot_wider(values_from=Value,names_from=Element) %>%
  arrange(Area,Item,Year)


# World plot
plot_world <- df %>%
  filter(Area=="World") %>%
  filter(!Item=="Agriculture") %>%
  rename(Production=`Gross Production Value (current thousand US$)`) %>%
  mutate(Consumption=Production +`Import Value`-`Export Value`) %>%
  filter(Year>1990) %>%
  #mutate(Local=100*(`Import Value`/Consumption)) %>%
  mutate(Local=100*(`Import Value`/Production)) %>%
  ggplot(aes(x=Year,y=Local,color=factor(Item,levels=c("Food","Cereals","Fruit and Vegetables")))) +
    geom_line(size=1.5) +
    labs(color="Item", x = "Year", y = "Percent",caption="US$ value of global imports/US$ value of global production") +
    ggtitle(paste0("Food Imports as Percent of Production (World)"))+
    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
ggsave(paste0("world_local.png"))



# World plot
plot_continent <- df %>%
  filter(Area %in% c("Europe","Asia","Americas","Africa")) %>%
  filter(!Item=="Agriculture") %>%
  rename(Production=`Gross Production Value (current thousand US$)`) %>%
  mutate(Consumption=Production +`Import Value`-`Export Value`) %>%
  filter(Year>1990) %>%
  #mutate(Local=100*(`Export Value`/Consumption)) %>%
  mutate(Local=100*(`Import Value`/Production)) %>%
  ggplot(aes(x=Year,y=Local,color=factor(Item,levels=c("Food","Cereals","Fruit and Vegetables")))) +
    geom_line(size=1.5) +
    facet_wrap(~Area) +
    labs(color="Item", x = "Year", y = "Percent",caption="US$ value of continent imports/ US$ value of continent production") +
    ggtitle(paste0("Food Imports as Percent of 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_continent
ggsave(paste0("continent_local.png"))



# World plot
plot_usa <- df %>%
  filter(Area=="United States of America") %>%
  filter(!Item=="Agriculture") %>%
  rename(Production=`Gross Production Value (current thousand US$)`) %>%
  mutate(Consumption=Production +`Import Value`-`Export Value`) %>%
  filter(Year>1990) %>%
  #mutate(Local=100*(`Import Value`/Consumption)) %>%
  mutate(Local=100*(`Import Value`/Production)) %>%
  ggplot(aes(x=Year,y=Local,color=factor(Item,levels=c("Food","Cereals","Fruit and Vegetables")))) +
    geom_line(size=1.5) +
    labs(color="Item", x = "Year", y = "Percent",caption="US$ value of US imports/US$ value of US production") +
    ggtitle(paste0("Food Imports as Percent of Production (USA)"))+
    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_usa
ggsave(paste0("usa_local.png"))




########################################################
##  NASS data on US production
########################################################

# 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
nassqs_param_values("statisticcat_desc") # get commodity list


# set years for survey data
# you will receive an "bad request - invalid query" error if there are no data for your chosen years.
years <- 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)
params <-list(source_desc="SURVEY",agg_level_desc = c("STATE","NATIONAL"),statisticcat_desc="PRODUCTION",freq_desc="ANNUAL") 

# create national data based on the parameters and years using 'nassqs'
df_survey <- map(years, function(x){params[["year"]] <- x; nassqs(params)}) %>% bind_rows() %>% as.data.table() %>%
  mutate(Value = as.numeric(gsub(",", "", Value))) %>%
  filter(location_desc %in% c("US TOTAL","CALIFORNIA")) %>%
  filter(unit_desc =="$") %>%
  select(reference_period_desc,short_desc,commodity_desc,state_alpha,Value) %>%
  pivot_wider(values_from=Value,names_from=state_alpha)


comms<- c("BEANS, DRY EDIBLE, (EXCL CHICKPEAS) - PRODUCTION, MEASURED IN $"
          ,"FIELD CROP TOTALS, INCL POTATOES & OTHER - PRODUCTION, MEASURED IN $"
          ,"VEGETABLE TOTALS - PRODUCTION, MEASURED IN $"
          ,"MILK - PRODUCTION, MEASURED IN $"
          ,"CATTLE, INCL CALVES - PRODUCTION, MEASURED IN $"
          ,"TURKEYS - PRODUCTION, MEASURED IN $"
          ,"ALMONDS, UTILIZED, SHELLED - PRODUCTION, MEASURED IN $"
          ,"GRAPES, UTILIZED - PRODUCTION, MEASURED IN $"
          ,"APPLES, UTILIZED - PRODUCTION, MEASURED IN $"
          ,"TOMATOES, IN THE OPEN - PRODUCTION, MEASURED IN $"
          ,"STRAWBERRIES - PRODUCTION, MEASURED IN $"
          ,"CHERRIES, SWEET, UTILIZED - PRODUCTION, MEASURED IN $"
          ,"BLUEBERRIES, TAME, UTILIZED - PRODUCTION, MEASURED IN $"
          ,"ONIONS, DRY - PRODUCTION, MEASURED IN $"
          ,"BROCCOLI - PRODUCTION, MEASURED IN $"
          ,"CARROTS - PRODUCTION, MEASURED IN $"
          ,"AVOCADOS, UTILIZED - PRODUCTION, MEASURED IN $"
          ,"PEPPERS, BELL - PRODUCTION, MEASURED IN $"
          ,"RASPBERRIES, UTILIZED - PRODUCTION, MEASURED IN $"
          ,"SPINACH - PRODUCTION, MEASURED IN $"
          ,"LETTUCE, HEAD - PRODUCTION, MEASURED IN $"
          ,"SUGARCANE, SUGAR - PRODUCTION, MEASURED IN $"
          ,"PEANUTS - PRODUCTION, MEASURED IN $"
          ,"RICE - PRODUCTION, MEASURED IN $"
          ,"POTATOES - PRODUCTION, MEASURED IN $"
          ,"WHEAT - PRODUCTION, MEASURED IN $"
          ,"CORN, GRAIN - PRODUCTION, MEASURED IN $"
          ,"SOYBEANS - PRODUCTION, MEASURED IN $"
          ,"EGGS - PRODUCTION, MEASURED IN $")
  
comms2 <- c("Population" 
        ,"GDP"
        ,"Field Crop Totals"
        ,"Vegetable Totals"
        ,"Cattle"
        ,"Eggs"
        ,"Milk" 
        ,"Turkeys" 
        ,"Corn"
        ,"Rice"              
        ,"Soybeans"
        ,"Wheat"  
        ,"Almonds"           
        ,"Apples"            
        ,"Avocados"          
        ,"Beans"             
        ,"Blueberries"      
        ,"Broccoli"          
        ,"Carrots"            
        ,"Cherries"          
        ,"Grapes"            
        ,"Lettuce"           
        ,"Onions"           
        ,"Peanuts"           
        ,"Peppers"                  
        ,"Potatoes"          
        ,"Raspberries"      
        ,"Spinach"           
        ,"Strawberries"      
        ,"Sugarcane"         
        ,"Tomatoes")     


# Make CA production plot 
plot_CA <- df_survey %>%
  filter(short_desc %in% comms) %>%
  select(-short_desc) %>%
  mutate(commodity_desc=str_to_title(commodity_desc)) %>%
  mutate(CA = replace_na(CA, 0)) %>%
  add_row(reference_period_desc=c("Other","Other"),commodity_desc=c("Population","GDP"),US=c(329.5,20.94),CA=c(39.51,3.007)) %>%
  mutate(reference_period_desc=ifelse(commodity_desc %in% c("Field Crop Totals","Vegetable Totals"),"Totals",reference_period_desc)) %>%
  mutate(reference_period_desc=ifelse(commodity_desc %in% c("Cattle","Eggs","Milk","Turkeys" ),"Animal",reference_period_desc)) %>%
  mutate(reference_period_desc=ifelse(commodity_desc %in% c("Corn","Rice","Soybeans","Wheat"),"Grains",reference_period_desc)) %>%
  mutate(proportion=100*CA/US) %>%
  ggplot(aes(y=proportion,x=factor(commodity_desc,levels=comms2),fill=reference_period_desc)) +
    geom_bar(stat="identity",position="stack",show.legend = FALSE) +
    labs(x = "Year",y="Percent") +
    ggtitle(paste0("California Production as a Percent of the Nation")) +
    theme_minimal()+
    ylim(0,100)+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=12),axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# Draw and save plot
plot_CA
ggsave("CA_prod.png")



# set years for census data
# you will receive an "bad request - invalid query" error if there are no data for your chosen years.
years <- c(2002,2007,2012,2017)

# 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="CENSUS",domain_desc="ECONOMIC CLASS",agg_level_desc="NATIONAL") 

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






##############################################################
# Plot Number of Operations, Sales and Land by Economic Class
##############################################################
plot_bar_operations2 <- df_census %>%
  filter(short_desc %in% c("COMMODITY TOTALS - OPERATIONS WITH SALES","COMMODITY TOTALS - SALES, MEASURED IN $","AG LAND, CROPLAND, HARVESTED - ACRES","AG LAND, CROPLAND, (EXCL HARVESTED & PASTURED) - ACRES","AG LAND, PASTURELAND - ACRES","AG LAND, WOODLAND, (EXCL PASTURED) - ACRES","AG LAND, (EXCL CROPLAND & PASTURELAND & WOODLAND) - ACRES")) %>%
  mutate(short_desc=ifelse(short_desc %in% c("AG LAND, CROPLAND, HARVESTED - ACRES","AG LAND, CROPLAND, (EXCL HARVESTED & PASTURED) - ACRES","AG LAND, PASTURELAND - ACRES","AG LAND, WOODLAND, (EXCL PASTURED) - ACRES","AG LAND, (EXCL CROPLAND & PASTURELAND & WOODLAND) - ACRES"),"Land",short_desc)) %>%
  mutate(short_desc=ifelse(short_desc=="COMMODITY TOTALS - OPERATIONS WITH SALES","Farms",short_desc)) %>%
  mutate(short_desc=ifelse(short_desc=="COMMODITY TOTALS - SALES, MEASURED IN $","Sales",short_desc)) %>%
  filter(!short_desc=="Land") %>%
  filter(!domaincat_desc=="ECONOMIC CLASS: (1,000,000 TO 2,499,999 $)") %>%
  filter(!domaincat_desc=="ECONOMIC CLASS: (2,500,000 TO 4,999,999 $)") %>%
  filter(!domaincat_desc=="ECONOMIC CLASS: (5,000,000 OR MORE $)") %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (LESS THAN 1,000 $)","<1",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (1,000 TO 2,499 $)","1-5",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (2,500 TO 4,999 $)","1-5",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (5,000 TO 9,999 $)","5-10",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (10,000 TO 24,999 $)","10-25",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (25,000 TO 49,999 $)","25-50",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (50,000 TO 99,999 $)","50-100",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (100,000 TO 249,999 $)","100-500",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (250,000 TO 499,999 $)","100-500",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (500,000 TO 999,999 $)","500-1,000",domaincat_desc)) %>%
  mutate(domaincat_desc=ifelse(domaincat_desc=="ECONOMIC CLASS: (1,000,000 OR MORE $)",">1,000",domaincat_desc)) %>%
  mutate(Value=as.numeric(gsub(",","",Value))/1000) %>%
  select(year,domaincat_desc,short_desc,Value)%>%
  group_by(year,domaincat_desc,short_desc) %>%
  summarise(Value=sum(Value)) %>%
  ggplot(aes(fill=factor(domaincat_desc,levels=c("<1","1-5","5-10","10-25","25-50","50-100","100-500","500-1,000",">1,000"
  )),y=Value,x=factor(year,levels=c(2017,2012,2007,2002)))) +
    geom_bar(stat="identity",position = position_fill(),width=.9) +
    facet_wrap(~short_desc)+
    #geom_text(size = 3, position = position_stack(vjust = .5))+
    ggtitle(paste0("Distributions of US Farms and Farm Sales"))+
    coord_flip() +
    theme_minimal()+
    scale_fill_brewer(palette="Spectral")+
    labs(x=NULL,y="Percent",fill="Sales ($000)")+
    scale_y_continuous(labels = scales::percent_format())+
    theme(plot.title = element_text(hjust = 0.5,size=16),legend.text = element_text(hjust = 1,size=12), text = element_text(size=14))

# Draw and save plot
plot_bar_operations2
ggsave("farm_size.png")