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

#############################################################
## Code to create plots in Ag Data News article            ##
##  Title: Big Drop in California Rice Acres             ##
##  Link:  https://asmith.ucdavis.edu/news/big-drop-california-rice-acres    ##
#############################################################

# 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,RColorBrewer,data.table,dplyr,purrr,lubridate,tmap,tigris,pals,httr,tidyr,stringr,furrr,readxl)

########################################################################################################
# Read FSA data. Code by Seunghyun Lee

folder_to_save <- "Crop Acreage" # specify directory
dir.create(folder_to_save)

# extract relevant URLs and titles ---------------------------------------------------
url <- "https://www.fsa.usda.gov/news-room/efoia/electronic-reading-room/frequently-requested-information/crop-acreage-data/index"

text <- readLines("https://www.fsa.usda.gov/news-room/efoia/electronic-reading-room/frequently-requested-information/crop-acreage-data/index") %>%
  grep("/Assets/USDA-FSA-Public/usdafiles/NewsRoom/eFOIA/crop-acre-data/zips/", ., value = T) %>%
  str_split("/Assets/USDA-FSA-Public/usdafiles/NewsRoom/eFOIA/crop-acre-data/zips")

df_url_title <- data.frame(text = text[[1]]) %>%
  dplyr::filter(str_detect(text, "\" title")) %>%
  tidyr::separate("text", "\" title", into = c("url", "title")) %>%
  dplyr::filter(str_ends(url, ".zip"))

saveRDS(df_url_title, file = paste0(folder_to_save, "/", "df_url_title.rds"))


# parallel processing  ----------------------------------------------------
plan(multisession)

future_map(df_url_title$url, function(url) {
  
  #--- download zip file using url ---#
  temp <- tempfile()
  endpoint <- paste0(
    "https://www.fsa.usda.gov/Assets/USDA-FSA-Public/usdafiles/NewsRoom/eFOIA/crop-acre-data/zips/",
    url
  )
  
  zipfile <- paste0(folder_to_save, "/", basename(url))
  download.file(endpoint,
                destfile = zipfile
  )
  
  #--- unzip and extract .xlsx file ---#
  unzip(zipfile, exdir = folder_to_save)
  
  #--- remove zipfile ---#
  file.remove(zipfile)
})




###############################################################################

# Read downloaded FSA files
files <- c("2012_fsa_acres_jan_2013.xlsx","2013_fsa_acres_jan_2014.xlsx","2014_fsa_acres_jan2014.xlsx","2015_fsa_acres_01052016.xlsx",
          "2016_fsa_acres_010417.xlsx","2017_fsa_acres_010418.xlsx","2018_fsa_acres_012819.xlsx","2019_fsa_acres_web_010220.xlsx",
          "2020_fsa_acres_web_010521.xlsx")
df_FSA <- map(files, function(x){read_excel(paste0("./Crop Acreage/",x),sheet="county_data") %>% 
                                 mutate(Year = as.numeric(str_sub(x,1,4))) }) %>% 
          bind_rows() 
df_FSA_CA <- df_FSA %>%
  select(-'State County Code') %>%
  bind_rows(read_excel(paste0("./Crop Acreage/","2021_fsa_acres_web_010322.xlsx"),sheet="county_data") %>% mutate_at(c('Planted Acres','Volunteer Acres','Prevented Acres','Failed Acres','Not Planted Acres','Planted and Failed Acres'), as.numeric) %>% mutate(Year = 2021) %>% select(-'State County Code')) %>%
  bind_rows(read_excel(paste0("./Crop Acreage/","2022_fsa_acres_web_082222.xlsx"),sheet="county_data") %>% mutate_at(c('Planted Acres','Volunteer Acres','Prevented Acres','Failed Acres','Not Planted Acres','Planted and Failed Acres'), as.numeric) %>% mutate(Year = 2022) %>% select(-'State County Code')) %>%
  filter(State=="California"|State=="CALIFORNIA") %>%
  mutate(Crop=str_to_title(Crop)) %>%
  mutate(County=str_to_title(County)) %>%
  mutate(Crop=ifelse(str_detect(Crop,"Rice|rice"),"Rice",Crop)) %>%
  mutate(Crop=ifelse(str_detect(Crop,"Cotton|cotton"),"Cotton",Crop)) %>%
  mutate(fips=as.integer(paste0(`State Code`,`County Code`))) %>%
  group_by(Year,fips,County,Crop) %>%
  summarize(Prevented=sum(`Prevented Acres`),Planted=sum(`Planted Acres`)) %>%
  ungroup() 


files_prov <- c("2012_fsa_acres_Aug_2012.xlsx","2013_fsa_acres_aug_2013.xlsx","2014_fsa_acres_aug2014.xlsx","2015_fsa_acres_08032015.xlsx",
           "2016_fsa_acres_ddbk6438.xlsx","2017_fsa_acres_080117.xlsx","2018_fsa_acres_080118.xlsx","2019_fsa_acres_080119.xlsx",
           "2020_fsa_acres_web_073120.xlsx","2021_fsa_acres_web_080221.xlsx")
df_FSA_prov <- map(files_prov, function(x){read_excel(paste0("./Crop Acreage/",x),sheet="county_data") %>% 
    mutate(Year = as.numeric(str_sub(x,1,4))) }) %>% 
  bind_rows() 
df_FSA_CA_prov <- df_FSA_prov %>%
  select(-'State County Code') %>%
  mutate(Year2=as.character(Year)) %>%
  filter(State=="California"|State=="CALIFORNIA"|State=="CALIFRNIA") %>%
  mutate(Crop=str_to_title(Crop)) %>%
  mutate(County=str_to_title(County)) %>%
  mutate(Crop=ifelse(str_detect(Crop,"Rice|rice"),"Rice",Crop)) %>%
  mutate(Crop=ifelse(str_detect(Crop,"Cotton|cotton"),"Cotton",Crop)) %>%
  mutate(fips=as.integer(paste0(`State Code`,`County Code`))) %>%
  group_by(Year,Year2,fips,County,Crop) %>%
  summarize(Prevented=sum(`Prevented Acres`),Planted=sum(`Planted Acres`)) %>%
  ungroup() 




# Read NASS/CDFA acreage data (source: https://asmith.ucdavis.edu/data/california-crops)
df_NASS_CDFA <- read_csv("https://files.asmith.ucdavis.edu/CDFA_vs_NASS.csv") %>%
  filter(CA_variable=="Harvested_Acres") %>%
  select(Year,CA_commodity,value,source) %>%
  rename(Crop=CA_commodity) %>%
  mutate(Crop=ifelse(Crop=="Pr. Tomatoes","Tomatoes",Crop)) %>%
  mutate(Crop=ifelse(str_detect(Crop,"Grapes"),"Grapes",Crop)) %>%
  mutate(Crop=ifelse(str_detect(Crop,"Alfalfa"),"Alfalfa",Crop)) %>%
  group_by(Year,Crop,source) %>%
  summarise(value=sum(value)) %>%
  arrange(Year,source,Crop) %>%
  ungroup() %>%
  filter(Year>2011) %>%
  pivot_wider(values_from = value, names_from=source)


## Read CDL fallow fields data (compiled by Seunghyun Lee)
GET("https://files.asmith.ucdavis.edu/df_fallow.Rds", write_disk(tf1 <- tempfile(fileext = ".Rds")))
df_cdl <- readRDS(tf1) %>%
  arrange(year) %>%
  select(year,county_name,acr_fallow)%>%
  group_by(year) %>%
  summarise(acr_fallow=sum(acr_fallow)) %>%
  filter(year>2011)
  



##########################################
### Make prevented and fallowed bar chart
##########################################
comms <- c("Fallow","Idle","Prevented Cotton","Prevented Rice","Prevented Other")

plot_bar <- df_FSA_CA %>%
  pivot_longer(cols=c(Prevented,Planted)) %>%
  filter(name=="Prevented"|Crop=="Fallow"|Crop=="Idle") %>%
  mutate(Crop=ifelse(Crop=="Cotton","Prevented Cotton",Crop)) %>%
  mutate(Crop=ifelse(Crop=="Rice","Prevented Rice",Crop)) %>%
  mutate(Crop=ifelse(Crop %in% comms,Crop,"Prevented Other")) %>%
  group_by(Year,Crop) %>%
  summarize(value=sum(value)) %>%
  ungroup() %>%
  mutate(value=value/1000) %>%
  ggplot() +
    geom_bar(aes(x=factor(Year),y=value,fill=factor(Crop,levels=comms)),position="stack", stat="identity") +
    theme_minimal()+
    ggtitle(paste0("California Fallow, Idle, and Prevented-Planting Acres"))+
    labs(x="Year",y="Thousand Acres",fill="Crop",caption="Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))


 
# Draw plot
plot_bar
ggsave("fallowidleprev_bar_2022.png",bg="white")




##########################################
### Make prevented and fallowed bar chart -- provisional vs final
##########################################
every_nth = function(n) {
  return(function(x) {x[c(TRUE, rep(FALSE, n - 1))]})
}

comms <- c("Fallow","Idle","Prevented Cotton","Prevented Rice","Prevented Other")

plot_bar_prov <- df_FSA_CA_prov %>%
  bind_rows(df_FSA_CA%>%mutate(Year2=paste0(Year,"(final)")))%>%
  bind_rows(df_FSA_CA%>%mutate(Year2=paste0(Year,"(null)"),Prevented=0,Planted=0))%>%
  mutate(Year2=ifelse(Year==2022,"2022",Year2))%>%
  pivot_longer(cols=c(Prevented,Planted)) %>%
  filter(name=="Prevented"|Crop=="Fallow"|Crop=="Idle") %>%
  mutate(Crop=ifelse(Crop=="Cotton","Prevented Cotton",Crop)) %>%
  mutate(Crop=ifelse(Crop=="Rice","Prevented Rice",Crop)) %>%
  mutate(Crop=ifelse(Crop %in% comms,Crop,"Prevented Other")) %>%
  group_by(Year2,Crop) %>%
  summarize(value=sum(value)) %>%
  ungroup() %>%
  mutate(value=value/1000) %>%
  ggplot() +
  geom_bar(aes(x=Year2,y=value,fill=factor(Crop,levels=comms)),position="stack", stat="identity") +
    theme_minimal()+
    ggtitle(paste0("California Fallow, Idle, and Prevented-Planting Acres: Provisional vs Final"))+
    labs(x="Year",y="Thousand Acres",fill="Crop",caption=" Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n first bar is provisional data (Aug); second bar is final data (Jan)  \n https://agdatanews.substack.com")+
    scale_x_discrete(breaks = every_nth(n = 3))+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14), axis.text.x = element_text(hjust=0.25))



# Draw plot
plot_bar_prov
ggsave("fallowidleprev_bar_2022_prov.png",bg="white")




###################################
### Make planted/prevented bar chart
###################################
comms_cr <- c("Cotton","Rice","Barley","Corn","Wheat")
plot_bar_cr <- df_FSA_CA %>%
  pivot_longer(cols=c(Prevented,Planted)) %>%
  filter(Crop %in% comms_cr) %>%
  group_by(Year,Crop,name) %>%
  summarize(value=sum(value)) %>%
  ungroup() %>%
  mutate(value=value/1000) %>%
  ggplot(aes(x=factor(Year),y=value,fill=factor(name,levels=c("Prevented","Planted")))) +
    geom_bar(position="stack", stat="identity") +
    facet_grid(rows=vars(factor(Crop,levels=comms_cr))) +
    theme_minimal()+
    ggtitle(paste0("California Prevented and Planted Acres"))+
    labs(x="Year",y="Thousand Acres",fill="Status",caption="Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw plot
plot_bar_cr
ggsave("cottonrice_bar_2022.png",bg="white")




###################################
### Make plant acres charts
###################################

# Draw Plot 1
comms_plant1 <- c("Alfalfa","Rice","Tomatoes")
plot_bar_pla1 <- df_FSA_CA %>%
  filter(Crop %in% comms_plant1) %>%
  group_by(Year,Crop) %>%
  summarize(Planted=sum(Planted)) %>%
  ungroup() %>%
  left_join(df_NASS_CDFA) %>%
  rename(`FSA Planted`=Planted) %>%
  rename(`NASS Harvested`=NASS) %>%
  rename(`CDFA Harvested`=CDFA) %>%
  pivot_longer(cols=c(`FSA Planted`,`CDFA Harvested`,`NASS Harvested`)) %>%
  mutate(value=value/1000) %>%
  ggplot(aes(x=Year,y=value,color=name)) +
    geom_line(size=1) +
    facet_grid(cols=vars(Crop),scales="free_y") +
    theme_minimal()+
    scale_x_continuous(breaks=c(2012,2015,2018,2021))+
    ggtitle(paste0("California Crop Acres by Source"))+
    labs(x="Year",y="Thousand Acres",color="Source",caption="Source: https://asmith.ucdavis.edu/data/california-crops \n Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw plot 1
plot_bar_pla1
ggsave("plant1_bar_2022.png",bg="white")


# Make Plot 2
comms_plant2 <- c("Almonds","Pistachios","Grapes")
plot_bar_pla2 <- df_FSA_CA %>%
  filter(Crop %in% comms_plant2) %>%
  group_by(Year,Crop) %>%
  summarize(Planted=sum(Planted)) %>%
  ungroup() %>%
  left_join(df_NASS_CDFA) %>%
  rename(`FSA Planted`=Planted) %>%
  rename(`NASS Harvested`=NASS) %>%
  rename(`CDFA Harvested`=CDFA) %>%
  pivot_longer(cols=c(`FSA Planted`,`CDFA Harvested`,`NASS Harvested`)) %>%
  mutate(value=value/1000) %>%
  ggplot(aes(x=Year,y=value,color=name)) +
    geom_line(size=1) +
    facet_grid(cols=vars(Crop),scales="free_y") +
    theme_minimal()+
    scale_x_continuous(breaks=c(2012,2015,2018,2021))+
    ggtitle(paste0("California Crop Acres by Source"))+
    labs(x="Year",y="Thousand Acres",color="Source",color="Source",caption="Source: https://asmith.ucdavis.edu/data/california-crops \n Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw plot 2
plot_bar_pla2
ggsave("plant2_bar_2022.png",bg="white")


# Make Plot 3
comms_plant3 <- c("Oranges","Lettuce","Strawberries")
plot_bar_pla3 <- df_FSA_CA %>%
  filter(Crop %in% comms_plant3) %>%
  group_by(Year,Crop) %>%
  summarize(Planted=sum(Planted)) %>%
  ungroup() %>%
  left_join(df_NASS_CDFA) %>%
  rename(`FSA Planted`=Planted) %>%
  rename(`NASS Harvested`=NASS) %>%
  rename(`CDFA Harvested`=CDFA) %>%
  pivot_longer(cols=c(`FSA Planted`,`CDFA Harvested`,`NASS Harvested`)) %>%
  mutate(value=value/1000) %>%
  ggplot(aes(x=Year,y=value,color=name)) +
    geom_line(size=1) +
    facet_grid(cols=vars(Crop),scales="free_y") +
    theme_minimal()+
    ggtitle(paste0("California Crop Acres by Source"))+
    labs(x="Year",y="Thousand Acres",color="Source",color="Source",caption="Source: https://asmith.ucdavis.edu/data/california-crops \n Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw plot 3
plot_bar_pla3
ggsave("plant3_bar_2022.png",bg="white")





##########################################
### Make prevented and fallowed bar chart (with CDL)
##########################################
comms <- c("Fallow","Idle","Prevented Cotton","Prevented Rice","Prevented Other")

plot_bar_CDL <- df_FSA_CA %>%
  pivot_longer(cols=c(Prevented,Planted)) %>%
  filter(name=="Prevented"|Crop=="Fallow"|Crop=="Idle") %>%
  mutate(Crop=ifelse(Crop=="Cotton","Prevented Cotton",Crop)) %>%
  mutate(Crop=ifelse(Crop=="Rice","Prevented Rice",Crop)) %>%
  mutate(Crop=ifelse(Crop %in% comms,Crop,"Prevented Other")) %>%
  group_by(Year,Crop) %>%
  summarize(value=sum(value)) %>%
  ungroup() %>%
  mutate(value=value/1000) %>%
  ggplot() +
    geom_bar(aes(x=as.integer(Year),y=value,fill=factor(Crop,levels=comms)),position="stack", stat="identity") +
    geom_line(data=df_cdl,aes(x=as.integer(year),y=as.integer(acr_fallow)/1000),color="black",size=1)+
    geom_text(aes(label="CDL Fallow/Idle",x=2019,y=1700), 
              position=position_nudge(0.1), color="black", size=4, hjust=0.5, show.legend=FALSE) +
    theme_minimal()+
    ggtitle(paste0("Two Sources of California Fallow, Idle, and Prevented-Planting Acres"))+
    labs(x="Year",y="Thousand Acres",fill="FSA Crop")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))



# Draw plot
plot_bar_CDL
ggsave("fallowidleprevCDL_bar.png")






###################################
### Make prevented plant bar chart
###################################
comms_prev <- c("Fallow","Barley","Corn","Wheat","Cotton","Rice","Other")
plot_bar_prev <- df_FSA_CA %>%
  mutate(Crop=ifelse(Crop %in% comms_prev,Crop,"Other"))%>%
  group_by(Year,Crop) %>%
  summarize(Prevented=sum(Prevented)) %>%
  ungroup() %>%
  mutate(Prevented=Prevented/1000) %>% 
  ggplot(aes(x=factor(Year),y=Prevented,fill=factor(Crop,levels=comms_prev))) +
    geom_bar(position="stack", stat="identity") +
    theme_minimal()+
    ggtitle(paste0("California Prevented Planting Acres"))+
    labs(x="Year",y="Thousand Acres",fill="Crop",color="Source",caption="Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com")+
    scale_fill_brewer(palette = "Set2") +
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw plot
plot_bar_prev
ggsave("prevent_bar.png")








###################################
### Make fallow/idle bar chart
###################################

# Make plot
plot_bar_idle <- df_FSA_CA %>%
  filter(Crop %in% c("Fallow","Idle")) %>%
  group_by(Year,Crop) %>%
  summarize(Planted=sum(Planted)) %>%
  ungroup() %>%
  mutate(Planted=Planted/1000) %>%
  ggplot(aes(x=factor(Year),y=Planted,fill=factor(Crop,levels=c("Fallow","Idle")))) +
    geom_bar(position="stack", stat="identity") +
    theme_minimal()+
    ggtitle(paste0("California Fallowed and Idle Acres"))+
    labs(x="Year",y="Thousand Acres",fill="Status")+
    theme(plot.title = element_text(hjust = 0.5,size=16), text = element_text(size=14))

# Draw plot
plot_bar_idle
ggsave("fallow_bar_2022.png",bg="white")






#######
# map #
#######


# get county boundary, merge with df, and prep data
CA_counties <- counties(cb=T) %>% # for state map, use `states(cb=T)`
  filter(STATEFP %in% c("06")) %>%
  mutate(fips = as.integer(GEOID)) 


cb_rice <- df_FSA_CA %>%
  mutate(Prevented=ifelse(is.na(Prevented),0,Prevented/1000)) %>%
  mutate(Planted=ifelse(is.na(Planted),0,Planted/1000)) %>%
  mutate(Percent=100*Prevented/(Prevented+Planted)) %>%  
  filter(Year==2022)%>% 
  filter(Crop=="Rice") %>%
  right_join(CA_counties,by="fips") %>%
  mutate(county_ag = ifelse(Planted+Prevented>0,NAME,NA)) %>%
  st_as_sf()
  


# draw map of 2022 prevented rice planting
cb_rice <- CA_counties %>%
  left_join(df_FSA_CA %>% filter(Year==2022) %>%  filter(Crop=="Rice"),by="fips") %>%
  mutate(Prevented=ifelse(is.na(Prevented),0,Prevented/1000)) %>%
  mutate(Planted=ifelse(is.na(Planted),0,Planted/1000)) %>%
  mutate(Planned=Planted+Prevented) %>%
  mutate(Percent=100*Prevented/(Prevented+Planted)) %>%  
  mutate(county_ag = ifelse(Planted+Prevented>1,NAME,NA)) 

CA_map_prevented_rice <- tm_shape(cb_rice)+ 
  tm_polygons("Prevented",style="fixed",breaks=c(0,1,10,50,100,150),colorNA="white", title="000 Acres", palette =  brewer.greens(10))+ 
  tm_legend(position = c("left", "bottom"), 
            frame = TRUE,
            #  legend.outside=T,
            bg.color="white")+
  tm_layout(main.title = "Prevented Rice Acres in 2022",
            main.title.size=1.2,
            main.title.position = "center")+
  tm_credits("Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com",position = c("RIGHT", "BOTTOM")) +
  tm_text("county_ag", scale=.5) 


CA_map_prevented_rice
tmap_save(CA_map_prevented_rice,"CA_prevented_rice_map.png")



CA_map_planned_rice <- tm_shape(cb_rice)+ 
  tm_polygons("Planned",style="fixed",breaks=c(0,1,10,50,100,150),colorNA="white", title="000 Acres", palette =  brewer.greens(10))+ 
  tm_legend(position = c("left", "bottom"), 
            frame = TRUE,
            #  legend.outside=T,
            bg.color="white")+
  tm_layout(main.title = "Planned Rice Acres in 2022",
            main.title.size=1.2,
            main.title.position = "center")+
  tm_credits("Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com",position = c("RIGHT", "BOTTOM")) +
  tm_text("county_ag", scale=.5) 


CA_map_planned_rice
tmap_save(CA_map_planned_rice,"CA_planned_rice_map.png")



CA_map_percent_prevented_rice <- tm_shape(cb_rice)+ 
  tm_polygons("Percent",style="fixed",breaks=c(0,50,75,90,100),colorNA="white", title="Percent", palette =  brewer.greens(10))+ 
  tm_legend(position = c("left", "bottom"), 
            frame = TRUE,
            #  legend.outside=T,
            bg.color="white")+
  tm_layout(main.title = "Percent of Planned Rice Acres Prevented from Planting in 2022",
            main.title.size=1,
            main.title.position = "center")+
  tm_credits("Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com",position = c("RIGHT", "BOTTOM")) +
  tm_text("county_ag", scale=.5) 


CA_map_percent_prevented_rice
tmap_save(CA_map_percent_prevented_rice,"CA_percent_prevented_rice_map.png")





# draw map of 2022 prevented planting
cb_prevent <- CA_counties %>%
  left_join(df_FSA_CA %>% filter(Year==2022) %>% group_by(fips) %>% summarise(Prevented=sum(Prevented),Planted=sum(Planted),na.rm=T),by="fips") %>%
  mutate(Prevented=ifelse(is.na(Prevented),0,Prevented/1000)) %>%
  mutate(Planted=ifelse(is.na(Planted),0,Planted/1000)) %>%
  mutate(Percent=100*Prevented/(Prevented+Planted)) %>%  
  mutate(county_ag = ifelse(Prevented>1,NAME,NA)) 

CA_map_prevent <- tm_shape(cb_prevent)+ 
  tm_polygons("Prevented",style="fixed",breaks=c(0,10,50,100,150),colorNA="white", title="000 Acres", palette =  brewer.greens(10))+ 
  tm_legend(position = c("left", "bottom"), 
            frame = TRUE,
            #  legend.outside=T,
            bg.color="white")+
  tm_layout(main.title = "Prevented Planting Acres in 2022",
            main.title.size=1.2,
            main.title.position = "center")+
  tm_credits("Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com",position = c("RIGHT", "BOTTOM")) +
  tm_text("county_ag", scale=.5) 


CA_map_prevent
tmap_save(CA_map_prevent,"CA_prevented_map.png")



# draw map of 2022 crop most prevented from planting
cb_prevent_most <- CA_counties %>%
  left_join(df_FSA_CA %>% filter(Year==2022) %>% group_by(fips) %>% summarise(Maxprevented=Crop[which.max(Prevented)],Prevented=sum(Prevented),Planted=sum(Planted),na.rm=T),by="fips") %>%
  mutate(Prevented=ifelse(is.na(Prevented),0,Prevented/1000)) %>%
  mutate(Planted=ifelse(is.na(Planted),0,Planted/1000)) %>%
  mutate(Percent=100*Prevented/(Prevented+Planted)) %>%  
  mutate(Maxprevented=ifelse(Prevented<=1,NA,Maxprevented)) %>%
  mutate(county_ag = ifelse(Prevented>1,NAME,NA)) 

CA_map_prevent_most <- tm_shape(cb_prevent_most)+ 
  tm_polygons("Maxprevented",style="fixed",colorNA="white", title="Crop")+ 
  tm_legend(position = c("left", "bottom"), 
            frame = TRUE,
            #  legend.outside=T,
            bg.color="white")+
  tm_layout(main.title = "Crop Most Prevented from Planting in 2022 (>1000 acres)",
            main.title.size=1.2,
            main.title.position = "center")+
  tm_credits("Source: https://www.fsa.usda.gov/online-services/fsa-online-data-resources/index \n https://agdatanews.substack.com",position = c("RIGHT", "BOTTOM")) +
  tm_text("county_ag", scale=.5) 


CA_map_prevent_most
tmap_save(CA_map_prevent,"CA_prevented_map.png")